 # Chapter 7. Simulation Methods for VaR for Options and Bonds (in R/Python)

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 R Last updated 2011

x=seq(-3,3,by=0.1)
plot(x,pnorm(x),type="l")

##### Listing 7.1/7.2: Transformation in Python Last updated June 2018

import numpy as np
import matplotlib.pyplot as plt
x = np.arange(-3,3.1, step = 0.1) # Python's arange excludes the last value
plt.plot(x, stats.norm.cdf(x))
plt.show()
plt.close()


##### Listing 7.3/7.4: Various RNs in R Last updated 2011

set.seed(12) # set seed
S=10
runif(S)
rnorm(S)
rt(S,4)

##### Listing 7.3/7.4: Various RNs in Python Last updated June 2018

import numpy as np
np.random.seed(12)                        # set seed
S = 10
print (np.random.uniform(size=S))
print (np.random.normal(size=S))
print (np.random.standard_t(df=4,size=S))


##### Listing 7.5/7.6: Price bond in R Last updated August 2019

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

##### Listing 7.5/7.6: Price bond in Python Last updated June 2018

import numpy as np
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 = len(yield_c)
r = 0.07                                          # initial yield rate
Par = 10                                          # par value
coupon = r * Par                                  # coupon payments
cc = [coupon] * 10                                # vector of cash flows
cc += Par                                      # add par to cash flows
P=np.sum(cc/(np.power((1+np.divide(yield_c,100)),
list(range(1,T+1)))))       # calc price
print(P)


##### Listing 7.7/7.8: Simulate yields in R Last updated August 2019

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')

##### Listing 7.7/7.8: Simulate yields in Python Last updated June 2018

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(12)                   # set seed
sigma = 1.5                          # daily yield volatility
S = 8                                # number of simulations
r = np.random.normal(0,sigma,size=S) # generate random numbers
ysim = np.zeros([T,S])
for i in range(S):
ysim[:,i] = yield_c + r[i]
plt.plot(ysim)
plt.show()
plt.close()


##### Listing 7.9/7.10: Simulate bond prices in R Last updated August 2019

SP = vector(length=S)
for (i in 1:S){                                      # S simulations
SP[i] = sum(cc/((1+ysim[,i]/100)^(length(yield))))
}
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)))

##### Listing 7.9/7.10: Simulate bond prices in Python Last updated June 2018

import numpy as np
import matplotlib.pyplot as plt
S = 8
SP = np.zeros([S])
for i in range(S):                                        # S simulations
SP[i] = np.sum(cc/((1+ysim[:,i]/100)**T))
SP -= (np.mean(SP) - P)                                   # correct for mean
plt.bar(range(1,S+1), SP)
plt.show()
plt.close()
S = 50000
r = np.random.normal(0, sigma, size = S)
ysim = np.zeros([T,S])
for i in range(S):
ysim[:,i] = yield_c + r[i]
SP = np.zeros([S])
for i in range(S):
SP[i] = np.sum(cc/((1+ysim[:,i]/100)**T))
SP -= (np.mean(SP) - P)
plt.hist(SP, bins = 30, range = (7, 13), density = True)
fitted_norm=stats.norm.pdf(np.linspace(7,13,30),
np.mean(SP),np.std(SP,ddof=1))
plt.plot(np.linspace(7,13,30), fitted_norm)
plt.show()
plt.close()


##### Listing 7.11/7.12: Black-Scholes valuation in R Last updated August 2019

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,P0,r,sigma,length(yield)) # analytical call price
## this calculation uses the Black-Scholes pricing function (Listing 6.1/6.2)
print(f)

##### Listing 7.11/7.12: Black-Scholes valuation in Python 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)
print (f)


##### Listing 7.13/7.14: Black-Scholes simulation in R Last updated August 2016

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

##### Listing 7.13/7.14: Black-Scholes simulation in Python Last updated June 2018

