Chapter 7. Simulation Methods for VaR for Options and Bonds

x=-3:0.1:3;
plot(x,normcdf(x))
title("CDF of Normal Distribution")

x = -3:0.1:3
using Distributions, Plots;
return plot(x, cdf.(Normal(0,1), x), title = "CDF of Standard Normal", legend = false)

rng default; % set seed
S=10;
rand(S,1)
randn(S,1)
trnd(4,S,1)

using Random;
Random.seed!(12);     # set seed
S = 10;
println(S, " draws from Uniform(0,1): \n", rand(Uniform(0,1), S), "\n") # alternatively, rand(S)
println(S, " draws from Normal(0,1): \n", rand(Normal(0,1), S), "\n")   # alternatively, randn(S)
println(S, " draws from Student-t(4): \n", rand(TDist(4), S), "\n")

yield = [5.00 5.69 6.09 6.38 6.61...
         6.79 6.94 7.07 7.19 7.30]; % yield curve
T = length(yield); 
r=0.07;                             % initial yield rate
Par=10;                             % par value
coupon=r*Par;                       % coupon payments
cc=zeros(1,T)+coupon;               % vector of cash flows 
cc(T)=cc(T)+Par;                    % add par to cash flows
P=sum(cc./((1+yield./100).^(1:T)))  % calculate price

yield_c = [5.00, 5.69, 6.09, 6.38, 6.61,
           6.79, 6.94, 7.07, 7.19, 7.30] # yield curve
T = length(yield_c)
r = 0.07                                 # initial yield rate
Par = 10                                 # par value
coupon = r * Par                         # coupon payments
cc = repeat([coupon], outer = 10)        # vector of cash flows
cc[T] += Par                             # add par to cash flows
P = sum(cc./((1 .+ yield_c/100).^(1:T))) # calc price
println("Price of the bond: ", round(P, digits = 3), " USD")

randn('state',123);         % set the seed
sigma = 1.5;                % daily yield volatility
S = 8;                      % number of simulations
r = randn(1,S)*sigma;       % generate random numbers
ysim=nan(T,S);
for i=1:S
    ysim(:,i)=yield+r(i);
