Mathematics and Finance Summer School at CMI
Sourish Das
Chennai Mathematical Institute25 May 2018 to 25 June, 2018
Risk Analysis
Volatility Risk, GARCH Model with R and Value at Risk
We consider the log-return of an asset. That is, \[ r_t=\log(P_t)-\log(P_{t-1}), \] where \(P_t\) is the price of the asset at time \(t\).
The mean return \[ \bar{r}=\frac{1}{m}\sum_{t=1}^m r_{n-t}, \] where \(m\) is the number of observations leading upto present period.
If we assume that the mean return is zero, we obtain the maximum likelihood estimator of variance: \[ \sigma_n^2=\frac{1}{m}\sum_{t=1}^m r_{n-t}^2. \]
This formula to estimate volatility assign weght equally to each observations.
So more remote older return have the same influence on the estimated volatility, as returns that are more current.
As our aim is to estimate the present level of volatility, we would like to assign weight on current return more heavily than the older one.
The generic scheme can be presented as \[ \sigma_n^2=\sum_{t=1}^m \alpha_t r_{n-t}^2, ~~ \sum_{t=1}^m \alpha_t =1, \] \(\alpha_t\) is the weight of the return \(t\) days ago.
As our aim is to have a greater influence on current returns, then weight should decline in value for older returns, i.e., \(\alpha_t \geq \alpha_{t+1}\)
An addition to this representation is to assume a long-run variance.
The most commonly used model is the autoregressive conditional heteroskedasticity model, ARCH(m), which can be represented by: \[ \sigma_n^2=\gamma \sigma_L+ \sum_{t=1}^m \alpha_t r_{n-t}^2,~~ ~~ \gamma+ \sum_{t=1}^m \alpha_t = 1, \] where \(\sigma_L\) is the long-run variance weighted by parameter \(\gamma\).
We consider ARCH(1) model, which can be presented as \[ \sigma_n^2=\gamma \sigma_L+ \alpha_1 r_{n-1}^2,~~ ~~ \gamma+ \alpha_1 = 1, \]
We replace \(\alpha_1\) by \((1-\gamma)\) we have the model as \[ \sigma_n^2=\gamma \sigma_L+ (1-\gamma) r_{n-1}^2, \] where \(\gamma\) is weight on long-run variance.
Though ARCH model consider long-run variance, however it does not consider the previous volatility estimate.
This disadvantage of ARCH model can be overcome by the generalized autoregressive conditional heteroskedasticity, more popularly known as the GARCH model.
The GARCH(1,1) model can be presented as as \[ \sigma_n^2=\gamma \sigma_L+ \alpha_1 r_{n-1}^2 + \beta_1 \sigma_{n-1}^2,~~ ~~ \gamma+ \alpha_1 +\beta_1 = 1, \] where \(\alpha_1\) is the weight on previous period’s return, \(\beta_1\) is the weight on previous volatility estimate, \(\gamma\) is the weight on long-run variance.
We can write \(\omega=\gamma \sigma_L^2\), such that \[ \sigma_n^2=\omega+ \alpha_1 r_{n-1}^2 + \beta_1 \sigma_{n-1}^2, \] where \(\sigma_L^2=\frac{\omega}{1-\alpha_1-\beta_1}\) is the long-run variances.
If \(\omega=0\), \(\alpha=(1-\lambda)\) and \(\beta_1=\lambda\) then we have Exponentially Weighted Moving Average (EWMA) model, i.e., \[ \sigma_n^2=(1-\lambda)r_{n-1}^2 + \lambda \sigma_{n-1}^2 \] where \(\lambda\) is exponential decay rate. So EWMA is the special case of GARCH model.
Fitting GARCH Model
library(tseries)
VIX<-get.hist.quote(instrument = "^VIX"
,start="2015-01-01"
,end=Sys.Date()
,quote="AdjClose"
,provider = "yahoo")## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.##
## WARNING: There have been significant changes to Yahoo Finance data.
## Please see the Warning section of '?getSymbols.yahoo' for details.
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.yahoo.warning"=FALSE).## time series starts 2015-01-02
## time series ends 2018-05-23par(mfrow=c(1,2))
plot(ts(VIX$Adjusted),ylab="VIX")
## Why are we using SPY and not ^GSPC?
SnP.500<-get.hist.quote(instrument = "SPY"
,start="2015-01-01"
,end=Sys.Date()
,quote="AdjClose"
,provider = "yahoo")## time series starts 2015-01-02
## time series ends 2018-05-23plot(ts(SnP.500$Adjusted),ylab="S&P 500")
## log-return
ln_rt<-as.numeric(diff(log(SnP.500)))
## Fit GARCH(1,1)
fit.garch<-garch(ln_rt,order=c(1,1))##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 6.132124e-05 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 -3.688e+03
## 1 7 -3.689e+03 4.29e-04 1.76e-03 1.0e-04 6.5e+10 1.0e-05 5.70e+07
## 2 8 -3.690e+03 1.90e-04 2.57e-04 9.1e-05 2.0e+00 1.0e-05 2.08e+01
## 3 9 -3.690e+03 1.08e-05 9.74e-06 9.9e-05 2.0e+00 1.0e-05 2.11e+01
## 4 17 -3.714e+03 6.46e-03 1.15e-02 5.2e-01 2.0e+00 1.1e-01 2.10e+01
## 5 20 -3.747e+03 8.72e-03 7.55e-03 7.0e-01 1.9e+00 2.8e-01 1.03e+00
## 6 28 -3.748e+03 2.97e-04 8.68e-04 4.5e-06 1.1e+01 3.1e-06 3.42e+00
## 7 29 -3.748e+03 2.83e-05 2.40e-05 4.3e-06 2.0e+00 3.1e-06 4.81e+00
## 8 30 -3.748e+03 8.34e-07 8.39e-07 4.4e-06 2.0e+00 3.1e-06 4.93e+00
## 9 40 -3.772e+03 6.28e-03 6.30e-03 2.3e-01 2.0e+00 2.0e-01 4.90e+00
## 10 42 -3.777e+03 1.40e-03 1.40e-03 3.6e-02 1.9e+00 4.0e-02 4.34e-02
## 11 44 -3.781e+03 1.08e-03 1.16e-03 3.3e-02 1.8e+00 4.0e-02 5.87e-02
## 12 45 -3.786e+03 1.27e-03 1.67e-03 6.1e-02 9.4e-01 8.0e-02 3.83e-03
## 13 47 -3.786e+03 1.82e-05 1.55e-04 8.3e-03 1.4e+00 1.2e-02 1.87e-04
## 14 48 -3.786e+03 9.68e-05 8.82e-05 7.6e-03 9.4e-01 1.2e-02 9.28e-05
## 15 49 -3.787e+03 5.67e-05 4.86e-05 7.7e-03 0.0e+00 1.2e-02 4.86e-05
## 16 51 -3.787e+03 3.81e-05 2.64e-05 7.9e-03 0.0e+00 1.6e-02 2.64e-05
## 17 52 -3.787e+03 1.41e-05 1.63e-05 5.9e-03 0.0e+00 1.2e-02 1.63e-05
## 18 53 -3.787e+03 5.82e-07 7.19e-07 5.4e-04 0.0e+00 1.1e-03 7.19e-07
## 19 54 -3.787e+03 3.92e-08 4.96e-08 3.2e-04 0.0e+00 4.8e-04 4.96e-08
## 20 55 -3.787e+03 7.50e-10 4.46e-10 3.5e-05 0.0e+00 6.1e-05 4.46e-10
## 21 56 -3.787e+03 1.85e-11 2.21e-12 1.5e-06 0.0e+00 2.8e-06 2.21e-12
## 22 57 -3.787e+03 -1.74e-12 2.13e-15 5.0e-08 0.0e+00 7.6e-08 2.13e-15
##
## ***** RELATIVE FUNCTION CONVERGENCE *****
##
## FUNCTION -3.786816e+03 RELDX 5.018e-08
## FUNC. EVALS 57 GRAD. EVALS 22
## PRELDF 2.132e-15 NPRELDF 2.132e-15
##
## I FINAL X(I) D(I) G(I)
##
## 1 4.276370e-06 1.000e+00 -5.952e+00
## 2 1.910014e-01 1.000e+00 1.083e-06
## 3 7.457199e-01 1.000e+00 -8.096e-05## Fitted Values
sigma_hat<-fit.garch$fitted.values[,"sigt"]
vol_hat<-sigma_hat*sqrt(252)*100
## plot GARCH fitted volatility vs. VIX
par(mfrow=c(1,1))
plot(ts(VIX$Adjusted),col="red",lwd=2,ylim=c(0,40))
lines(vol_hat,col="blue",lwd=2)
text<-c("VIX","GARCH(1,1)")
legend(350,40,text,col=c("red","blue"),lwd=c(2,2))
Value at Risk (VaR)
If Volatility is a measure of average risk or (average loss) of investment.
The Value at Risk (VaR) is a measure of the etreme risk of investments.
VaR is defined as: for a given portfolio, time horizon, and probability \(p\), the \(p\) VaR is defined as a threshold loss value, such that the probability that the loss on the portfolio over the given time horizon exceeds this value is \(p\).
For example, if a portfolio of stocks has a one-day 5% VaR of $100,000, that means that there is a 0.05 probability that the portfolio will fall in value by more than $100,000 over a one-day period.
If \(X\) is the underlying price of the portfolio, then \(VaR_{\alpha}(X)\) is the negative if the \(\alpha\)-quantile, i.e., \[ VaR_{\alpha}(X)=\inf\{x \in \mathbb{R}: \mathbb{P}(X+x<0)\leq 1-\alpha\}=\inf\{x\in \mathbb{R}: 1-F_X(-x)\geq \alpha\} \]
Value at Risk if log return follow Gaussian distribution \[ VaR_{\alpha}(X)= |V*Z_{1-\alpha}*\sigma*\sqrt{t}|, \] where \(V\) is portfolio value, \(Z_{1-\alpha}\) is the (\(1-\alpha\)) quantile, \(\sigma\) is the portfolio volatility and \(t\) is the holding period
# Daily volatility
vol <- sd(ln_rt)
# Daily average return
rbar <- mean(ln_rt)
###
n<-length(SnP.500)
unit <- 1000 # Number of units of S&P 500
Value<-SnP.500[n] # Current Value of S&P 500
Value## 2018-05-23
## 273.36value_p <- Value*unit # Value of portfolio
hp <- 1 # Holding period
a <- .95 # Confidence level (5%)
# Parametric VaR with Gaussian
parvar_Gauss <- abs(value_p*qnorm(1-a,0,1)*vol*sqrt(hp))
# Parametric VaR with t(15)
parvar_t <- abs(value_p*qt(1-a,df=15)*vol*sqrt(hp))
# Historical VaR
hvar <- abs(quantile(ln_rt,1-a)*value_p)
# Vector comparing the two VarS
varv <- cbind(value_p,parvar_Gauss,parvar_t,hvar)
names(varv)<- c("Portfolio Value","VaR Gaussian","VaR t(15)","Historical VaR")
print(varv)## Portfolio Value VaR Gaussian VaR t(15) Historical VaR
## 2018-05-23 273360 3711.474 3955.611 3713.435