R and Matlab Chapter 7. Simulation Methods for VaR for Options and Bonds

# Chapter 7. Simulation Methods for VaR for Options and Bonds

### R and Matlab

Copyright 2011 - 2023 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: www.gnu.org/licenses.

##### Transformation in R
x = seq(-3, 3, by = 0.1)
plot(x, pnorm(x), type = "l")

##### % Transformation in MATLAB
x=-3:0.1:3;
plot(x,normcdf(x))
title("CDF of Normal Distribution")


##### Various RNs in R
set.seed(12) # set seed
S = 10
runif(S)
rnorm(S)
rt(S,4)

##### % Various RNs in MATLAB
rng default; % set seed
S=10;
rand(S,1)
randn(S,1)
trnd(4,S,1)


##### Price bond in R
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)                       # number of time periods
r = 0.07                                # initial yield rate
Par = 10                                # par value
coupon = r * Par                        # coupon payments
cc = rep(coupon, T)                     # vector of cash flows
cc[T] = cc[T] + Par                     # add par to cash flows
P = sum(cc/((1+yield/100)^(1:T)))      # calculate price
print(P)

##### % Price bond in MATLAB
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


##### Simulate yields in R
set.seed(12)                           # set seed
sigma = 1.5                            # daily yield volatiltiy
S = 8                                  # number of simulations
r = rnorm(S, 0, sigma)                 # generate random numbers
ysim = matrix(nrow=length(yield),ncol=S)
for (i in 1:S) ysim[,i]=yield+r[i]
matplot(ysim,type='l')

##### % Simulate yields in MATLAB
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")


##### Simulate bond prices in R
SP = rep(NA, length = S)
for (i in 1:S){                            # S simulations
SP[i] = sum(cc/((1+ysim[,i]/100)^(1:T)))
}
SP = SP-(mean(SP) - P)                     # correct for mean
par(mfrow=c(1,2), pty="s")
barplot(SP)
hist(SP,probability=TRUE)
x=seq(6,16,length=100)
lines(x, dnorm(x, mean = mean(SP), sd = sd(SP)))
S = 50000
r = rnorm(S, 0, sigma)                 # generate random numbers
ysim = matrix(nrow=length(yield),ncol=S)
for (i in 1:S) ysim[,i]=yield+r[i]
SP = rep(NA, length = S)
for (i in 1:S){                            # S simulations
SP[i] = sum(cc/((1+ysim[,i]/100)^(1:T)))
}
SP = SP-(mean(SP) - P)                     # correct for mean
par(mfrow=c(1,2), pty="s")
barplot(SP)
hist(SP,probability=TRUE)
x=seq(6,16,length=100)
lines(x, dnorm(x, mean = mean(SP), sd = sd(SP)))

##### % Simulate bond prices in MATLAB
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")


##### Black-Scholes valuation in R
P0 = 50                          # initial spot price
sigma = 0.2                      # annual volatility
r = 0.05                         # annual interest
TT = 0.5                         # time to expiration
X = 40                           # strike price
f = bs(X = X, P = P0, r = r, sigma = sigma, T = TT)          # analytical call price
print(f)

##### % Black-Scholes valuation in MATLAB
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


##### Black-Scholes simulation in R
set.seed(12)                                   # set seed
S = 1e6                                        # number of simulations
F = P0*exp(r*TT)                               # futures price
ysim = rnorm(S,-0.5*sigma^2*TT,sigma*sqrt(TT)) # 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*TT)                           # discount
call_sim = mean(fsim)                          # simulated price
print(call_sim)

##### % Black-Scholes simulation in MATLAB
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


##### Option density plots in R
par(mfrow=c(1,2), pty="s")
hist(F, probability=TRUE, ylim=c(0,0.06))
x = seq(min(F), max(F), length=100)
lines(x, dnorm(x, mean = mean(F), sd = sd(SP)))
hist(fsim, nclass=100, probability=TRUE)

##### % Option density plots in MATLAB
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');


##### Simulate VaR in R
set.seed(1)                            # 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 = rnorm(S,r/365-0.5*s2,sqrt(s2))  # sim returns
Psim = P*exp(ysim)                     # sim future prices
q = sort(Psim-P)                       # simulated P/L
VaR1 = -q[p*S]
print(VaR1)

##### % Simulate VaR in MATLAB
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)


##### Simulate option VaR in R
TT = 0.25                                         # time to expiration
X = 100                                           # strike price
sigma = sqrt(s2*250)                              # annual volatility
f = bs(X, P, r, sigma, TT)                        # analytical call price
fsim = bs(X,Psim,r,sigma,TT-(1/365))              # sim option prices
q = sort(fsim$Call-f$Call)                        # simulated P/L
VaR2 = -q[p*S]
print(VaR2)

##### % Simulate option VaR in MATLAB
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)


##### Example 7.3 in R
X1 = 100
X2 = 110
f1 = bs(X1, P, r, sigma, TT)
f2 = bs(X2, P, r, sigma, TT)
f2sim = bs(X2, Psim, r, sigma, TT-(1/365))
f1sim = bs(X1, Psim, r, sigma, TT-(1/365))
q = sort(f1sim$Call + f2sim$Put + Psim-f1$Call - f2$Put-P);
VaR3 = -q[p*S]
print(VaR3)

##### % Example 7.3 in MATLAB
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)


##### Simulated two-asset returns in R
library (MASS)
set.seed(12)  # set seed
mu = rep(r/365, 2)                                    # return mean
Sigma = matrix(c(0.01, 0.0005, 0.0005, 0.02),ncol=2)  # covariance matrix
y = mvrnorm(S,mu,Sigma)  # simulated returns

##### % Simulated two-asset returns in MATLAB
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


##### Two-asset VaR in R
K=2
P = c(100, 50)                                 # prices
x = rep(1, 2)                                  # number of assets
Port = P %*% x                                 # portfolio at t
Psim = matrix(t(matrix(P,K,S)),ncol=K)*exp(y)  # simulated prices
PortSim = Psim %*% x                           # simulated portfolio value
q = sort(PortSim-Port[1,1])                    # simulated P/L
VaR4 = -q[S*p]
print(VaR4)

##### % Two-asset VaR in MATLAB
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)


##### A two-asset case in R with an option
f = bs(X = P[2], P = P[2], r = r, sigma = sigma, T = TT)
fsim = bs(X = P[2], P = Psim[,2], r = r, sigma = sigma, T = TT-(1/365))
q = sort(fsim$Call + Psim[,1] - f$Call - P[1]);
VaR5 = -q[p*S]
print(VaR5)

##### % A two-asset case in MATLAB with an option
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)


##### Financial Risk Forecasting
Market risk forecasting with R, Julia, Python and Matlab. Code, lecture slides, implementation notes, seminar assignments and questions.