Chapter 7. Simulation Methods for VaR for Options and Bonds

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

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.title("CDF of Normal distribution")
plt.show()
plt.close()

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

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

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)

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[T-1] += 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)

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

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.title("Simulated yield curves")
plt.show()
plt.close()

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

S = 8                   # simulation size
SP = np.zeros([S])      # initializing vector with zeros
for i in range(S):      # S simulations
    SP[i] = np.sum(cc/(np.power((1+ysim[:,i]/100), list(range(1,T+1)))))
SP -= (np.mean(SP) - P) # correct for mean
plt.bar(range(1,S+1), SP)
plt.title("Simulated bond prices, S = 8")
plt.show()
plt.close()
S = 50000
np.random.seed(12)
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/(np.power((1+ysim[:,i]/100), list(range(1,T+1)))))
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.title("Histogram for simulated bond prices, S = 50000")
plt.show()
plt.close()

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)

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
print(f)

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)

np.random.seed(12)          # set seed
S = 10**6                   # number of simulations
ysim=np.random.normal(-0.5*sigma**2*T,
                      sigma*np.sqrt(T),size=S) # sim returns, lognorm corrected
F = P0 * np.exp(r * T) * np.exp(ysim)          # sim future prices
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)

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)

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.title("Density of future prices")
plt.show()
plt.close()
plt.hist(fsim, bins = 60, range = (0, 35), density = True)
plt.axvline(x=f['Call'], color='k')
plt.title("Density of Call option prices")
plt.show()
plt.close()

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)

from math import ceil
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[ceil(p*S) - 1]
print(VaR1)

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)

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[ceil(p*S) - 1] 
print(VaR2)

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)

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[ceil(p*S) - 1]
print(VaR3)

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

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

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)

import numpy.matlib
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[ceil(p*S) - 1]
print(VaR4)

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)

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