Chapter 7. Simulation Methods for VaR for Options and Bonds


Copyright 2016 Jon Danielsson. Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at. http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.

Listing 7.1: Transformation in R
Last edited: 2011

x=seq(-3,3,by=0.1)
plot(x,pnorm(x),type="l")
		
Listing 7.2: Transformation in Matlab
Last edited: 2011

x=-3:0.1:3;
plot(x,normcdf(x))
		

Listing 7.3: Various RNs in R
Last edited: 2011

set.seed(12)     
S=10
runif(S)
rnorm(S)
rt(S,4)
		
Listing 7.4: Various RNs in Matlab
Last edited: December 2016

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

Listing 7.5: Price bond in R
Last edited: Augua 2016

yield=c(5.00, 5.69, 6.09, 6.38, 6.61, 6.79, 6.94, 7.07, 7.19, 7.30)                  
TT=length(yield)    
r=0.07                          
Par=10                          
coupon=r*Par                    
cc=1:10*0+coupon                
cc[10]=cc[10]+Par                
P=sum(cc/((1+yield/100)^(1:TT)))  
		
Listing 7.6: Price bond in Matlab
Last edited: 2011

yield = [5.00 5.69 6.09 6.38 6.61 6.79 6.94 7.07 7.19 7.30];                             
T = length(yield);           
r=0.07;                              
Par=10;                              
coupon=r*Par;                        
cc=zeros(1,T)+coupon;                
cc(10)=cc(10)+Par;                  
P=sum(cc./((1+yield./100).^(1:T)))  
		

Listing 7.7: Simulate yields in R
Last edited: August 2016

set.seed(12)           
sigma = 1.5            
S = 8                  
r = rnorm(S,0,sigma)   
ysim = matrix(nrow=TT,ncol=S)
for (i in 1:S) ysim[,i]=yield+r[i]
matplot(ysim,type='l')
		
Listing 7.8: Simulate yields in Matlab
Last edited: 2011

S = 8                   
sigma = 1.5             
randn('state',123);      
r = randn(1,S)*sigma;    
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.9: Simulate bond prices in R
Last edited: August 2016

SP = vector(length=S)
for (i in 1:S){                
  SP[i] = sum(cc/((1+ysim[,i]/100)^(TT)))
}
SP = SP-(mean(SP) - P)         
barplot(SP)                    
hist(SP,probability=TRUE)  
x=seq(6,16,length=100)
lines(x, dnorm(x, mean = mean(SP), sd = sd(SP)))
		
Listing 7.10: Simulate bond prices in Matlab
Last edited: 2011

