Chapter 7. Simulation Methods for VaR for Options and Bonds (in MATLAB/Julia)


Copyright 2011 - 2019 Jon Danielsson. This code is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This code is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. The GNU General Public License is available at: https://www.gnu.org/licenses/.


Listing 7.1/7.2: Transformation in MATLAB
Last updated 2011

x=-3:0.1:3;
plot(x,normcdf(x))
		
Listing 7.1/7.2: Transformation in Julia
Last updated June 2018

x = collect(linspace(-3, 3, 61))
using Plots;
using Distributions;
return plot(x, cdf.(Normal(0,1), x))
		

Listing 7.3/7.4: Various RNs in MATLAB
Last updated August 2016

rng default; % set seed
S=10;
rand(S,1)
randn(S,1)
trnd(4,S,1)
		
Listing 7.3/7.4: Various RNs in Julia
Last updated June 2018

srand(12);                     # set seed
S = 10;
println(rand(Uniform(0,1), S)) # alternatively, rand(S)
println(rand(Normal(0,1), S))  # alternatively, randn(S)
println(rand(TDist(4), S))
		

Listing 7.5/7.6: Price bond in MATLAB
Last updated 2011

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(10)=cc(10)+Par;                   % add par to cash flows
P=sum(cc./((1+yield./100).^(1:T)))   % calculate price
		
Listing 7.5/7.6: Price bond in Julia
Last updated June 2018

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[10] += Par                            # add par to cash flows
P = sum(cc./((1+yield_c/100).^(1:T)))    # calc price
		

Listing 7.7/7.8: Simulate yields in MATLAB
Last updated 2011

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)
		
Listing 7.7/7.8: Simulate yields in Julia
Last updated June 2018

srand(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!(Array{Float64}(T,S), NaN)
for i in range(1, S)
    ysim[:,i] = yield_c + r[i]
end
using Plots;
plot(ysim)
		

Listing 7.9/7.10: Simulate bond prices in MATLAB
Last updated 2011

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)
		
Listing 7.9/7.10: Simulate bond prices in Julia
Last updated June 2018

SP = fill!(Array{Float64}(S,1), NaN)
for i in range(1,S)                              # S simulations
    SP[i] = sum(cc./(1+ysim[:,i]./100).^T)
end
SP -= (mean(SP) - P)                             # correct for mean
using Plots;
bar(SP)
S = 50000
r = randn(S) * sigma
ysim = fill!(Array{Float64}(T,S), NaN)
for i in range(1, S)
    ysim[:,i] = yield_c + r[i]
end
SP = fill!(Array{Float64}(S,1), NaN)
for i in range(1,S)
    SP[i] = sum(cc./(1+ysim[:,i]./100).^T)
end
SP -= (mean(SP) - P)
using Plots;
histogram(SP,nbins=100,normed=true,xlims=(7,13))
res = fit_mle(Normal, SP)
plot!(Normal(res.μ, res.σ), linewidth = 4)
		

Listing 7.11/7.12: Black-Scholes valuation in MATLAB
Last updated 2011

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
%% this calculation uses the Black-Scholes pricing function (Listing 6.1/6.2)
		
Listing 7.11/7.12: Simulate bond prices in Julia
Last updated June 2018

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
## this calculation uses the Black-Scholes pricing function (Listing 6.1/6.2)
		

Listing 7.13/7.14: Black-Scholes simulation in MATLAB
Last updated 2011

randn('state',0);                            % set seed
S = 1e6;                                     % number of simulations
F = P0*exp(r*T);                             % futures price
ysim=randn(S,1)*sigma*sqrt(T)-0.5*T*sigma^2; % sim returns, lognorm corrected
F=F*exp(ysim);                               % sim futures price
SP = F-X;                                    % payoff
SP(find(SP < 0)) = 0;                        % set negative outcomes to zero
fsim =SP*exp(-r*T) ;                         % discount
mean(fsim)                                   % simulated price
		
Listing 7.13/7.14: Black-Scholes simulation in Julia
Last updated June 2018

srand(12)                                   # set seed
S = 10^6                                    # number of simulations
F = P0 * exp(r * T)                         # futures price
ysim = randn(S)*sigma*sqrt(T)-0.5*sigma^2*T # sim returns, lognorm corrected
F = F * exp.(ysim)                          # sim futures price
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
		

Listing 7.15/7.16: Option density plots in MATLAB
Last updated 2011

subplot(1,2,1)
histfit(F);
subplot(1,2,2)
hist(fsim,100);
		
Listing 7.15/7.16: Option density plots in Julia
Last updated June 2018

using Plots;
histogram(F, bins = 100, normed = true, xlims = (20,80))
res = fit_mle(Normal, F)
plot!(Normal(res.μ, res.σ), linewidth = 4)
vline!([X], linewidth = 4, color = "black")
histogram(fsim, bins = 110, normed = true, xlims = (0,35))
vline!([f["Call"]], linewidth = 4, color = "black")
		

Listing 7.17/7.18: Simulate VaR in MATLAB
Last updated 2011

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)
		
Listing 7.17/7.18: Simulate VaR in Julia
Last updated June 2018

srand(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[convert(Int, p*S)]
		

Listing 7.19/7.20: Simulate option VaR in MATLAB
Last updated 2011

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)
		
Listing 7.19/7.20: Simulate option VaR in Julia
Last updated June 2018

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[convert(Int, p*S)]
		

Listing 7.21/7.22: Example 7.3 in MATLAB
Last updated 2011

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)
		
Listing 7.21/7.22: Example 7.3 in Julia
Last updated June 2018

X1 = 100
X2 = 110
f1 = bs(X1,P,r,sigma,T)
f2 = bs(X2,P,r,sigma,T)
f2sim = bs(X2,Psim,r,sigma,T-(1/365))
f1sim = bs(X1,Psim,r,sigma,T-(1/365))
q = sort(f1sim["Call"] + f2sim["Put"] + Psim - f1["Call"] - f2["Put"] - P)
VaR2 = -q[convert(Int, p*S)]
		

Listing 7.23/7.24: Simulated two-asset returns in MATLAB
Last updated 2011

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
		
Listing 7.23/7.24: Simulated two-asset returns in Julia
Last updated June 2018

using Distributions;
srand(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
		

Listing 7.25/7.26: Two-asset VaR in MATLAB
Last updated 2011

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)
		
Listing 7.25/7.26: Two-asset VaR in Julia
Last updated June 2018

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

Listing 7.27/7.28: A two-asset case in MATLAB with an option
Last updated 2011

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)
		
Listing 7.27/7.28: A two-asset case in Julia with an option
Last updated June 2018

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[convert(Int, p * S)]