function xxx = lect04varsim
N = 30000; % number of simulations
NC = 3; % number of options
T = 1/252;
X = [1.2 1.0 0.8];
M = [ 60  60  30]/252;
H = [-2    4  -1];
sigma = 0.30;
r     = 0.0;
mu    = 0.0;
H0 = bs(X, r, sigma, M) * H'; % initial portfolio value S0=1
errs = randn(N,1); % normal errors
%errs = trnd(3,N,1)/sqrt(3); % t-distributed errors
ST = exp( mu*T + sigma*sqrt(T)*errs); % final asset price
C  = zeros(N,NC);
for i=1:NC;
    C(:,i) = ST.*bs(X(i)./ST,r,sigma,M(i)-T); % final option values
end
sims = C*H';
SimVaR5 = prctile(sims,5); % Monte Carlo VaR
SimVaR1 = prctile(sims,1); % Monte Carlo VaR
%
%H = bs(X, r, sigma, M) .* H / H0; % H is now portfolio _weights_
% Portfolio delta
H0Delta = bsdelta(X, r, sigma, M) * H'; % initial portfolio delta
mP = H0 + H0 * mu * H0Delta * T;
sP = sqrt( H0^2 * sigma^2 * H0Delta^2 * T);
DVaR5 = mP + norminv(0.05) * sP;
DVaR1 = mP + norminv(0.01) * sP;
%
% Portfolio delta-gamma
H0Gamma = bsgamma(X, r, sigma, M) * H'; % initial portfolio gamma
mPG = mP + .5 * H0^2 * sigma^2 * H0Gamma *T;
sPG = sqrt( sP^2 + .75 * H0^4 * sigma^4 * H0Gamma^2 * T^2 );
DGVaR5 = mPG + norminv(0.05)*sPG;
DGVaR1 = mPG + norminv(0.01)*sPG;
%
disp('   initial   MC5%VaR   MC1%VaR')
disp([H0, SimVaR5, SimVaR1])
disp('   initial    D5%VaR    D1%VaR')
disp([H0, DVaR5, DVaR1])
disp('   initial   DG5%VaR   DG1%VaR')
disp([H0, DGVaR5, DGVaR1])
subplot(2,1,1);
normplot(sims); % normal distribution plot
%b = mP-3*sP:6*sP/100:mP+3*sP;
subplot(2,1,2);
hist(sims,100);

function y = bs(x, r, sig, t)
% The simple Black-Scholes European call option Price
d1 = ( -log(x) + ( r + .5*sig.^2 ) .*t ) ./ sig ./sqrt(t);
d2 = d1 - sig.*sqrt(t);
n1 = normcdf(d1);
n2 = normcdf(d2);
y = n1 - x.*n2.*exp(-r.*t);

function y = bsdelta(x, r, sig, t)
% The simple Black-Scholes European call option Delta
d1 = ( -log(x) + ( r + .5*sig.^2 ) .*t ) ./ sig ./sqrt(t);
y = normcdf(d1);

function y = bsgamma(x, r, sig, t)
% The simple Black-Scholes European call option Price
d1 = ( -log(x) + ( r + .5*sig.^2 ) .*t ) ./ sig ./sqrt(t);
n1 = normpdf(d1);
y = n1 ./ ( sig.*sqrt(t) );