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


Copyright 2011 - 2020 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 Python
Last updated July 2020

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()
		
Listing 7.1/7.2: Transformation in Julia
Last updated July 2020

x = -3:0.1:3
using Distributions, Plots;
return plot(x, cdf.(Normal(0,1), x), title = "CDF of Standard Normal", legend = false)
		

Listing 7.3/7.4: Various RNs in Python
Last updated July 2020

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))
		
Listing 7.3/7.4: Various RNs in Julia
Last updated July 2020

using Random;
Random.seed!(12);    # set seed
S = 10;
println(S, " draws from Uniform(0,1): \n", rand(Uniform(0,1), S), "\n") # alternatively, rand(S)
println(S, " draws from Normal(0,1): \n", rand(Normal(0,1), S), "\n")   # alternatively, randn(S)
println(S, " draws from Student-t(4): \n", rand(TDist(4), S), "\n")
		

Listing 7.5/7.6: Price bond in Python
Last updated July 2020

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)
		
Listing 7.5/7.6: Price bond in Julia
Last updated July 2020

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[T] += Par         # add par to cash flows
P = sum(cc./((1 .+ yield_c/100).^(1:T))) # calc price
println("Price of the bond: ", round(P, digits = 3), " USD")
		

Listing 7.7/7.8: Simulate yields in Python
Last updated July 2020

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

Random.seed!(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(NaN, (T,S))
for i in 1:S
    ysim[:,i] = yield_c .+ r[i]
end
plot(ysim, title = "Simulated yield curves")
		

Listing 7.9/7.10: Simulate bond prices in Python
Last updated November 2020

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

using Statistics, StatsPlots;
SP = fill(NaN, S)
for i in 1:S         # S simulations
    SP[i] = sum(cc./(1 .+ ysim[:,i]./100).^(1:T))
end
SP .-= (mean(SP) - P) # correct for mean
bar(SP)
S = 50000
r = randn(S) * sigma
ysim = fill(NaN, (T,S))
for i in 1:S
    ysim[:,i] = yield_c .+ r[i]
end
SP = fill(NaN, S)
for i in 1:S
    SP[i] = sum(cc./(1 .+ ysim[:,i]./100).^(1:T))
end
SP .-= (mean(SP) - P)
histogram(SP,nbins=100,normed=true,xlims=(7,13), legend = false, title = "Histogram of simulated bond prices with fitted normal")
res = fit_mle(Normal, SP)
plot!(Normal(res.μ, res.σ), linewidth = 4, legend = false)
		

Listing 7.11/7.12: Black-Scholes valuation in Python
Last updated July 2020

## This calculation uses the Black-Scholes pricing function (Listing 6.1/6.2)
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)
		
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 Python
Last updated July 2020

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)
		
Listing 7.13/7.14: Black-Scholes simulation in Julia
Last updated July 2020

Random.seed!(12)     # set seed
S = 10^6   # number of simulations
ysim = randn(S) * sigma * sqrt(T) .- 0.5 * sigma^2 * T # sim returns, lognorm corrected
F = P0 * exp(r * T) * exp.(ysim)         # simulated future prices
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
println("Simulated call price: ", round(call_sim, digits = 3))
		

Listing 7.15/7.16: Option density plots in Python
Last updated July 2020

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()
		
Listing 7.15/7.16: Option density plots in Julia
Last updated July 2020

histogram(F, bins = 100, normed = true, xlims = (20,80), title = "Simulated prices density", label = false)
res = fit_mle(Normal, F)
plot!(Normal(res.μ, res.σ), linewidth = 4, label = "Fitted normal")
vline!([X], linewidth = 4, color = "black", label = "Strike price")
histogram(fsim, bins = 110, normed = true, xlims = (0,35), title = "Option density", label = false)
vline!([f["Call"]], linewidth = 4, color = "black", label = "Call price")
		

Listing 7.17/7.18: Simulate VaR in Python
Last updated July 2020

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

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 = 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[ceil(Int, p*S)]
println("Simulated VaR", ceil(Int, p*100), "%: ", round(VaR1, digits = 3), " USD")
		

Listing 7.19/7.20: Simulate option VaR in Python
Last updated July 2020

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

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[ceil(Int, p*S)]
println("Simulated Option VaR", ceil(Int, p*100), "%: ", round(VaR2, digits = 3), " USD")
		

Listing 7.21/7.22: Example 7.3 in Python
Last updated July 2020

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

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)
VaR3 = -q[ceil(Int, p*S)]
println("Portfolio Option VaR", ceil(Int, p*100), "%: ", round(VaR3, digits = 3), " USD")
		

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

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.23/7.24: Simulated two-asset returns in Julia
Last updated July 2020

Random.seed!(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 Python
Last updated July 2020

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

K = 2
P = [100 50]         # prices
x = [1 1]  # number of assets
Port = reshape(P * x', 1)[1]   # portfolio at t
Psim = repeat(P, outer = [S,1]).*exp.(y)' # simulated prices
PortSim = reshape(Psim * x', S)          # simulated portfolio value
q = sort(PortSim .- Port)      # simulated P/L
VaR4 = -q[ceil(Int, p * S)]
println("Two-asset VaR", ceil(Int, p*100), "%: ", round(VaR4, digits = 3), " USD")
		

Listing 7.27/7.28: A two-asset case in Python with an option
Last updated July 2020

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

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[ceil(Int, p * S)]
println("Two-asset with option VaR", ceil(Int, p*100), "%: ", round(VaR5, digits = 3), " USD")