end
ysim=repmat(yield',1,S)+repmat(r,T,1);
plot(ysim)
title("Simulated yield curves")

Random.seed!(12)                # set seed
sigma = 1.5              # daily yield volatility
S = 8                    # number of simulations
r = rand(Normal(0,1), S) # generate random numbers
ysim = fill(NaN, (T,S))
for i in 1:S
    ysim[:,i] = yield_c .+ r[i]
end
plot(ysim, title = "Simulated yield curves")

SP = nan(S,1);
for s = 1:S                                        % S simulations
    SP(s) = sum(cc./(1+ysim(:,s)'./100).^((1:T)));
end
SP = SP-(mean(SP) - P);                            % correct for mean
bar(SP)
S = 50000;
rng("default")
r = randn(S,1) * sigma;
ysim = nan(T,S);
for i = 1:S
    ysim(:,i) = yield' + r(i);
end
SP = nan(S,1);
for i = 1:S
    SP(i) = sum(cc./(1+ysim(:,i)'./100).^((1:T)));
end
SP = SP  - (mean(SP)-P);
histfit(SP)
title("Histogram of simulated bond prices with fitted normal")

using Statistics, StatsPlots;
SP = fill(NaN, S)
for i in 1:S  # S simulations
    SP[i] = sum(cc./(1 .+ ysim[:,i]./100).^(1:T))
end
SP .-= (mean(SP) - P) # correct for mean
bar(SP)
S = 50000
r = randn(S) * sigma
ysim = fill(NaN, (T,S))
for i in 1:S
    ysim[:,i] = yield_c .+ r[i]
end
SP = fill(NaN, S)
for i in 1:S
    SP[i] = sum(cc./(1 .+ ysim[:,i]./100).^(1:T))
end
SP .-= (mean(SP) - P)
histogram(SP,nbins=100,normed=true,xlims=(7,13), legend = false, title = "Histogram of simulated bond prices with fitted normal")
res = fit_mle(Normal, SP)
plot!(Normal(res.μ, res.σ), linewidth = 4, legend = false)

P0 = 50;                      % initial spot price
sigma = 0.2;                  % annual volatility
r = 0.05;                     % annual interest
T = 0.5;                      % time to expiration
X = 40;                       % strike price
f = bs(X,P0,r,sigma,T)        % analytical call price

P0 = 50     # initial spot price
sigma = 0.2 # annual volatility
r = 0.05    # annual interest
T = 0.5     # time to expiration
X = 40      # strike price
f = bs(X = X, P = P0, r = r, sigma = sigma, T = T) # analytical call price

randn('state',0);            % set seed
S = 1e6;                     % number of simulations
ysim = randn(S,1)*sigma*sqrt(T)-0.5*T*sigma^2; % sim returns, lognorm corrected
F = P0*exp(r*T)*exp(ysim);   % sim future prices
SP = F-X;                    % payoff
SP(find(SP < 0)) = 0;        % set negative outcomes to zero
fsim = SP * exp(-r*T) ;      % discount
mean(fsim)                   % simulated price  

Random.seed!(12)        # set seed
S = 10^6                # number of simulations
ysim = randn(S) * sigma * sqrt(T) .- 0.5 * sigma^2 * T # sim returns, lognorm corrected
F = P0 * exp(r * T) * exp.(ysim)            # simulated future prices
SP = F .- X              # payoff
SP[SP.<0] .= 0           # set negative outcomes to zero
fsim = SP * exp(-r * T) # discount
call_sim = mean(fsim)     # simulated price
println("Simulated call price: ", round(call_sim, digits = 3))

subplot(1,2,1)
histfit(F);
title("Simulated prices");
xline(X, 'LineWidth', 1, 'label', 'Strike');
subplot(1,2,2)
hist(fsim,100);
title("Option price density");
xline(mean(fsim), 'LineWidth', 1, 'label', 'Call');

histogram(F, bins = 100, normed = true, xlims = (20,80), title = "Simulated prices density", label = false)
res = fit_mle(Normal, F)
plot!(Normal(res.μ, res.σ), linewidth = 4, label = "Fitted normal")
vline!([X], linewidth = 4, color = "black", label = "Strike price")
histogram(fsim, bins = 110, normed = true, xlims = (0,35), title = "Option density", label = false)
vline!([f["Call"]], linewidth = 4, color = "black", label = "Call price")

randn('state',0);   % set seed
S = 1e7;            % number of simulations
s2 = 0.01^2;        % daily variance
p = 0.01;           % probability
r = 0.05;           % annual riskfree rate
P = 100;            % price today
ysim = randn(S,1)*sqrt(s2)+r/365-0.5*s2; % sim returns
Psim = P*exp(ysim);  % sim future prices 
q = sort(Psim-P);  % simulated P/L
VaR1 = -q(S*p)

Random.seed!(1)                                     # set seed
S = 10^7                                            # number of simulations
s2 = 0.01^2                                         # daily variance
p = 0.01                                            # probability
r = 0.05                                            # annual riskfree rate
P = 100                                             # price today
ysim = randn(S) * sqrt(s2) .+ r/365 .- 0.5 * s2     # sim returns
Psim = P * exp.(ysim)                               # sim future prices
q = sort(Psim .- P)                                 # simulated P/L
VaR1 = -q[ceil(Int, p*S)]
println("Simulated VaR", ceil(Int, p*100), "%: ", round(VaR1, digits = 3), " USD")

T = 0.25;                           % time to expiration
X = 100;                            % strike price
sigma = sqrt(s2*250);               % annual volatility
f = bs(X,P,r,sigma,T);              % analytical call price
fsim=bs(X,Psim,r,sigma,T-(1/365));  % sim option prices
q = sort(fsim.Call-f.Call);         % simulated P/L
VaR2 = -q(p*S)

T = 0.25                            # time to expiration
X = 100                             # strike price
sigma = sqrt(s2 * 250)              # annual volatility
f = bs(X = X,P = P,r = r,sigma = sigma,T = T)               # analytical call price
fsim = bs(X = X,P = Psim, r = r, sigma = sigma, T = T-(1/365)) # sim option prices
q = sort(fsim["Call"] .- f["Call"])  # simulated P/L
VaR2 = -q[ceil(Int, p*S)]
println("Simulated Option VaR", ceil(Int, p*100), "%: ", round(VaR2, digits = 3), " USD")

X1 = 100;
X2 = 110;
f1 = bs(X1,P,r,sigma,T);
f2 = bs(X2,P,r,sigma,T);  
f1sim=bs(X1,Psim,r,sigma,T-(1/365));
f2sim=bs(X2,Psim,r,sigma,T-(1/365));
q = sort(f1sim.Call+f2sim.Put+Psim-f1.Call-f2.Put-P); 
VaR3 = -q(p*S)

X1 = 100
X2 = 110
f1 = bs(X = X1,P = P,r = r,sigma = sigma,T = T)
f2 = bs(X = X2,P = P,r = r,sigma = sigma,T = T)
f2sim = bs(X = X2,P = Psim,r = r,sigma = sigma,T = T-(1/365))
f1sim = bs(X = X1,P = Psim,r = r,sigma = sigma,T = T-(1/365))
q = sort(f1sim["Call"] .+ f2sim["Put"] .+ Psim .- f1["Call"] .- f2["Put"] .- P)
VaR3 = -q[ceil(Int, p*S)]
println("Portfolio Option VaR", ceil(Int, p*100), "%: ", round(VaR3, digits = 3), " USD")

randn('state',12)                 % set seed
mu = [r/365 r/365]';              % return mean
Sigma=[0.01 0.0005; 0.0005 0.02]; % covariance matrix
y = mvnrnd(mu,Sigma,S);           % simulated returns

Random.seed!(12);                          # set seed
mu = Vector([r/365, r/365]);        # return mean
Sigma = [0.01 0.0005; 0.0005 0.02]; # covariance matrix
y = rand(MvNormal(mu,Sigma), S);    # simulated returns

K = 2; 
P = [100 50]';                        % prices
x = [1 1]';                           % number of assets
Port = P'*x;                          % portfolio at t
Psim = repmat(P,1,S)' .*exp(y);       % simulated prices
PortSim=Psim * x;                     % simulated portfolio value
q = sort(PortSim-Port);               % simulated P/L
VaR4 = -q(S*p)

K = 2
P = [100 50]                       # prices
x = [1 1]                          # number of assets
Port = reshape(P * x', 1)[1]       # portfolio at t
Psim = repeat(P, outer = [S,1]).*exp.(y)'   # simulated prices
PortSim = reshape(Psim * x', S)    # simulated portfolio value
q = sort(PortSim .- Port)           # simulated P/L
VaR4 = -q[ceil(Int, p * S)]
println("Two-asset VaR", ceil(Int, p*100), "%: ", round(VaR4, digits = 3), " USD")

f = bs(P(2),P(2),r,sigma,T);
fsim=bs(P(2),Psim(:,2),r,sigma,T-(1/365));
q = sort(fsim.Call+Psim(:,1)-f.Call-P(1)); 
VaR5 = -q(p*S)

f = bs(X = P[2], P = P[2], r = r, sigma = sigma, T = T)
fsim = bs(X = P[2], P = Psim[:,2], r = r, sigma = sigma, T = T-(1/365))
q = sort(fsim["Call"] .+ Psim[:,1] .- f["Call"] .- P[1])
VaR5 = -q[ceil(Int, p * S)]
println("Two-asset with option VaR", ceil(Int, p*100), "%: ", round(VaR5, digits = 3), " USD")