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

Chapter 7. Simulation Methods for VaR for Options and Bonds

R and Julia

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.

Listing 7.1/7.2
Transformation in R
x = seq(-3, 3, by = 0.1)
plot(x, pnorm(x), type = "l")
Listing 7.1/7.2
Transformation in Julia
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 R
set.seed(12) # set seed
S = 10
runif(S)
rnorm(S)
rt(S,4)
Listing 7.3/7.4
Various RNs in Julia
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 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)
Listing 7.5/7.6
Price bond in Julia
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 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')
Listing 7.7/7.8
Simulate yields in Julia
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 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)))
Listing 7.9/7.10
Simulate bond prices in Julia
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 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)
Listing 7.11/7.12
Simulate bond prices in Julia
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 = X, P = P0, r = r, sigma = sigma, T = T) # analytical call price

Listing 7.13/7.14
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)
Listing 7.13/7.14
Black-Scholes simulation in Julia
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 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)
Listing 7.15/7.16
Option density plots in Julia
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 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)
Listing 7.17/7.18
Simulate VaR in Julia
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 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)
Listing 7.19/7.20
Simulate option VaR in Julia
T = 0.25                            # time to expiration
X = 100                             # strike price
sigma = sqrt(s2 * 250)              # annual volatility
f = bs(X = X,P = P,r = r,sigma = sigma,T = T)               # analytical call price
fsim = bs(X = X,P = Psim, r = r, sigma = sigma, T = 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 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)
Listing 7.21/7.22
Example 7.3 in Julia
X1 = 100
X2 = 110
f1 = bs(X = X1,P = P,r = r,sigma = sigma,T = T)
f2 = bs(X = X2,P = P,r = r,sigma = sigma,T = T)
f2sim = bs(X = X2,P = Psim,r = r,sigma = sigma,T = T-(1/365))
f1sim = bs(X = X1,P = Psim,r = r,sigma = sigma,T = 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 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
Listing 7.23/7.24
Simulated two-asset returns in Julia
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 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)
Listing 7.25/7.26
Two-asset VaR in Julia
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 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)
Listing 7.27/7.28
A two-asset case in Julia with an option
f = bs(X = P[2], P = P[2], r = r, sigma = sigma, T = T)
fsim = bs(X = P[2], P = Psim[:,2], r = r, sigma = sigma, T = 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")


Financial Risk Forecasting
Market risk forecasting with R, Julia, Python and Matlab. Code, lecture slides, implementation notes, seminar assignments and questions.
© All rights reserved, Jon Danielsson,