function [call,put,DeltaC,DeltaP] = BlackScholes(S,K,r,sigma,T)
%[call,put,DeltaC,DeltaP] = BlackScholes(S0,K,r,sigma,T)
%European call & put option prices under Black-Scholes model
%
%inputs:
%     S0: current stock price 
%      K: strike 
%      r: annual interest rate for continuous compounding
%  sigma: annual volatility
%      T: time to maturity in years
%outputs:
%   call: price of European call option
%    put: price of European put option
% DeltaC: Delta (deriv. w.r.t. S0) of call price 
% DeltaP: Delta (deriv. w.r.t. S0) of put price
%
%The arguments S,K,r,sigma,T can be one of the following
%   * all scalars
%   * one argument is a vector, the others are scalars
%   * all are vectors of the same size
%
%Example: 
%    r=.05; sigma=.3; S0=10; K=11; T=.5;
%    [call,put] = BlackScholes(S0,K,r,sigma,T)
%% plot option prices vs. S0 for S0=6,6.1,...,16:
%    Sv = 6:.1:16;  [call,put] = BlackScholes(Sv,K,r,sigma,T);
%    plot(Sv,call,Sv,put); 
%    xlabel('initial stock price S0'); legend('call','put')
%% plot option prices vs. K for K=6,6.1,...,16:
%    Kv = 6:.1:16;  [call,put] = BlackScholes(S0,Kv,r,sigma,T);
%    plot(Kv,call,Kv,put); 
%    xlabel('strike K'); legend('call','put')

Phi = @(x) erfc(-x/sqrt(2))/2;   % CDF for std. normal distr.

d1 = (log(S./K) + (r+sigma.^2/2).*T)./(sigma.*sqrt(T));
d2 = d1 - sigma.*sqrt(T);
call = S.*Phi(d1) - K.*exp(-r.*T).*Phi(d2);
put = call + K.*exp(-r.*T) - S;   % use put-call parity
DeltaC = Phi(d1);
DeltaP = Phi(d1) - 1;


