【代码实践】使用Garch模型估计VaR

---

title: “Value at Risk estimation using GARCH model”
author: “Ionas Kelepouris & Dimos Kelepouris”
date: “July 6, 2019”
output:
html_document:
fig_height: 7
fig_width: 10
highlight: tango
toc: yes

knitr::opts_chunk$set(echo = TRUE , warning = FALSE , error = FALSE , message = FALSE)

Introduction

本分析的目的：通过给定的时变波动率构建VaR。VaR被广泛地运用于计算市场风险。
我们的时间序列包括2668天的黄金期货收益。为了解释日度收益的部分波动率，我们
使用Box-Jenkins方法拟合ARIMA模型，并测试其假设。
其次，我们检测了收益的概率密度分布是否符合正态分布，以先择最合适的分布形式。
再次，我们使用Garch模型估算了残差的条件方差，并与delta-normal方法及进行对比。
最后，我们进行了1-step的VaR预测，并且回测数据以检验模型是否完善。

Data & Libraries

为了建模，我们收集了10年的黄金期货日度收益，共计2668个观测值。

# Load libraries
library(tidyverse)
library(ggthemes)
library(forecast)
library(tseries)
library(gridExtra)
library(rugarch)

# Load data
stocks = read.csv('E:\\BaiduNetdiskWorkspace\\Gold_CSS\\Model\\R.csv' , header = T)
stocks = stocks %>% select(Date , R , Css_e)

rets = stocks$R

p1 = qplot(x = 1:length(rets) , y = rets , geom = 'line') + geom_line(color = 'darkblue') + 
    geom_hline(yintercept = mean(rets) , color = 'red' , size = 1) + 
    labs(x = '' , y = 'Daily Returns')

p2 = qplot(rets , geom = 'density') + coord_flip() + geom_vline(xintercept = mean(rets) , color = 'red' , size = 1) +
    geom_density(fill = 'lightblue' , alpha = 0.4) + labs(x = '')

grid.arrange(p1 , p2 , ncol = 2)

为了识别收益的平稳性，我们使用Augmented Dickey-Fuller检验。
其零假设为：不平稳。

adf.test(rets)

Small P-value (<0.01) 表明可以拒绝零假设。我们的收益率是平稳的。

Box-Jenkins Methodology

Box-Jenkins approach可应用ARIMA模型以发现时间序列最好的拟合，该方法主要有三个
阶段：
a) identification, b) estimation, c) diagnostic checking.

Identification

前期我们已经检验时间序列的平稳性。基于Autocorelation Function (ACF) and Partial Autocorrelation Function (PACF)，以选择合适的参数p,d,q。Akaike Information Criterion (AICc)同样可以用于选取参数。

$AIC = ln\frac{\sum\hat{u}^2}{T} + \frac{2k}{T}$

where

• $\sum {\hat{u}}^2$ = Sum of Squared Residuals
• $T$ = number of observations
• $k$ = number of model parameters (p + q + 1)

AIC可以解决过度拟合以及拟合不足的风险。AIC值最低的模型参数将被选择。

model.arima = auto.arima(rets , max.order = c(3 , 0 ,3) , stationary = TRUE , trace = T , ic = 'aicc')

根据上述结果，选择ARIMA(0,0,2)

Estimation

可以估计所选参数的coefficients，我们使用Maximum Likelihood。

model.arima

Therefore the process can be described as:

$r_{t} = -0.1238*\epsilon_{t-1} - 0.1611*\epsilon_{t-2} + \epsilon_{t}$ where $\epsilon_{t}$ is White Noise

Diagnostics Checking

该过程包含observing residual plot以及ACF & PACF diagra,以及Ljung-Box的检验
结果。如果ACF & PACF表明模型残差不存在显著滞后项，说明模型合适。

model.arima$residuals %>% ggtsdisplay(plot.type = 'hist' , lag.max = 14)

ACF 和 PACF结果相似，其相关性接近于0。右侧残差接近正态分布 N(0 , ${\sigma }^2$).

为了进一步验证残差不相关，we perform Ljung-Box test.

