source("common/functions.r",chdir=TRUE)19 Simulation methods for risk
The risk forecasting we have done so far in these notes is for basic assets, such as equities, foreign exchange and commodities. In that case, all we need to do is to estimate the stochastic processes governing returns and use that to forecast risk.
We cannot do that for bonds or options, assets whose value derives from others. For example, bonds depend on the yield curve and options on the underlying assets, such as equities. The reason we cannot just estimate the stochastic process of such assets is that their intrinsic characteristics continually change as they move towards exploration. An option today might have a 30-day maturity but only 29 tomorrow, which means that it is converging to a fixed value so that the risk in it is decreasing. While we can use analytical methods to approximate the risk, as we do in Chapter 6 of Financial Risk Forecasting, it is usually better to use simulation methods.
The basic idea is that we have some estimate of the stochastic process of returns and use that to simulate prices one day into the future. We then can use today’s price and the Black-Scholes equation to get today’s portfolio value and then apply that to the simulated prices tomorrow. The vector of simulated profit and loss, P/L, is then the difference between the two, and we can read off the risk estimates directly from the P/L vector just like we did with historical simulation in Section Section 17.6.
19.1 Libraries
library(mvtnorm)19.2 Black-Sholes equation
bs = function(K, P, r, sigma, T){
d1 = (log(P/K) + (r + 0.5*sigma^2)*(T))/(sigma*sqrt(T))
d2 = d1 - sigma*sqrt(T)
Call = P*pnorm(d1) - K*exp(-r*(T))*pnorm(d2)
Put = K*exp(-r*(T))*pnorm(-d2) - P*pnorm(-d1)
return(list(Call=Call,Put=Put))
}19.3 Simulating returns
It is important to recognise that we are not pricing options and, therefore, don’t have to worry about simulating prices under the risk-neutral distribution, as we did in Section 10.6. We can use any stochastic process for returns we want and whichever definition of returns is most convenient.
In most of the below, we assume normality and use simple returns, but he could just as easily use any of the more complex models. It is important to recognise that we are not pricing options and, therefore, don’t have to worry about simulating prices under the risk-neutral distribution, as we did in Section Section 10.6. We can use any stochastic process for returns we want and whichever definition of returns is most convenient.
In most of the below, we assume normality and use simple returns. Still, whiche could just as easily use any of the more complex models from Chapter 17, including those with nonnormal conditional distributions and even historical simulation.
Simple returns are defined by \[ R_{t+1} = \frac{P_{t+1}}{P_t}-1 \] We know \(P_t\) and simulate \(R_{t+1}\) to get simulated price in the future \(P_{t+1}\) so \[P_{t+1,s} = P_t \times (1+R_{t+1,s})\] where \(s\) indicates a particular simulaiton.
Note \(P_{t+1}\) is not the simulated futures price.
19.4 Number of assets
In the examples below, we assume one holds one unit of each asset, which simplifies the implementation. In that case, if one holds one stock, the current portfolio value is P.
19.5 One asset — Univariate
19.5.1 Setup
par = list()
par$probability = 0.05
par$S = 1e3
par$P = 100
par$r = 0.05
par$K = 99
par$Mat = 1
par$sigma = 0.01 19.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 returns. If necessary, we could modify it to also return the Returns. It would also be straightforward to modify it to include other distributions besides the normal and even more sophisticated mean.
While S is already inside par, in some cases, we may want to input our own. For that purpose, we have a default argument, S=NULL, which does nothing. However, we can optionally input our own, which is then set by if(!is.null(S)) par$S=S. Also, we always set the seed, but if you don’t want to, we can pass NULL as the seed.
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
)
Prices = par$P*(1+Returns)
return(Prices)
}Here is one example.
plot(Sim.Prices(par,S=100),
pch=16,
col='red',
las=1,
bty='l'
)
hist(
Sim.Prices(par,S=10000),
freq=FALSE,
las=1
)
19.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 simply a quantile on the PL vector, while ES is the average from that quantile to the minimum.
19.5.4 VaR for one stock
cat("Number of simulations:",par$S,"\n")Number of simulations: 1000 today = par$P
tomorrow = Sim.Prices(par,seed=8888)
PL = tomorrow - today
VaR = -sort(PL)[par$probability*par$S]
cat("VaR = $",VaR,"\n",sep="")VaR = $1.808545tomorrow = Sim.Prices(par,seed=666)
PL = tomorrow - today
VaR = -sort(PL)[par$probability*par$S]
cat("VaR = $",VaR,"\n",sep="")VaR = $1.682279We did two simulations with different seeds, and the two VaR numbers are quite different. We analyse the number of simulations below.
19.5.5 VaR 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=sqrt(250)*par$sigma,
T=par$Mat
)[[type]]
tomorrow = bs(
K=par$K,
P=x,
r=par$r,
sigma=sqrt(250)*par$sigma,
T=par$Mat-1/365
)[[type]]
PL = tomorrow - today
VaR=-sort(PL)[par$probability*par$S]
cat("VaR = $",VaR,"\n",sep="")VaR = $1.09491919.5.6 VaR for one stock and option
When we have one asset, and both hold it and have an option on it, we only simulate the prices once and use the same simulated price vector for both.
x=Sim.Prices(par)
type = "Call"
today = par$P
today =
today+
bs(
K=par$K,
P=par$P,
r=par$r,
sigma=sqrt(250)*par$sigma,
T=par$Mat
)[[type]]
tomorrow = x
tomorrow =
tomorrow +
bs(
K=par$K,
P=x,
r=par$r,
sigma=sqrt(250)*par$sigma,
T=par$Mat-1/365
)[[type]]
PL= tomorrow - today
VaR=-sort(PL)[par$probability*par$S]
cat("VaR = $",VaR,"\n",sep="")VaR = $2.73594719.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 call rmvnorm(), but it would be straightforward to do it manually with a Choleski decomposition.
19.6.1 Setup
par = list()
par$probability = 0.05
par$S = 1e3
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)19.6.2 Simulate bivariate normal
It is then easy to use the number of simulations and the covariance matrix, par$Sigma, to simulate correlated returns and hence prices. Note that par$Sigma is the covariance matrix, while par$sigma is the standard deviation.
print(par$Sigma) [,1] [,2]
[1,] 0.010 0.005
[2,] 0.005 0.020r=rmvnorm(10,sigma=par$Sigma)
print(r) [,1] [,2]
[1,] 0.008934439 -0.15707369
[2,] -0.089471235 0.07925927
[3,] 0.076655490 -0.13516996
[4,] -0.116082299 -0.10491110
[5,] -0.096791561 0.02192440
[6,] 0.016549428 -0.12714536
[7,] -0.123844029 0.05159013
[8,] -0.086154555 -0.14738669
[9,] -0.152007705 -0.12858155
[10,] -0.004291241 -0.22524615print(cov(r)) [,1] [,2]
[1,] 0.005652299 -0.003724937
[2,] -0.003724937 0.010262649print(cor(r)) [,1] [,2]
[1,] 1.0000000 -0.4890764
[2,] -0.4890764 1.000000019.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))
}19.6.4 Two stocks
today = sum(par$P)
tomorrow = rowSums(Sim.BivarPrices(par))
PL = tomorrow - today
VaR = -sort(PL)[par$probability*par$S]
cat("VaR = $",VaR,"\n",sep="") VaR = $19.2699319.6.5 Two stocks and two options
19.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]),
T=par$Mat[i]
)[[type[i]]]
}19.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]),
T=par$Mat[i]-1/365
)[[type[i]]]
}
PL = tomorrow - today
VaR=-sort(PL)[par$probability*par$S]
cat("VaR = $",VaR,"\n",sep="") VaR = $30.5895219.7 Sim VaR vs. true VaR
Since we know the true VaR in cases where we don’t have the option, we can easily compare the true number to the one obtained by simulations. As the number of simulations increases, these two values should converge. And if they don’t, we have done something wrong. We can also use this to ascertain how many simulations we need.
19.7.1 Setup
par = list()
par$probability = 0.01
par$P = 100
par$sigma = 0.01 19.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 number of simulations, we can get an idea of the accuracy.
simVaR=function(par,S,seed=14){
set.seed(seed)
today = par$P
par$S = S
tomorrow = Sim.Prices(par,seed=seed)
PL = tomorrow - par$P
VaR = -sort(PL)[par$probability*par$S]
return(VaR)
}Here, we show the simulated and true VaR.
VaR = simVaR(par=par,S=100)
cat("VaR = $",VaR,"\n",sep="")VaR = $2.136976trueVaR = - par$P * qnorm(par$probability) * par$sigma
cat("trueVaR = $",trueVaR,"\n",sep="")trueVaR = $2.326348We can repeat this:
trueVaR = - par$P * qnorm(par$probability) * par$sigma
N=10
VaR=vector(length=10)
for(i in 1:N) VaR[i] = simVaR(par=par,S=100,seed=i)
Error=VaR- trueVaR
True=0*VaR +trueVaR
df=as.data.frame(cbind(True,VaR,Error))
print(df) True VaR Error
1 2.326348 2.214700 -0.1116480
2 2.326348 2.451706 0.1253585
3 2.326348 2.265401 -0.0609468
4 2.326348 1.797382 -0.5289659
5 2.326348 2.183967 -0.1423811
6 2.326348 1.952349 -0.3739992
7 2.326348 1.785893 -0.5404544
8 2.326348 3.014527 0.6881794
9 2.326348 2.617706 0.2913577
10 2.326348 2.185287 -0.1410610We can also vary the number of simulations.
S=c(1e3,1e4,1e5,1e6,1e7)
VaR=S
for(i in 1:length(S)) VaR[i] = simVaR(par=par,S=S[i])
df=as.data.frame(cbind(S,VaR))
df$error=df$VaR - trueVaR
print(df) S VaR error
1 1e+03 2.327881 0.0015331722
2 1e+04 2.392073 0.0657246966
3 1e+05 2.315741 -0.0106066797
4 1e+06 2.325955 -0.0003930874
5 1e+07 2.326765 0.0004172930S=c(seq(1e4,1e4,by=100),seq(1e4,1e6,by=1000))
VaR=S
for(i in 1:length(S)) VaR[i] = simVaR(par=par,S=S[i])
df=as.data.frame(cbind(S,VaR))
df$error=df$VaR - trueVaR
plot(
df[,1],df[,2],
type='l',
lty=1,
col="red",
bty='l',
las=1,
xlab="simulations",
ylab="VaR")
segments(min(S),trueVaR,max(S),trueVaR,col="blue",lwd=3)
19.8 Extensions
It would be straightforward to extend the examples above to allow for the following:
- A variable number of assets
- Holding multiple options on each stock, with puts and calls, a range of strikes and maturities
- Variable number of units of its assets held
19.9 Exercise
- Make a function that allows an arbitrarily sized number of assets
- Reimplement the examples above using historical simulation
- Use GARCH(1,1) to obtain time-varying volatilities and implement backtesting