SP = nan(S,1);                      
for s = 1:S                          
  SP(s) = sum(cc./(1+ysim(:,s)'./100).^((1:T))) ;
end
SP = SP-(mean(SP) - P);             
bar(SP)                              
histfit(SP)                          
		

Listing 7.11: Black-Scholes valuation in R
Last edited: August 2016

P0 = 50                  
sigma = 0.2              
r = 0.05                
TT = 0.5                  
X = 40                  
f = bs(X,P0,r,sigma,TT)  
		
Listing 7.12: Black-Scholes valuation in Matlab
Last edited: 2011

P0 = 50;                   
sigma = 0.2;               
r = 0.05;                 
T = 0.5;                   
X = 40;                   
f = bs(X,P0,r,sigma,T);    
		

Listing 7.13: Black-Scholes simulation in R, continued
Last edited: August 2016

S = 1e6                    
set.seed(12)              
F = P0*exp(r*TT)            
ysim = rnorm(S,-0.5*sigma*sigma*TT,sigma*sqrt(TT)) 
F=F*exp(ysim)              
SP = F-X                  
SP[SP<0] = 0               
fsim = SP*exp(-r*TT)       
		
Listing 7.14: Black-Scholes simulation in Matlab, continued
Last edited: 2011

randn('state',0);         
S = 1e6;                   
F = P0*exp(r*T);           
ysim = randn(S,1)*sigma*sqrt(T)-0.5*T*sigma^2; 
F=F*exp(ysim);
SP = F-X;                 
SP(find(SP < 0)) = 0;     
fsim =SP*exp(-r*T) ;      
		

Listing 7.15: Option density plots in R, continued
Last edited: 2011

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)  
		
Listing 7.16: Option density plots in Matlab, continued
Last edited: 2011

histfit(F);      
hist(fsim,100);  
		

Listing 7.17: Simulate VaR in R
Last edited: 2011

set.seed(1)                
S = 1e7                    
s2 = 0.01^2                
p = 0.01                  
r = 0.05                  
P = 100                    
ysim = rnorm(S,r/365-0.5*s2,sqrt(s2))  
Psim = P*exp(ysim)        
q = sort(Psim-P)          
VaR1 = -q[p*S]            
		
Listing 7.18: Simulate VaR in Matlab
Last edited: 2011

randn('state',0);         
S = 1e7;                   
s2 = 0.01^2;               
p = 0.01;                 
r = 0.05;                 
P = 100;                   
ysim = randn(S,1)*sqrt(s2)+r/365-0.5*s2; 
Psim = P*exp(ysim);       
q = sort(Psim-P);         
VaR1 = -q(S*p)
		

Listing 7.19: Simulate option VaR in R
Last edited: August 2016

TT = 0.25;                   
X = 100;                     
sigma = sqrt(s2*250);        
f = bs(X,P,r,sigma,TT)        
fsim = bs(X,Psim,r,sigma,TT-(1/365))  
q = sort(fsim$Call-f$Call)  
VaR2 = -q[p*S]              
		
Listing 7.20: Simulate option VaR in Matlab
Last edited: 2011

T = 0.25;                           
X = 100;                             
sigma = sqrt(s2*250);                 
f = bs(X,P,r,sigma,T);              
fsim=bs(X,Psim,r,sigma,T-(1/365));  
q = sort(fsim.Call-f.Call);         
VaR2 = -q(p*S)                      
		

Listing 7.21: 2 Options
Last edited: August 2016

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]
		
Listing 7.22: 2 Options
Last edited: 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.23: Simulated two asset returns in R
Last edited: 2011

library (MASS)
mu = c(r/365,r/365)             
Sigma = matrix(c(0.01, 0.0005, 0.0005, 0.02),ncol=2)   
set.seed(12)                     
y = mvrnorm(S,mu,Sigma)         
		
Listing 7.24: Simulated two asset returns in Matlab
Last edited: 2011

mu = [r/365 r/365]';                   
Sigma = [0.01 0.0005; 0.0005 0.02];   
randn('state',12)                     
y = mvnrnd(mu,Sigma,S);               
		

Listing 7.25: Two asset VaR in R
Last edited: 2011

K=2
P = c(100,50)                         
x = c(1,1)                             
Port = P %*% x                         
Psim = matrix(t(matrix(P,K,S)),ncol=K)*exp(y)
PortSim = Psim %*% x                  
q = sort(PortSim-Port[1,1])            
VaR4 = -q[S*p]
		
Listing 7.26: Two asset VaR Matlab
Last edited: 2011

K = 2;                                   
P = [100 50]';                           
x = [1 1]';                             
Port = P'*x;                             
Psim = repmat(P,1,S)' .*exp(y);         
PortSim=Psim * x;                        
q = sort(PortSim-Port);                  
VaR4 = -q(S*p)
		

Listing 7.27: A 2 asset case in R with an option
Last edited: August 2016

f = bs(P[2],P[2],r,sigma,TT)   
fsim = bs(P[2],Psim[,2],r,sigma,TT-(1/365))
q = sort(fsim$Call+Psim[,1]-f$Call-P[1]); 
VaR5 = -q[p*S]
		
Listing 7.28: A 2 asset case in Matlab with an option
Last edited: 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)