library(mvtnorm)
library(ggplot2)
source("common/risk.r") # EmpiricalES only; no data or network dependencies25 Simulation methods for risk
Analytical VaR works for stocks. For bonds it relies on a linear duration approximation, and once we add options or any other non-linear payoff, closed-form methods stop being enough.
Every simulated risk figure in this chapter is built the same way. Draw the future price, revalue the position at that price, subtract what it is worth today, and read a quantile off the resulting profit and loss. What changes from section to section is only the portfolio — one stock, then one option, then both, then two of each.
25.1 Data and libraries
The multivariate normal generator is the only addition to the libraries used so far.
import numpy as np
from scipy import stats
import sys
sys.path.insert(0, 'common')
from risk import EmpiricalES # risk helpers only; avoids the data/network importsusing Random, Statistics, Distributions, Printf
include("common/risk.jl"); # EmpiricalES only; functions.jl would pull in HTTP/CSV25.2 Black-Scholes equation
To compute profit and loss for a portfolio containing options, we need a way to reprice each instrument under simulated market scenarios. For stocks, the repricing is trivial — the simulated price is the new value. For options, we must feed each simulated underlying price back through a pricing model to obtain the option’s new value. The Black-Scholes formula does this. It maps a given underlying price, volatility and maturity to an option price, so we can revalue the option under every simulated scenario.
bs = function(K, P, r, sigma, Mat) {
d1 = (log(P/K) + (r + 0.5*sigma^2)*(Mat)) / (sigma*sqrt(Mat))
d2 = d1 - sigma*sqrt(Mat)
Call = P*pnorm(d1) - K*exp(-r*(Mat))*pnorm(d2)
Put = K*exp(-r*(Mat))*pnorm(-d2) - P*pnorm(-d1)
return(list(Call = Call, Put = Put))
}def bs(K, P, r, sigma, Mat):
d1 = (np.log(P/K) + (r + 0.5*sigma**2)*Mat) / (sigma*np.sqrt(Mat))
d2 = d1 - sigma*np.sqrt(Mat)
Call = P*stats.norm.cdf(d1) - K*np.exp(-r*Mat)*stats.norm.cdf(d2)
Put = K*np.exp(-r*Mat)*stats.norm.cdf(-d2) - P*stats.norm.cdf(-d1)
return {'Call': Call, 'Put': Put}function bs(K, P, r, sigma, Mat)
d1 = (log(P/K) + (r + 0.5*sigma^2)*Mat) / (sigma*sqrt(Mat))
d2 = d1 - sigma*sqrt(Mat)
Call = P*cdf(Normal(), d1) - K*exp(-r*Mat)*cdf(Normal(), d2)
Put = K*exp(-r*Mat)*cdf(Normal(), -d2) - P*cdf(Normal(), -d1)
return Dict("Call" => Call, "Put" => Put)
end;25.3 Simulating returns
VaR uses the real-world distribution of returns, not the risk-neutral distribution used in option pricing.
In most of the examples below, we assume normality and use simple returns, but we could just as easily use any of the more complex models from Chapter 23, including those with non-normal conditional distributions and even historical simulation. Under normality the simulated simple return can in theory fall below -1, so the simulated price can go negative, although this is negligible over a one-day horizon. Log returns avoid the problem entirely, since they cannot take the price below zero.
Simple returns are defined by \[ \SimpleReturns_{t+1} = \frac{\Price_{t+1}}{\Price_t}-1 \] We know \(\Price_t\) and simulate \(\SimpleReturns_{t+1}\) to get the simulated future price \(\Price_{t+1}\), so \[\Price_{t+1,s} = \Price_t \times (1+\SimpleReturns_{t+1,s})\] where \(s\) indicates a particular simulation.
Note that \(\Price_{t+1,s}\) is the simulated one-day-ahead price, not the price of a futures contract.
Two volatilities appear in this chapter, and they are not the same quantity. sigma_return is the daily volatility of the real-world return we simulate. sigma_implied is the volatility we hand to Black-Scholes when repricing the option, and in a pricing model that input is an implied volatility, backed out of market option prices under the risk-neutral measure. Annualising an estimate of physical return volatility does not turn it into an implied volatility.
We set the two equal here, sigma_implied = sqrt(250) * sigma_return, because one number keeps the example readable. That is a simplification we are choosing, not an identity, and the separate names are there so it stays visible. In a real system the two are estimated from different data and generally differ, with the gap between them the variance risk premium.
Holding sigma_implied fixed across scenarios costs us something specific. The option is repriced only for the simulated move in the underlying, so the simulation contains no vega. It cannot produce a loss from implied volatility rising while the price is unchanged. That matters because the two usually move together in exactly the direction that hurts — equity falls and implied volatility jumps — so a portfolio with material option exposure will have its risk understated here. The extension is to simulate the underlying price and a shock to the implied-volatility surface jointly, then reprice at the shocked surface.
Otherwise, the interest rate r is flat and unchanging. Where a bivariate portfolio is used, returns are drawn from a fixed covariance matrix rather than a time-varying one. Options are simulated one day ahead, so maturity decreases by 1/365 of a year (calendar days), while volatility is annualised using sqrt(250) (trading days).
25.4 Holdings
We start by holding one unit of each asset. For a single stock, the current portfolio value is P. This assumption can easily be relaxed to handle arbitrary position sizes and portfolio weights.
25.5 One asset
25.5.1 Setup
par = list()
par$probability = 0.05
par$S = 1000
par$P = 100
par$r = 0.05
par$K = 99
par$Mat = 1
par$sigma_return = 0.01 # daily volatility of the simulated return
par$sigma_implied = sqrt(250) * par$sigma_return # annualised volatility used to price the optionpar_py = {
'probability': 0.05,
'S': 1000,
'P': 100,
'r': 0.05,
'K': 99,
'Mat': 1,
'sigma_return': 0.01, # daily volatility of the simulated return
'sigma_implied': np.sqrt(250) * 0.01
}par_jl = Dict(
"probability" => 0.05,
"S" => 1000,
"P" => 100,
"r" => 0.05,
"K" => 99,
"Mat" => 1,
"sigma_return" => 0.01, # daily volatility of the simulated return
"sigma_implied" => sqrt(250) * 0.01
);25.5.2 Simulate prices
We create the function Sim.Prices(par) to implement the simulation. It takes the par list as an argument and returns a vector of simulated prices. The function could be modified to return the Returns vector, to draw from distributions other than the normal or to incorporate a mean model.
The number of simulations S is read from par by default. An optional S argument overrides this value when provided, via if(!is.null(S)) par$S=S. The seed argument controls reproducibility. Passing NULL disables seeding.
Sim.Prices = function(par, seed = 888, S = NULL) {
if (!is.null(seed)) set.seed(seed)
if (!is.null(S)) par$S = S
Returns = rnorm(n = par$S, mean = 0, sd = par$sigma_return)
Prices = par$P * (1 + Returns)
return(Prices)
}def Sim_Prices(par, seed=888, S=None):
if seed is not None:
np.random.seed(seed)
n = S if S is not None else par['S']
Returns = np.random.normal(0, par['sigma_return'], n)
Prices = par['P'] * (1 + Returns)
return Pricesfunction Sim_Prices(par; seed=888, S=nothing)
if seed !== nothing
Random.seed!(seed)
end
n = S !== nothing ? S : par["S"]
Returns = randn(n) * par["sigma_return"]
Prices = par["P"] .* (1 .+ Returns)
return Prices
end;Here is one example.
25.5.3 Design
The general setup below is that we first calculate the current portfolio value, called today, and then do the simulations to get the simulated vector of tomorrow’s portfolio value, called tomorrow. The difference between the two is profit and loss, PL. The VaR is then the negative of a quantile of the PL vector, while ES is the negative of the average from that quantile to the minimum.
25.5.4 VaR and ES for one stock
cat("Number of simulations:", par$S, "\n")
today = par$P
tomorrow = Sim.Prices(par, seed = 8888)
PL = tomorrow - today
k = ceiling(par$probability * par$S)
VaR = -sort(PL)[k]
cat("VaR = $", round(VaR, 4), "\n", sep = "")
ES = -EmpiricalES(sort(PL), par$probability)
cat("ES = $", round(ES, 4), "\n", sep = "")
tomorrow = Sim.Prices(par, seed = 666)
PL = tomorrow - today
k = ceiling(par$probability * par$S)
VaR = -sort(PL)[k]
cat("VaR = $", round(VaR, 4), "\n", sep = "")
ES = -EmpiricalES(sort(PL), par$probability)
cat("ES = $", round(ES, 4), "\n", sep = "")Number of simulations: 1000
VaR = $1.8085
ES = $2.0912
VaR = $1.6823
ES = $2.0817print(f"Number of simulations: {par_py['S']}")
today_py = par_py['P']
tomorrow_py = Sim_Prices(par_py, seed=8888)
PL_py = tomorrow_py - today_py
k = int(np.ceil(par_py['probability'] * par_py['S']))
VaR_py = -np.sort(PL_py)[k - 1]
print(f"VaR = ${VaR_py:.4f}")
ES_py = -EmpiricalES(np.sort(PL_py), par_py['probability'])
print(f"ES = ${ES_py:.4f}")
tomorrow_py = Sim_Prices(par_py, seed=666)
PL_py = tomorrow_py - today_py
k = int(np.ceil(par_py['probability'] * par_py['S']))
VaR_py = -np.sort(PL_py)[k - 1]
print(f"VaR = ${VaR_py:.4f}")
ES_py = -EmpiricalES(np.sort(PL_py), par_py['probability'])
print(f"ES = ${ES_py:.4f}")Number of simulations: 1000
VaR = $1.7325
ES = $2.2062
VaR = $1.7044
ES = $2.0813using Statistics
println("Number of simulations: ", par_jl["S"]);
today_jl = par_jl["P"];
tomorrow_jl = Sim_Prices(par_jl, seed=8888);
PL_jl = tomorrow_jl .- today_jl;
k = ceil(Int, par_jl["probability"] * par_jl["S"]);
VaR_jl = -sort(PL_jl)[k];
println("VaR = \$", round(VaR_jl, digits=4));
ES_jl = -EmpiricalES(sort(PL_jl), par_jl["probability"]);
println("ES = \$", round(ES_jl, digits=4));
tomorrow_jl = Sim_Prices(par_jl, seed=666);
PL_jl = tomorrow_jl .- today_jl;
k = ceil(Int, par_jl["probability"] * par_jl["S"]);
VaR_jl = -sort(PL_jl)[k];
println("VaR = \$", round(VaR_jl, digits=4));
ES_jl = -EmpiricalES(sort(PL_jl), par_jl["probability"]);
println("ES = \$", round(ES_jl, digits=4))Number of simulations: 1000
VaR = $1.6665
ES = $1.9869
VaR = $1.5942
ES = $1.9746Two seeds already give noticeably different VaR estimates. With only 1,000 draws, Monte Carlo error is still large. We analyse the number of simulations below.
25.5.5 VaR and ES for one option
When we have one option, we need to run the Black-Scholes equation for both today’s price and tomorrow’s simulated price.
x = Sim.Prices(par)
type = "Call"
today = bs(
K = par$K,
P = par$P,
r = par$r,
sigma = par$sigma_implied,
Mat = par$Mat
)[[type]]
tomorrow = bs(
K = par$K,
P = x,
r = par$r,
sigma = par$sigma_implied,
Mat = par$Mat - 1/365
)[[type]]
PL = tomorrow - today
k = ceiling(par$probability * par$S)
VaR = -sort(PL)[k]
cat("VaR = $", round(VaR, 4), "\n", sep = "")
ES = -EmpiricalES(sort(PL), par$probability)
cat("ES = $", round(ES, 4), "\n", sep = "")VaR = $1.0949
ES = $1.3439x_py = Sim_Prices(par_py)
type_py = 'Call'
today_py = bs(
K=par_py['K'],
P=par_py['P'],
r=par_py['r'],
sigma=par_py['sigma_implied'],
Mat=par_py['Mat']
)[type_py]
tomorrow_py = bs(
K=par_py['K'],
P=x_py,
r=par_py['r'],
sigma=par_py['sigma_implied'],
Mat=par_py['Mat'] - 1/365
)[type_py]
PL_py = tomorrow_py - today_py
k = int(np.ceil(par_py['probability'] * par_py['S']))
VaR_py = -np.sort(PL_py)[k - 1]
print(f"VaR = ${VaR_py:.4f}")
ES_py = -EmpiricalES(np.sort(PL_py), par_py['probability'])
print(f"ES = ${ES_py:.4f}")VaR = $1.0586
ES = $1.3123using Statistics
x_jl = Sim_Prices(par_jl);
type_jl = "Call";
today_jl = bs(
par_jl["K"],
par_jl["P"],
par_jl["r"],
par_jl["sigma_implied"],
par_jl["Mat"]
)[type_jl];
tomorrow_jl = bs.(
par_jl["K"],
x_jl,
par_jl["r"],
par_jl["sigma_implied"],
par_jl["Mat"] - 1/365
);
tomorrow_jl = [t[type_jl] for t in tomorrow_jl];
PL_jl = tomorrow_jl .- today_jl;
k = ceil(Int, par_jl["probability"] * par_jl["S"]);
VaR_jl = -sort(PL_jl)[k];
println("VaR = \$", round(VaR_jl, digits=4))
ES_jl = -EmpiricalES(sort(PL_jl), par_jl["probability"]);
println("ES = \$", round(ES_jl, digits=4))VaR = $1.0903
ES = $1.3063The option carries markedly less risk than the stock it is written on. A call’s delta is below one, so it loses less than the underlying when the price falls, and the most it can lose is the premium paid for it.
25.5.6 VaR and ES for one stock and option
If we hold the stock and an option written on it, one simulated price vector is enough for both positions.
x = Sim.Prices(par)
type = "Call"
today = par$P + bs(
K = par$K,
P = par$P,
r = par$r,
sigma = par$sigma_implied,
Mat = par$Mat
)[[type]]
tomorrow = x + bs(
K = par$K,
P = x,
r = par$r,
sigma = par$sigma_implied,
Mat = par$Mat - 1/365
)[[type]]
PL = tomorrow - today
k = ceiling(par$probability * par$S)
VaR = -sort(PL)[k]
cat("VaR = $", round(VaR, 4), "\n", sep = "")
ES = -EmpiricalES(sort(PL), par$probability)
cat("ES = $", round(ES, 4), "\n", sep = "")VaR = $2.7359
ES = $3.3803x_py = Sim_Prices(par_py)
type_py = 'Call'
today_py = par_py['P'] + bs(
K=par_py['K'],
P=par_py['P'],
r=par_py['r'],
sigma=par_py['sigma_implied'],
Mat=par_py['Mat']
)[type_py]
tomorrow_py = x_py + bs(
K=par_py['K'],
P=x_py,
r=par_py['r'],
sigma=par_py['sigma_implied'],
Mat=par_py['Mat'] - 1/365
)[type_py]
PL_py = tomorrow_py - today_py
k = int(np.ceil(par_py['probability'] * par_py['S']))
VaR_py = -np.sort(PL_py)[k - 1]
print(f"VaR = ${VaR_py:.4f}")
ES_py = -EmpiricalES(np.sort(PL_py), par_py['probability'])
print(f"ES = ${ES_py:.4f}")VaR = $2.6429
ES = $3.2987using Statistics
x_jl = Sim_Prices(par_jl);
type_jl = "Call";
today_jl = par_jl["P"] + bs(
par_jl["K"],
par_jl["P"],
par_jl["r"],
par_jl["sigma_implied"],
par_jl["Mat"]
)[type_jl];
tomorrow_bs = bs.(
par_jl["K"],
x_jl,
par_jl["r"],
par_jl["sigma_implied"],
par_jl["Mat"] - 1/365
);
tomorrow_jl = x_jl .+ [t[type_jl] for t in tomorrow_bs];
PL_jl = tomorrow_jl .- today_jl;
k = ceil(Int, par_jl["probability"] * par_jl["S"]);
VaR_jl = -sort(PL_jl)[k];
println("VaR = \$", round(VaR_jl, digits=4))
ES_jl = -EmpiricalES(sort(PL_jl), par_jl["probability"]);
println("ES = \$", round(ES_jl, digits=4))VaR = $2.7242
ES = $3.2827Holding both gives a figure close to the sum of the two separately. Both positions are written on the same underlying, so they fall together and there is nothing to diversify. Diversification needs a second risk factor, and this portfolio has only one.
25.6 Bivariate simulation
If we hold a portfolio of two stocks and options on each, we need to simulate from the bivariate normal distribution. In this case, we draw from the bivariate normal with each language’s multivariate normal generator, but it would be straightforward to do it manually with a Cholesky decomposition.
25.6.1 Setup
par = list()
par$probability = 0.05
par$S = 1000
par$r = 0.05
par$P = c(100, 25)
par$Sigma = rbind(
c(0.01, 0.005),
c(0.005, 0.02)
)
par$Mat = c(0.5, 1)
par$K = c(90, 30)par_py = {
'probability': 0.05,
'S': 1000,
'r': 0.05,
'P': np.array([100, 25]),
'Sigma': np.array([[0.01, 0.005], [0.005, 0.02]]),
'Mat': np.array([0.5, 1]),
'K': np.array([90, 30])
}par_jl = Dict(
"probability" => 0.05,
"S" => 1000,
"r" => 0.05,
"P" => [100.0, 25.0],
"Sigma" => [0.01 0.005; 0.005 0.02],
"Mat" => [0.5, 1.0],
"K" => [90.0, 30.0]
);25.6.2 Simulate bivariate normal
Use the number of simulations and the covariance matrix, par$Sigma, to generate correlated returns and then prices. Note that par$Sigma is the covariance matrix, unlike the scalar standard deviation par$sigma of the univariate examples. The marginal standard deviations are the square roots of its diagonal.
print(par$Sigma)
set.seed(42)
r = rmvnorm(10, sigma = par$Sigma)
print(r)
print(cov(r))
print(cor(r)) [,1] [,2]
[1,] 0.010 0.005
[2,] 0.005 0.020
[,1] [,2]
[1,] 0.1221431166 -0.05012218
[2,] 0.0488171252 0.09614529
[3,] 0.0372884402 -0.00633415
[4,] 0.1457757997 0.01856928
[5,] 0.1960031758 0.03370345
[6,] 0.1756832449 0.34723901
[7,] -0.1416428126 -0.06821367
[8,] 0.0003487389 0.08613028
[9,] -0.0836887965 -0.37747893
[10,] -0.2108029484 0.13325900
[,1] [,2]
[1,] 0.019147085 0.008407915
[2,] 0.008407915 0.033515472
[,1] [,2]
[1,] 1.0000000 0.3319054
[2,] 0.3319054 1.0000000print(par_py['Sigma'])
np.random.seed(42)
r_py = np.random.multivariate_normal([0, 0], par_py['Sigma'], 10)
print(r_py)
print(np.cov(r_py, rowvar=False))
print(np.corrcoef(r_py, rowvar=False))[[0.01 0.005]
[0.005 0.02 ]]
[[ 0.01686504 0.0728878 ]
[ 0.1621171 0.03699968]
[-0.03257387 -0.02416021]
[ 0.15291668 0.19060327]
[ 0.01794358 -0.08292574]
[-0.06466049 -0.04773608]
[-0.14364254 0.09840714]
[-0.14432368 -0.21759253]
[-0.03173021 -0.14972402]
[-0.16780884 -0.07650503]]
[[0.01340282 0.00765833]
[0.00765833 0.0148957 ]]
[[1. 0.54200756]
[0.54200756 1. ]]println(par_jl["Sigma"])
Random.seed!(42)
mvdist = MvNormal(par_jl["Sigma"])
r_jl = rand(mvdist, 10)'
println(r_jl)
println(cov(r_jl))
println(cor(r_jl))[0.01 0.005; 0.005 0.02]
[-0.03633574814517775 0.015133829334874442; -0.031498797116895606 -0.05692422100225221; 0.08163067649323275 0.10388189792036637; -0.08595553820616213 -0.23734632890561533; -0.2114334831130985 -0.09992497207064548; -0.08253345499750689 0.06989299560375428; 0.04338858743048611 -0.030617480140826672; 0.051713088145035904 0.2173062695003553; -0.020661311448119266 -0.05143820620488855; -0.00404734002483549 0.01183620297442292]
[0.00712605952489706 0.0060762374903670195; 0.0060762374903670195 0.015118984958544492]
[1.0 0.5853945833699821; 0.5853945833699821 1.0]25.6.3 Simulate bivariate prices
Sim.BivarPrices = function(par, seed = 666, S = NULL) {
if (!is.null(seed)) set.seed(seed)
if (!is.null(S)) par$S = S
Returns = rmvnorm(par$S, sigma = par$Sigma)
Prices1 = par$P[1] * (1 + Returns[, 1])
Prices2 = par$P[2] * (1 + Returns[, 2])
return(cbind(Prices1, Prices2))
}def Sim_BivarPrices(par, seed=666, S=None):
if seed is not None:
np.random.seed(seed)
n = S if S is not None else par['S']
Returns = np.random.multivariate_normal([0, 0], par['Sigma'], n)
Prices1 = par['P'][0] * (1 + Returns[:, 0])
Prices2 = par['P'][1] * (1 + Returns[:, 1])
return np.column_stack([Prices1, Prices2])function Sim_BivarPrices(par; seed=666, S=nothing)
if seed !== nothing
Random.seed!(seed)
end
n = S !== nothing ? S : par["S"]
mvdist = MvNormal(par["Sigma"])
Returns = rand(mvdist, n)'
Prices1 = par["P"][1] .* (1 .+ Returns[:, 1])
Prices2 = par["P"][2] .* (1 .+ Returns[:, 2])
return hcat(Prices1, Prices2)
end;25.6.4 Two stocks
today = sum(par$P)
tomorrow = rowSums(Sim.BivarPrices(par))
PL = tomorrow - today
k = ceiling(par$probability * par$S)
VaR = -sort(PL)[k]
cat("VaR = $", round(VaR, 4), "\n", sep = "")
ES = -EmpiricalES(sort(PL), par$probability)
cat("ES = $", round(ES, 4), "\n", sep = "")VaR = $18.5912
ES = $23.0376today_py = np.sum(par_py['P'])
tomorrow_py = np.sum(Sim_BivarPrices(par_py), axis=1)
PL_py = tomorrow_py - today_py
k = int(np.ceil(par_py['probability'] * par_py['S']))
VaR_py = -np.sort(PL_py)[k - 1]
print(f"VaR = ${VaR_py:.4f}")
ES_py = -EmpiricalES(np.sort(PL_py), par_py['probability'])
print(f"ES = ${ES_py:.4f}")VaR = $19.0316
ES = $24.6802using Statistics
today_jl = sum(par_jl["P"]);
tomorrow_jl = sum(Sim_BivarPrices(par_jl), dims=2);
PL_jl = tomorrow_jl .- today_jl;
k = ceil(Int, par_jl["probability"] * par_jl["S"]);
VaR_jl = -sort(vec(PL_jl))[k];
println("VaR = \$", round(VaR_jl, digits=4))
ES_jl = -mean(sort(vec(PL_jl))[1:k]);
println("ES = \$", round(ES_jl, digits=4))VaR = $17.3065
ES = $22.829This is the first portfolio with two risk factors, and the covariance matrix now does real work. The two return series are drawn together rather than independently, so how closely JPM and INTC move decides the portfolio loss. Setting the off-diagonal to zero would give a different answer.
25.6.5 Two stocks and two options
25.6.5.1 Today
type = c("Call", "Put")
today = sum(par$P)
for (i in 1:2) {
today = today + bs(
K = par$K[i],
P = par$P[i],
r = par$r,
sigma = sqrt(250) * sqrt(par$Sigma[i, i]),
Mat = par$Mat[i]
)[[type[i]]]
}
cat("Today's portfolio value: $", round(today, 4), "\n", sep = "")Today's portfolio value: $192.653type_py = ['Call', 'Put']
today_py = np.sum(par_py['P'])
for i in range(2):
today_py = today_py + bs(
K=par_py['K'][i],
P=par_py['P'][i],
r=par_py['r'],
sigma=np.sqrt(250) * np.sqrt(par_py['Sigma'][i, i]),
Mat=par_py['Mat'][i]
)[type_py[i]]
print(f"Today's portfolio value: ${today_py:.4f}")Today's portfolio value: $192.6530type_jl = ["Call", "Put"];
today_jl = sum(par_jl["P"]);
for i in 1:2
global today_jl = today_jl + bs(
par_jl["K"][i],
par_jl["P"][i],
par_jl["r"],
sqrt(250) * sqrt(par_jl["Sigma"][i, i]),
par_jl["Mat"][i]
)[type_jl[i]]
end
println("Today's portfolio value: \$", round(today_jl, digits=4))Today's portfolio value: $192.65325.6.5.2 Tomorrow
x = Sim.BivarPrices(par)
tomorrow = 0
for (i in 1:2) {
tomorrow = tomorrow + x[, i] + bs(
K = par$K[i],
P = x[, i],
r = par$r,
sigma = sqrt(250) * sqrt(par$Sigma[i, i]),
Mat = par$Mat[i] - 1/365
)[[type[i]]]
}
PL = tomorrow - today
k = ceiling(par$probability * par$S)
VaR = -sort(PL)[k]
cat("VaR = $", round(VaR, 4), "\n", sep = "")
ES = -EmpiricalES(sort(PL), par$probability)
cat("ES = $", round(ES, 4), "\n", sep = "")VaR = $30.0299
ES = $36.3202x_py = Sim_BivarPrices(par_py)
tomorrow_py = np.zeros(par_py['S'])
for i in range(2):
tomorrow_py = tomorrow_py + x_py[:, i] + bs(
K=par_py['K'][i],
P=x_py[:, i],
r=par_py['r'],
sigma=np.sqrt(250) * np.sqrt(par_py['Sigma'][i, i]),
Mat=par_py['Mat'][i] - 1/365
)[type_py[i]]
PL_py = tomorrow_py - today_py
k = int(np.ceil(par_py['probability'] * par_py['S']))
VaR_py = -np.sort(PL_py)[k - 1]
print(f"VaR = ${VaR_py:.4f}")
ES_py = -EmpiricalES(np.sort(PL_py), par_py['probability'])
print(f"ES = ${ES_py:.4f}")VaR = $30.0951
ES = $38.5074using Statistics
x_jl = Sim_BivarPrices(par_jl);
tomorrow_jl = zeros(par_jl["S"]);
for i in 1:2
bs_vals = bs.(
par_jl["K"][i],
x_jl[:, i],
par_jl["r"],
sqrt(250) * sqrt(par_jl["Sigma"][i, i]),
par_jl["Mat"][i] - 1/365
)
global tomorrow_jl = tomorrow_jl .+ x_jl[:, i] .+ [b[type_jl[i]] for b in bs_vals]
end
PL_jl = tomorrow_jl .- today_jl;
k = ceil(Int, par_jl["probability"] * par_jl["S"]);
VaR_jl = -sort(PL_jl)[k];
println("VaR = \$", round(VaR_jl, digits=4))
ES_jl = -EmpiricalES(sort(PL_jl), par_jl["probability"]);
println("ES = \$", round(ES_jl, digits=4))VaR = $28.3792
ES = $35.7204Adding a call on each stock raises VaR by more than half. The options bring no new risk factor — each is written on a stock we already hold — so what they add is leverage on the same two sources of risk rather than exposure to a third.
25.7 Sim VaR vs. true VaR
When the portfolio contains only the stock, the true VaR is known. That gives us a clean benchmark for the simulated VaR. As the number of simulations increases, these two values should converge. We can also use this to ascertain how many simulations we need.
25.7.1 Setup
par = list()
par$probability = 0.01
par$P = 100
par$sigma_return = 0.01 # daily volatility of the simulated return
par$sigma_implied = sqrt(250) * par$sigma_return # annualised volatility used to price the optionpar_py = {
'probability': 0.01,
'P': 100,
'sigma_return': 0.01, # daily volatility of the simulated return
'sigma_implied': np.sqrt(250) * 0.01
}par_jl = Dict(
"probability" => 0.01,
"P" => 100,
"sigma_return" => 0.01, # daily volatility of the simulated return
"sigma_implied" => sqrt(250) * 0.01
);25.7.2 VaR for one stock
To facilitate the analysis, we create a function simVaR() that returns simulated VaR. Because it allows us to vary the seed and number of simulations, we can get an idea of the accuracy.
simVaR = function(par, S, seed = 14) {
set.seed(seed)
par$S = S
tomorrow = Sim.Prices(par, seed = seed)
PL = tomorrow - par$P
VaR = -sort(PL)[ceiling(par$probability * par$S)]
return(VaR)
}def simVaR_py(par, S, seed=14):
par_copy = par.copy()
par_copy['S'] = S
tomorrow = Sim_Prices(par_copy, seed=seed, S=S)
PL = tomorrow - par['P']
VaR = -np.sort(PL)[int(np.ceil(par['probability'] * S)) - 1]
return VaRfunction simVaR_jl(par, S; seed=14)
par_copy = copy(par)
par_copy["S"] = S
tomorrow = Sim_Prices(par_copy, seed=seed, S=S)
PL = tomorrow .- par["P"]
VaR = -sort(PL)[ceil(Int, par["probability"] * S)]
return VaR
end;Here, we show the simulated and true VaR.
VaR = simVaR(par = par, S = 100)
cat("VaR = $", round(VaR, 4), "\n", sep = "")
trueVaR = -par$P * qnorm(par$probability) * par$sigma_return
cat("trueVaR = $", round(trueVaR, 4), "\n", sep = "")VaR = $2.137
trueVaR = $2.3263VaR_py = simVaR_py(par_py, S=100)
print(f"VaR = ${VaR_py:.4f}")
trueVaR_py = -par_py['P'] * stats.norm.ppf(par_py['probability']) * par_py['sigma_return']
print(f"trueVaR = ${trueVaR_py:.4f}")VaR = $2.4385
trueVaR = $2.3263VaR_jl = simVaR_jl(par_jl, 100);
println("VaR = \$", round(VaR_jl, digits=4))
trueVaR_jl = -par_jl["P"] * quantile(Normal(), par_jl["probability"]) * par_jl["sigma_return"];
println("trueVaR = \$", round(trueVaR_jl, digits=4))VaR = $2.232
trueVaR = $2.3263We can repeat this:
trueVaR = -par$P * qnorm(par$probability) * par$sigma_return
N = 10
VaR = vector(length = N)
for (i in 1:N) VaR[i] = simVaR(par = par, S = 100, seed = i)
Error = VaR - trueVaR
cat("True VaR:", round(trueVaR, 4), "\n")
for (i in 1:N) cat("VaR[", i, "]: ", round(VaR[i], 4), ", Error: ", round(Error[i], 4), "\n", sep = "")True VaR: 2.3263
VaR[1]: 2.2147, Error: -0.1116
VaR[2]: 2.4517, Error: 0.1254
VaR[3]: 2.2654, Error: -0.0609
VaR[4]: 1.7974, Error: -0.529
VaR[5]: 2.184, Error: -0.1424
VaR[6]: 1.9523, Error: -0.374
VaR[7]: 1.7859, Error: -0.5405
VaR[8]: 3.0145, Error: 0.6882
VaR[9]: 2.6177, Error: 0.2914
VaR[10]: 2.1853, Error: -0.1411trueVaR_py = -par_py['P'] * stats.norm.ppf(par_py['probability']) * par_py['sigma_return']
N = 10
VaR_arr = np.zeros(N)
for i in range(N):
VaR_arr[i] = simVaR_py(par_py, S=100, seed=i+1)
Error = VaR_arr - trueVaR_py
print(f"True VaR: {trueVaR_py:.4f}")
for i in range(N):
print(f"VaR[{i+1}]: {VaR_arr[i]:.4f}, Error: {Error[i]:.4f}")True VaR: 2.3263
VaR[1]: 2.3015, Error: -0.0248
VaR[2]: 2.6594, Error: 0.3331
VaR[3]: 2.9157, Error: 0.5894
VaR[4]: 2.3798, Error: 0.0534
VaR[5]: 2.8597, Error: 0.5333
VaR[6]: 2.4868, Error: 0.1604
VaR[7]: 2.2883, Error: -0.0380
VaR[8]: 3.1349, Error: 0.8086
VaR[9]: 2.7934, Error: 0.4670
VaR[10]: 2.1317, Error: -0.1946trueVaR_jl = -par_jl["P"] * quantile(Normal(), par_jl["probability"]) * par_jl["sigma_return"];
N = 10;
VaR_arr = zeros(N);
for i in 1:N
VaR_arr[i] = simVaR_jl(par_jl, 100, seed=i)
end
Error = VaR_arr .- trueVaR_jl;
println("True VaR: ", round(trueVaR_jl, digits=4))
for i in 1:N
println("VaR[$i]: ", round(VaR_arr[i], digits=4), ", Error: ", round(Error[i], digits=4))
endTrue VaR: 2.3263
VaR[1]: 2.7903, Error: 0.4639
VaR[2]: 2.4537, Error: 0.1274
VaR[3]: 2.2304, Error: -0.096
VaR[4]: 2.9877, Error: 0.6614
VaR[5]: 2.7083, Error: 0.382
VaR[6]: 2.6721, Error: 0.3458
VaR[7]: 2.5861, Error: 0.2597
VaR[8]: 2.3542, Error: 0.0278
VaR[9]: 2.711, Error: 0.3846
VaR[10]: 2.3845, Error: 0.0582We can also vary the number of simulations.
S_vals = c(1e3, 1e4, 1e5, 1e6, 1e7)
VaR_r = S_vals
for (i in 1:length(S_vals)) VaR_r[i] = simVaR(par = par, S = S_vals[i])
cat("S\t\tVaR\t\tError\n")
for (i in 1:length(S_vals)) cat(S_vals[i], "\t\t", round(VaR_r[i], 4), "\t\t", round(VaR_r[i] - trueVaR, 4), "\n", sep = "")S VaR Error
1000 2.3279 0.0015
10000 2.3921 0.0657
1e+05 2.3157 -0.0106
1e+06 2.326 -4e-04
1e+07 2.3268 4e-04S_vals_py = [1000, 10000, 100000, 1000000, 10000000]
VaR_py_arr = np.zeros(len(S_vals_py))
for i, s in enumerate(S_vals_py):
VaR_py_arr[i] = simVaR_py(par_py, S=s)
print("S\t\tVaR\t\tError")
for i, s in enumerate(S_vals_py):
print(f"{s}\t\t{VaR_py_arr[i]:.4f}\t\t{VaR_py_arr[i] - trueVaR_py:.4f}")S VaR Error
1000 2.5206 0.1943
10000 2.4015 0.0752
100000 2.3329 0.0065
1000000 2.3337 0.0074
10000000 2.3298 0.0035S_vals_jl = [1000, 10000, 100000, 1000000, 10000000];
VaR_jl_arr = zeros(length(S_vals_jl));
for (i, s) in enumerate(S_vals_jl)
VaR_jl_arr[i] = simVaR_jl(par_jl, s)
end
println("S\t\tVaR\t\tError")
for (i, s) in enumerate(S_vals_jl)
println(s, "\t\t", round(VaR_jl_arr[i], digits=4), "\t\t", round(VaR_jl_arr[i] - trueVaR_jl, digits=4))
endS VaR Error
1000 2.4567 0.1304
10000 2.3646 0.0382
100000 2.3375 0.0112
1000000 2.3254 -0.0009
10000000 2.3256 -0.0007The figure below tracks VaR against the number of simulations on a finer grid, in R only, and shows how the estimate settles towards the true VaR as the sample grows. A single cumulative path like this is a visual check, not a measure of Monte Carlo error, since it can flatten out by chance.
25.7.3 Cross-language comparison
The table below collects the VaR estimate from each language at the five simulation counts used above, alongside the true VaR. All three languages compute the same order statistic, so agreement between the R, Python and Julia columns checks the implementations rather than delivering a new result.
S VaR_R VaR_Python VaR_Julia True_VaR
1 1e+03 2.3279 2.5206 2.4567 2.3263
2 1e+04 2.3921 2.4015 2.3646 2.3263
3 1e+05 2.3157 2.3329 2.3375 2.3263
4 1e+06 2.3260 2.3337 2.3254 2.3263
5 1e+07 2.3268 2.3298 2.3256 2.326325.8 How precise is a simulated VaR?
The convergence path above shows the estimate settling down, but a single path is not an error measure. It shows one realisation wandering toward the answer, and cannot say how far a fresh run of the same size would land from it. For that we need the sampling distribution of the estimator, not one draw from it.
For a quantile estimated from \(\NumberSims\) independent draws there is a standard large-sample result,
\[\operatorname{SE}\left(\widehat{\VaR}\right) \approx \frac{\sqrt{\probability(1-\probability)}}{\PDF(\Quantile_\probability)\sqrt{\NumberSims}},\]
where \(\PDF\) is the density of the P&L distribution and \(\Quantile_\probability\) its \(\probability\)-quantile. The \(\sqrt{\NumberSims}\) in the denominator is the familiar part. Cutting the standard error by a factor of ten costs a hundred times the simulations.
The density term does the less obvious work. It sits in the denominator, so the flatter the distribution is where we are trying to read it, the worse the estimate. Deeper in the tail the density is smaller, so precision degrades faster than the number of tail observations alone would suggest — doubling \(\NumberSims\) at \(\probability = 0.1\%\) does not buy what it buys at \(\probability = 5\%\).
The formula rests on three conditions that are easy to lose. The draws must be independent, which rules out antithetic or control variates, common random numbers, and resampling schemes such as filtered historical simulation. The distribution must be continuous with a positive density at the quantile. Where any of those fail, estimate the precision empirically, by running independent batches and looking at the spread of their estimates.
Batches also work for ES, which has no comparably simple closed form. Below we do both — the analytical standard error, and the spread across independent batches, which needs no conditions beyond the batches being independent of each other.
S_prec = 10000
B = 200
# Analytical. P&L is P * Returns, so its standard deviation is sd_PL and its
# p-quantile sits at q_p * sd_PL. The density there is dnorm(q_p) / sd_PL:
# the standard normal density at the standardised point, rescaled. Evaluating
# dnorm at the standardised q_p under scale sd_PL would be wrong, and happens
# to coincide only when sd_PL is 1, which it is for these parameters.
q_p = qnorm(par$probability)
sd_PL = par$P * par$sigma_return
dens_at_q = dnorm(q_p) / sd_PL
se_analytic = sqrt(par$probability * (1 - par$probability)) /
(dens_at_q * sqrt(S_prec))
# Empirical: B independent batches, each a fresh simulation of size S_prec
VaR_batch = sapply(1:B, function(b) {
PL = Sim.Prices(par, seed = 1000 + b, S = S_prec) - par$P
-sort(PL)[ceiling(par$probability * S_prec)]
})
ES_batch = sapply(1:B, function(b) {
PL = Sim.Prices(par, seed = 1000 + b, S = S_prec) - par$P
-EmpiricalES(sort(PL), par$probability)
})
cat(sprintf("VaR standard error, analytical: %.4f\n", se_analytic))
cat(sprintf("VaR standard error, %d batches: %.4f\n", B, sd(VaR_batch)))
cat(sprintf("ES standard error, %d batches: %.4f\n", B, sd(ES_batch)))
cat(sprintf("VaR mean over batches: %.4f (true %.4f)\n", mean(VaR_batch), trueVaR))VaR standard error, analytical: 0.0373
VaR standard error, 200 batches: 0.0409
ES standard error, 200 batches: 0.0475
VaR mean over batches: 2.3277 (true 2.3263)S_prec = 10000
B = 200
q_p = stats.norm.ppf(par_py['probability'])
sd_PL = par_py['P'] * par_py['sigma_return']
dens_at_q = stats.norm.pdf(q_p) / sd_PL # density of P&L at its p-quantile
se_analytic = (np.sqrt(par_py['probability'] * (1 - par_py['probability'])) /
(dens_at_q * np.sqrt(S_prec)))
k_prec = int(np.ceil(par_py['probability'] * S_prec))
VaR_batch = np.empty(B)
ES_batch = np.empty(B)
for b in range(B):
PL = Sim_Prices(par_py, seed=1000 + b, S=S_prec) - par_py['P']
PLs = np.sort(PL)
VaR_batch[b] = -PLs[k_prec - 1]
ES_batch[b] = -EmpiricalES(PLs, par_py['probability'])
print(f"VaR standard error, analytical: {se_analytic:.4f}")
print(f"VaR standard error, {B} batches: {np.std(VaR_batch, ddof=1):.4f}")
print(f"ES standard error, {B} batches: {np.std(ES_batch, ddof=1):.4f}")
print(f"VaR mean over batches: {np.mean(VaR_batch):.4f} (true {trueVaR_py:.4f})")VaR standard error, analytical: 0.0373
VaR standard error, 200 batches: 0.0401
ES standard error, 200 batches: 0.0437
VaR mean over batches: 2.3280 (true 2.3263)S_prec = 10000
B = 200
q_p = quantile(Normal(), par_jl["probability"])
sd_PL = par_jl["P"] * par_jl["sigma_return"]
dens_at_q = pdf(Normal(), q_p) / sd_PL # density of P&L at its p-quantile
se_analytic = sqrt(par_jl["probability"] * (1 - par_jl["probability"])) /
(dens_at_q * sqrt(S_prec))
k_prec = ceil(Int, par_jl["probability"] * S_prec)
VaR_batch = zeros(B)
ES_batch = zeros(B)
for b in 1:B
PL = Sim_Prices(par_jl, seed=1000 + b, S=S_prec) .- par_jl["P"]
PLs = sort(PL)
VaR_batch[b] = -PLs[k_prec]
ES_batch[b] = -EmpiricalES(PLs, par_jl["probability"])
end
@printf("VaR standard error, analytical: %.4f\n", se_analytic)
@printf("VaR standard error, %d batches: %.4f\n", B, std(VaR_batch))
@printf("ES standard error, %d batches: %.4f\n", B, std(ES_batch))
@printf("VaR mean over batches: %.4f (true %.4f)\n", mean(VaR_batch), trueVaR_jl)VaR standard error, analytical: 0.0373
VaR standard error, 200 batches: 0.0386
ES standard error, 200 batches: 0.0481
VaR mean over batches: 2.3274 (true 2.3263)The two VaR standard errors should land close to each other, and do. They will not coincide. The batch figure is itself an estimate, with a standard error of roughly \(1/\sqrt{2B}\) of its own size, so a gap of several percent between the two is ordinary sampling variation rather than a sign that something is wrong. They agree at all because this example meets every condition the formula requires — independent normal draws, a continuous distribution, and a positive density at the 1% quantile. That is the point of running both. Having seen them agree where the formula is valid, we know to reach for batches, and not the formula, in the many settings where it is not.
The ES standard error is the larger of the two. VaR depends on one order statistic, while ES averages several tail observations and so inherits the variability of each. Reporting a VaR to four decimal places from a simulation whose standard error sits in the second is a common way to imply precision that is not there.
25.9 Extensions
The same structure extends to many assets, to multiple options per stock with different strikes and maturities, to realistic position sizes and to the GARCH models of Chapter 23.
25.10 Exercise
- Make a function that allows an arbitrary number of assets
- Reimplement the examples above using historical simulation
- Use GARCH(1,1) to obtain time-varying volatilities and implement backtesting