import numpy as np
np.random.seed(12)                             # set seed
S = 10**6                                      # number of simulations
F = P0 * np.exp(r * T)                         # futures price
ysim=np.random.normal(-0.5*sigma**2*T,
sigma*np.sqrt(T),size=S) # sim returns, lognorm corrected
F = F * np.exp(ysim)                           # sim futures price
SP = F - X                                     # payoff
SP[SP < 0] = 0                                 # set negative outcomes to zero
fsim = SP * np.exp(-r * T)                     # discount
call_sim = np.mean(fsim)                       # simulated price
print(call_sim)


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

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)

##### Listing 7.15/7.16: Option density plots in Python Last updated June 2018

import numpy as np
import matplotlib.pyplot as plt
plt.hist(F, bins = 60, range = (20,80), density = True)
fitted_norm=stats.norm.pdf(np.linspace(20,80,60),np.mean(F),np.std(F,ddof=1))
plt.plot(np.linspace(20,80,60), fitted_norm)
plt.axvline(x=X, color='k')
plt.show()
plt.close()
plt.hist(fsim, bins = 60, range = (0, 35), density = True)
plt.axvline(x=f['Call'], color='k')
plt.show()
plt.close()


##### Listing 7.17/7.18: Simulate VaR in R Last updated 2011

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)

##### Listing 7.17/7.18: Simulate VaR in Python Last updated June 2018

import numpy as np
np.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=np.random.normal(r/365-0.5*s2,np.sqrt(s2),size=S) # sim returns
Psim = P * np.exp(ysim)                                # sim future prices
q = np.sort(Psim - P)                                  # simulated P/L
VaR1 = -q[int(p*S) - 1]
print(VaR1)


##### Listing 7.19/7.20: Simulate option VaR in R Last updated August 2016

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

##### Listing 7.19/7.20: Simulate option VaR in Python Last updated June 2018

import numpy as np
T = 0.25                                # time to expiration
X = 100                                 # strike price
sigma = np.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 = np.sort(fsim['Call']-f['Call'])     # simulated P/L
VaR2 = -q[int(p*S) - 1]
print(VaR2)


##### Listing 7.21/7.22: Example 7.3 in R Last updated 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]
print(VaR3)

##### Listing 7.21/7.22: Example 7.3 in Python Last updated June 2018

import numpy as np
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 = np.sort(f1sim['Call'] + f2sim['Put'] + Psim - f1['Call'] - f2['Put'] - P)
VaR3 = -q[int(p*S) - 1]
print(VaR3)


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

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

##### Listing 7.23/7.24: Simulated two-asset returns in Python Last updated June 2018

import numpy as np
np.random.seed(12)                                     # set seed
mu = np.transpose([r/365, r/365])                      # return mean
Sigma = np.matrix([[0.01, 0.0005],[0.0005, 0.02]])     # covariance matrix
y = np.random.multivariate_normal(mu, Sigma, size = S) # simulated returns


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

K=2
P = c(100,50)                                 # prices
x = c(1,1)                                    # 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)

##### Listing 7.25/7.26: Two-asset VaR in Python Last updated June 2018

import numpy as np
P = np.asarray([100, 50])              # prices
x = np.asarray([1, 1])                 # number of assets
Port = np.matmul(P, x)                 # portfolio at t
Psim=np.matlib.repmat(P,S,1)*np.exp(y) # simulated prices
PortSim = np.matmul(Psim, x)           # simulated portfolio value
q = np.sort(PortSim - Port)            # simulated P/L
VaR4 = -q[int(p*S) - 1]
print(VaR4)


##### Listing 7.27/7.28: A two-asset case in R with an option Last updated August 2016

f = bs(P,P,r,sigma,TT)
fsim = bs(P,Psim[,2],r,sigma,TT-(1/365))
q = sort(fsim$Call+Psim[,1]-f$Call-P);
VaR5 = -q[p*S]
print(VaR5)

##### Listing 7.27/7.28: A two-asset case in Python with an option Last updated June 2018

import numpy as np
f = bs(P, P, r, sigma, T)
fsim = bs(P, Psim[:,1], r, sigma, T-(1/365))
q = np.sort(fsim['Call'] + Psim[:,0] - f['Call'] - P)
VaR5 = -q[int(p*S) - 1]
print(VaR5)