function [Op,De,Ga,Th,Ve,Rh]=BSG(T,S0,K,r,v)
% Octave/MATLAB function to evaluate
% the partial derivatives, or "Greeks"
% for European-style Call and Put options
% using the Black-Scholes formulas.
% INPUTS:                      (Example)
%   T  = time to expiry        (1 year)
%   S0 = asset spot price      ($90)
%   K  = strike price          ($95)
%   r  = riskless APR          (0.02)
%   v  = volatility            (0.15)
% OUTPUTS:
%   Op = [C0; P0]         Option premiums
%   De = [DeltaC; DeltaP]  Deltas
%   Ga = [GammaC; GammaP]  Gammas
%   Th = [ThetaC; ThetaP]  Thetas
%   Ve = [VegaC; VegaP]    Vegas=kappas
%   Rh = [RhoC; RhoP]      rhos
% EXAMPLE:
%  [Op,De,Ga,Th,Ve,Rh]=BSG(1,90,95,0.02,0.15)
%
  normcdf = @(x) 0.5*(1.0+erf(x/sqrt(2.0)));
  normpdf = @(x) exp(-x^2/2)/sqrt(2*pi);
  d1 = (log(S0/K)+T*(r+v^2/2)) / (v*sqrt(T));
  d2 = (log(S0/K)+T*(r-v^2/2)) / (v*sqrt(T));
  C0 = S0*normcdf(d1)-K*exp(-r*T)*normcdf(d2);
  P0 = K*exp(-r*T)*normcdf(-d2)-S0*normcdf(-d1);
  DeltaC = normcdf(d1); DeltaP = normcdf(d1)-1;
  GammaP = GammaC = normpdf(d1)/(v*S0*sqrt(T));
  boTh = -v*S0*normpdf(d1)/(2*sqrt(T)); % common
  ThetaC = boTh - r*exp(-r*T)*K*normcdf(d2);
  ThetaP = boTh + r*exp(-r*T)*K*normcdf(-d2);
  VegaP = VegaC = S0*normpdf(d1)*sqrt(T);
  RhoC =  T*exp(-r*T)*K*normcdf(d2);
  RhoP = -T*exp(-r*T)*K*normcdf(-d2);
  Op = [C0; P0];  De = [DeltaC; DeltaP];
  Ga = [GammaC; GammaP]; Th = [ThetaC; ThetaP];
  Ve = [VegaC; VegaP];  Rh = [RhoC; RhoP];
  return
end