$$Q_{LB} = T(T+2)\sum_{s=1}^{m} \frac{\hat\rho_{s}^{2}}{T-s}$$

The $Q_{LB}$ statistic follows asymmetrically a $X^2$ distribution with m-p-q degrees of freedom. The null hypothesis refers to $H_0:{\rho }_1={\rho }_2=\cdots ={\rho }_m=0$

ar.res = model.arima$residuals
Box.test(model.arima$residuals , lag = 14 , fitdf = 2 , type = 'Ljung-Box')

我们不能拒绝零假设，因此残差的过程就像白噪声一样，所以没有可能被建模的模式的迹象。

GARCH Implementation

尽管ACF & PACF of residuals都表明没有显著的滞后性，但时间序列图仍然表明存在波动聚集。ARIMA对数据线性建模，无法反应新信息。为了对波动率建模，我们使用ARCH模型。
ARCH是一个时间序列数据的统计模型，它描述了当前误差项的方差作为前一时间段误差项的实际大小的函数。

We assume that the time series of interest, $r_t$, is decomposed into two parts, the predictable and unpredictable component, $r_t=E(r_t|I_{t-1})+{ϵ}_t$, where $I_{t-1}$ is the information set at time $t-1$ and $E(r_t|I_{t-1})=0.0437\ast r_{t-1}-0.0542\ast r_{t-2}$ and ${ϵ}_t$ is the unpredictable part, or innovation process.

The unpredictable component, can be expressed as a GARCH process in the following form:

$${ϵ}_t=z_t\ast {\sigma }_t$$

where $z_t$ is a sequence of independently and identically distributed random variables with zero mean and variance equal to 1. The conditional variance of ${ϵ}_t$ is ${\sigma }_t$, a time-varying function of the information set at time $t-1$.

Next step is to define the the second part of the error term decomposition which is the conditional variance, ${\sigma }_t$. For such a task, we can use a GARCH(1 , 1) model, expressed as:

$${\sigma }_t^2=\omega +a_1\ast {ϵ}_{t-1}^2+{\beta }_1\ast {\sigma }_{t-1}^2$$

当残差平方存在自相关时，可以使用Garch模型。ACF and PACF plots clearly indicate significant correlation.

tsdisplay(ar.res^2 , main = 'Squared Residuals')

另一种检测残差平方的异方差性的方法测试参数的显著性。
$a_1$ and ${\beta }_1$ parameters.

# Model specification
model.spec = ugarchspec(variance.model = list(model = 'sGARCH' , garchOrder = c(1 , 1)) , 
                        mean.model = list(armaOrder = c(0 , 2)))

model.fit = ugarchfit(spec = model.spec , data = ar.res , solver = 'solnp')

options(scipen = 999)
model.fit@fit$matcoef

Both $a_1$ and ${\beta }_1$ are significantly different from zero,
因此可以假设残差存在波动聚集现象。

With successive replacement of the ${\sigma }_{t-1}^2$ term, the GARCH equation can be written as:

$${\sigma }_t^2=\frac{\omega }{1-{\beta }_1}+a_1\ast \sum _{i=1}^{\infty }{\beta }_1^{i-1}\ast {ϵ}_{t-i}^2$$

When we replace with the coefficient estimates given by optimization we get the following equation:

$${\sigma }_t^2=0.000087+0.108\ast ({ϵ}_{t-1}^2+0.825\ast {ϵ}_{t-2}^2+0.680\ast {ϵ}_{t-3}^2+0.561\ast {ϵ}_{t-4}^2+\dots )$$

Given that $0<{\beta }_1<1$, as lag increases the effect of the squared residual decreases.

Value at Risk

VaR用于测度当下位置的下跌风险。
A VaR statistic has three components: a) time period, b) confidence level, c) loss ammount (or loss percentage). For 95% confidence level, we can say that the worst daily loss will not exceed VaR estimation. If we use historical data, we can estimate VaR by taking the 5% quantile value. For our data this estimation is:

quantile(rets , 0.05)
qplot(rets , geom = 'histogram') + geom_histogram(fill = 'lightblue' , bins = 30) +
    geom_histogram(aes(rets[rets < quantile(rets , 0.05)]) , fill = 'red' , bins = 30) +
    labs(x = 'Daily Returns')

Red bars refer to returns lower than 5% quantile.

Distributional Properties

为了估算VaR,我们需要适当定义假设分布的quantile。
对于正态分布而言，5%对应的quantile为-1.645。实证结果表明正态分布的假设往往
产生weak results。Jarque-Bera test can test the hypothesis that stock returns follow a normal distribution.

$$JB=\frac{n-k+1}{6}\ast (S^2+\frac{1}{4}\ast (C-3{)}^2)$$

where $S$ is skewness and $C$ is kurtosis. A normal distributed sample would return a JB score of zero. The low p-value indicates stock returns are not normally distributed.

jarque.bera.test(rets)

p2_1 = qplot(rets , geom = 'density') + geom_density(fill = 'blue' , alpha = 0.4) + 
    geom_density(aes(rnorm(200000 , 0 , sd(rets))) , fill = 'red' , alpha = 0.25) + 
    labs(x = '')

p2_2 = qplot(rets , geom = 'density') + geom_density(fill = 'blue' , alpha = 0.4) + 
    geom_density(aes(rnorm(200000 , 0 , sd(rets))) , fill = 'red' , alpha = 0.25) + 
    coord_cartesian(xlim = c(-0.07 , -0.02) , ylim = c(0 , 10)) + 
    geom_vline(xintercept = c(qnorm(p = c(0.01 , 0.05) , mean = mean(rets) , sd = sd(rets))) , 
               color = c('darkgreen' , 'green') , size = 1) + labs(x = 'Daily Returns')

grid.arrange(p2_1 , p2_2 , ncol = 1)

如图所示，Density plots are shown for 收益（蓝）以及正态分布（红）数据。
Vertical lines表明了5%和1%对应的quantile。
正态分布都会高估风险。

Student’s t-distribution

为了更好对收益建模，我们使用其他分布。
The t-distribution is symmetric and bell-shaped, like the normal distribution, but has heavier tails, meaning that it is more prone to producing values that fall far from its mean. We use the fitdist function from rugarch package to get the fitting parameters of t-distribution.

fitdist(distribution = 'std' , x = rets)$pars
cat("For a = 0.05 the quantile value of normal distribution is: " , 
    qnorm(p = 0.05) , "\n" ,
     "For a = 0.05 the quantile value of t-distribution is: " ,
    qdist(distribution = 'std' , shape = 3.7545967917 , p = 0.05) , "\n" , "\n" , 
    'For a = 0.01 the quantile value of normal distribution is: ' , 
    qnorm(p = 0.01) , "\n" , 
    "For a = 0.01 the quantile value of t-distribution is: " , 
    qdist(distribution = 'std' , shape = 3.7545967917 , p = 0.01) , sep = "")

95%置信区间下，正态分布高估风险；
99%置信区间下，正态分布未能捕捉outliers，低估风险。

Garch VaR vs Delta-normal approach

Delta-normal approach 假设收益正态分布。

$$VaR(a)=\mu +\sigma \ast N^{-1}(a)$$

where $\mu$ is the mean stock return, $\sigma$ is the standard deviation of returns, $a$ is the selected confidence level and $N^{-1}$ is the inverse PDF function, generating the corresponding quantile of a normal distribution given $a$.

The results of such a simple model is often disapointing and are rarely used in practice today. The assumption of normality and constant daily variance is usually wrong and that is the case for our data as well.

Previously we observed that returns exhibit time-varying volatility. Hence for the estimation of VaR we use the conditional variance given by GARCH(1,1) model. For the underlined asset’s distribution properties we use the student’s t-distribution. For this method Value at Risk is expressed as:

$$VaR(a)=\mu +{\hat{\sigma }}_{t|t-1}\ast F^{-1}(a)$$

where ${\hat{\sigma }}_{t|t-1}$ is the conditional standard deviation given the information at $t-1$ and $F^{-1}$ is the inverse PDF function of t-distribution.

qplot(y = rets , x = 1:2668 , geom = 'point') + geom_point(colour = 'lightgrey' , size = 2) + 
    geom_line(aes(y = model.fit@fit$sigma*(-1.485151) , x = 1:2668) , colour = 'red') +
    geom_hline(yintercept = sd(rets)*qnorm(0.05) , colour = 'darkgreen' , size = 1.2) + theme_light() + 
    labs(x = '' , y = 'Daily Returns' , title = 'Value at Risk Comparison')

Red line denotes VaR produced by GARCH model and green line refers to delta-normal VaR.

VaR forecasting

The ugarchroll method 可以进行滚动预测，it returns分布预测参数以计算forecasted density的必要测度。我们将最后500个观测设为测试集，对条件标准差${\hat{\sigma }}_{t+1|t}$进行1-step的预测。

model.roll = ugarchroll(spec = model.spec , data = rets , n.start = 2168 , refit.every = 50 ,
                        refit.window = 'moving')

# Test set 500 observations
VaR95_td = mean(rets) + model.roll@forecast$density[,'Sigma']*qdist(distribution='std', shape=3.7545967917, p=0.05)

Backtesting

设 $N={\sum }_{t=1}^T{I}_t$ 为收益低于VaR估计的个数。

VaR estimate, where $I_t$ is 1 if $y_{t+1}<Va{R}_{t+1|t}$ and 0 if $y_{t+1}\ge Va{R}_{t+1|t}$. Hence, N is the observed number of exceptions in the sample. As argued in Kupiec (1995), the failure number follows a binomial distribution, B(T, p).

p = c()
p[1] = pbinom(q = 0 , size = 500 , prob = 0.05)
for(i in 1:50){
    p[i] = (pbinom(q = (i-1) , size = 500 , prob = 0.05) - pbinom(q = (i-2) , size = 500 , prob = 0.05))
}
qplot(y = p , x = 1:50 , geom = 'line') + scale_x_continuous(breaks = seq(0 , 50 , 2)) + 
    annotate('segment' , x = c(16 , 35) , xend = c(16 , 35) , y = c(0 , 0) , yend = p[c(16 , 35)] , color = 'red' , 
             size = 1) + labs(y = 'Probability' , x = 'Number of Exceptions') + theme_light()

该图代表了异常的分布概率。The expected number为25个，两条红线表示95%的置信区间。 the lower being 16 and the upper 35.
我们期待16-35之间的数字，to state that GARCH model as successfully preditive.

qplot(y = VaR95_td , x = 1:500 , geom = 'line') +
    geom_point(aes(x = 1:500 , y = rets[2169:2668] , color = as.factor(rets[2169:2668] < VaR95_td)) , size = 2) + scale_color_manual(values = c('gray' , 'red')) + 
    labs(y = 'Daily Returns' , x = 'Test set Observation') + theme_light() + 
    theme(legend.position = 'none')

Black line为由Garch模型所预测VaR，red points refer to returns lower than VaR.
Final step is to count the number of exceptions and compare it with the one generated with delta-normal approach.

cat('Number of exceptions with delta-normal approach: ' , (sum(rets[2169:2668] < (mean(rets) + qnorm(p = 0.05)*sd(rets[1:2168])))) , '\n' , 'Number of exceptions with GARCH approach: ' , (sum(rets[2169:2668] < VaR95_td)) , sep = '')

正如前述，we expected that delta-normal approach会高估风险。
当回测时，只有12个收益低于VaR预测的95%下界（16），而 GARCH approach (23 exceptions) seems to be an effective predictive tool in this particular case.

References

Angelidis T., Benos A. and Degiannakis S. (December 2003). The Use of GARCH Models in VaR Estimation.

Ghalanos A. (August 2017). Introduction to the rugarch package (Version 1.3-8).

Montgomery D., Jennings C. and Kulahci M. (2015). Introduction to Time Series Analysis and Forecasting (Second Edition). New Jersey: Wiley.