19 改进的GARCH模型

本章讲GARCH模型的一些有针对性的改进。 来自(Tsay 2013)§4.9-4.11内容。

19.1 EGARCH模型

19.1.1 模型

(Nelson 1991)提出的指数GARCH(EGARCH)模型允许正负资产收益率对波动率有不对称的影响。 考虑如下变换 \[\begin{align} g(\varepsilon_t) = \alpha\varepsilon_t + \gamma\left[ |\varepsilon_t| - E|\varepsilon_t| \right] , \tag{19.1} \end{align}\] 其中\(\alpha\)\(\gamma\)是实常数。 \(\{\varepsilon_t \}\)\(\{ |\varepsilon_t| - E|\varepsilon_t| \}\)都分别是零均值独立同分布白噪声, 分布为连续分布。 易见\(Eg(\varepsilon_t)=0\)

由下式可见\(g(\varepsilon_t)\)的分布是非对称的: \[\begin{align} g(\varepsilon_t) = \begin{cases} (\alpha + \gamma) \varepsilon_t - \gamma E |\varepsilon_t|, & \varepsilon_t \geq 0 , \\ (\alpha - \gamma) \varepsilon_t - \gamma E |\varepsilon_t|, & \varepsilon_t < 0 . \end{cases} \tag{19.2} \end{align}\]

\(\varepsilon_t \sim \text{N}(0,1)\)\(E|\varepsilon_t| = \sqrt{\frac{2}{\pi}}\)。 对式(17.6)中的标准t分布,有 \[ E|\varepsilon_t| = \frac{2 \sqrt{v-2} \Gamma(\frac{v+1}{2})}{ (v-1)\Gamma(\frac{v}{2}) \sqrt{\pi}} . \]

EGARCH(\(m,s\))模型可以用滞后算子的形式写成

\[\begin{align} a_t = \sigma_t \varepsilon_t, \quad \ln\sigma_t^2 = \alpha_0' + \frac{1 + \alpha_2 B + \dots + \alpha_m B^{m-1}}{ 1 - \beta_1 B - \dots - \beta_s B^s} g(\varepsilon_{t-1}) \tag{19.3} \end{align}\] \(\alpha_0'\)为常数, 其中\(B\)是滞后算子, 多项式\(1 + \alpha_2 z + \dots + \alpha_m z^{m-1}\)\(1 - \beta_1 z - \dots - \beta_m z^m\)的根都在单位圆外且两个多项式没有公因子。

注意这里的模型阶相当于GARCH(\(m,s\))。

\(\xi_t = \ln\sigma_t^2\), 则(19.3)给出的\(\xi_t\)为一个平稳线性ARMA\((s, m-1)\)序列, 以零均值独立同分布白噪声\(g(\varepsilon_{t-1})\)为新息; 但是, \(\ln\sigma_t^2\)通过\(\varepsilon_{t-j} = a_{t-j} / \sigma_{t-j}\)\(\{ a_t \}\)序列依赖。 原始的GARCH模型的\(\sigma_t^2\)的方程是直接依赖于\(a_{t-j}^2\)的,\(\pm a_{t-j}\)\(\sigma_t^2\)影响相同。

易见\(E\ln\sigma_t^2 = \alpha_0'\)

更一般地, 可以令\(g(\cdot)\)中的\(\gamma\)随滞后\(j\)变化, 模型变成:

\[\begin{align} a_t =& \sigma_t \varepsilon_t, \\ \ln\sigma_t^2 =& \alpha_0 + \sum_{j=1}^m \left[\alpha_j \varepsilon_{t-j} + \gamma_j(|\varepsilon_{t-j}| - E |\varepsilon_{t-j}|) \right] + \sum_{i=1}^s \beta_i \ln \sigma_{t-i}^2 . \tag{19.4} \end{align}\]

(19.4)中, \(\alpha_j\)代表了对数收益率的正负扰动对波动率的不同影响, 如果\(\alpha_j=0\), 则正负扰动对波动率影响就相同了。

EGARCH与GARCH模型的区别还有:

  1. 使用条件方差的对数建模, 因为对数值可正可负, 这就取消了GARCH模型对系数必须非负的限制。
  2. \(g(\varepsilon_{t-j})=g(a_{t-j}/\sigma_{t-j})\)的使用使得波动率对\(a_{t-j}\)的依赖关系与\(a_{t-j}\)的正负号有关, 可以用来描述正负收益率对波动率的不同的影响, 即杠杆效应。

19.1.2 EGARCH(1,1)模型

以最简单的EGARCH(1,1)为例来讨论。 模型为

\[\begin{align} a_t =& \sigma_t \varepsilon_t, \\ \ln \sigma_t^2 =& \alpha_0 + \left[ \alpha_1 \varepsilon_{t-1} + \gamma_1 (|\varepsilon_{t-1}| - E|\varepsilon_{t-1}|) \right] + \beta_1 \ln\sigma_{t-1}^2 . \tag{19.5} \end{align}\]

其中\(\varepsilon_t\) iid N(0,1), 这个模型的阶类似GARCH(1,1), 模型实际上是关于\(\ln \sigma_t^2\)的一个AR(1)模型。

与一个GARCH(1,1)模型对比: \[\begin{aligned} \sigma_t^2 = \alpha_0 + \alpha_1 a_{t-1}^2 + \beta_1 \sigma_{t-1}^2 \end{aligned}\] 所以EGARCH模型用\(\ln\sigma_t^2\)替换了\(\sigma_t^2\), 用\(g(\varepsilon_{t-1})\)替换了\(\alpha_1 a_{t-1}^2\), 用\(\ln\sigma_{t-1}^2\)替换了\(\sigma_{t-1}^2\)。 注意\(g(\varepsilon_{t-1}) = g(a_{t-1} / \sigma_{t-1})\)

(19.5)中, 设\(\alpha_1 < 0\)\(\varepsilon_{t-1} \geq 0\)时, \(\ln\sigma_t^2\)变小, \(\varepsilon_{t-1} < 0\)时, \(\ln\sigma_t^2\)变大, 对于相同的\(|a_{t-1}|\), 负值比正值产生更大的\(\ln\sigma_t^2\)

波动率方程可以写成 \[\begin{aligned} \sigma_t^2 =\sigma_{t-1}^{2\beta_1} e^{\alpha_0''} \cdot \begin{cases} \exp\left[ (\gamma_1 + \alpha_1) \frac{|a_{t-1}|}{\sigma_{t-1}} \right], & a_{t-1} \geq 0 \\ \exp\left[ (\gamma_1 - \alpha_1) \frac{|a_{t-1}|}{\sigma_{t-1}} \right], & a_{t-1} < 0 \end{cases} \end{aligned}\] 其中\(\alpha_0'' = \alpha_0 - \gamma_1 E|\varepsilon_{t-1}|\)

系数\(\gamma_1 + \alpha_1\)\(\gamma_1 - \alpha_1\)表明当\(\alpha_1\neq 0\) 时模型对正的和负的\(a_{t-1}\)的响应是非对称的, 由于一般负的扰动造成的波动更大, 所以\(\alpha_1\)应该是负的才比较符合实际数据。

19.1.3 实例1:CRSP价值加权指数超额收益率

例子来自(Nelson 1991)的文章。 考虑CRSP价值加权指数的超额收益率日数据, 从1962年7月到1987年12月,共6408个观测值, 超额收益率是价值加权指数的收益率减去国债的月收益率, 假定一个月每天的国债收益率相等。

\(r_t\)表示超额收益率, 所用的模型为

\[\begin{aligned} r_t =& \phi_0 + \phi_1 r_{t-1} + c\sigma_{t}^2 + a_t, \quad a_t = \sigma_t \varepsilon_t,\\ \ln\sigma_t^2 =& \alpha_0 + \ln(1 + \omega N_t) + \frac{1 + \beta B}{1 - \alpha_1 B - \alpha_2 B^2} g(\varepsilon_{t-1}), \\ g(\varepsilon_{t-1}) =& \theta \varepsilon_{t-1} + \gamma (|\varepsilon_{t-1}| - E|\varepsilon_{t-1}|) . \end{aligned}\] 均值方程是带有回归自变量的AR(1)模型, 回归自变量是波动率平方, 这个回归自变量代表了风险溢价, \(c\)是风险溢价参数。 波动率方程中的回归自变量\(N_t\)是第\(t-1\)交易日到第\(t\)交易日之间没有交易的天数, 这一项解释了放假休市若干天以后增加的波动。 \(\varepsilon_t\)使用广义误差分布(GED), 波动率方程是是EGARCH(2,2)模型。 这里所用的符号与19.1.1有所不同。 \(\theta\)代表了杠杆效应。

(Nelson 1991)论文中参数估计值与标准误差列表如下:

\[ \begin{array}{crr} \text{参数} & \text{估计值} & \text{标准误差} \\ \hline \alpha_0 & -10.06 & 0.346 \\ \omega & 0.183 & 0.028 \\ \gamma & 0.156 & 0.023 \\ \alpha_1 & 1.929 & 0.015 \\ \alpha_2 & -0.929 & 0.015 \\ \beta & -0.978 & 0.006 \\ \theta & -0.118 & 0.009 \\ \phi_0 & 3.5 \times 10^{-4} & 9.9\times 10^{-5} \\ \phi_1 & 0.205 & 0.012 \\ c & -3.361 & 2.026 \\ v & 1.576 & 0.032 \end{array} \]

估计中只有风险溢价参数\(c\)是不显著, 其它参数, 包括代表了不对称性的\(\theta\),都是显著的, \(\theta\)估计为负, 这表明模型的波动率收到负的扰动影响更大。

19.1.4 实例2:IBM股票月对数收益率

考虑IBM股票月对数收益率建模, 从1967年1月到2009年12月, 561个观测值。 建立EGARCH(1,1)模型, 使用(19.5)的模型形式, 模型为 \[\begin{aligned} r_t =& \mu + a_t, \\ a_t =& \sigma_t \varepsilon_t, \\ \ln \sigma_t^2 =& \alpha_0 + \left[ \alpha_1 \varepsilon_{t-1} + \gamma_1 (|\varepsilon_{t-1}| - E|\varepsilon_{t-1}|) \right] + \beta_1 \ln\sigma_{t-1}^2 . \end{aligned}\]

读入数据:

da <- read_table(
  "m-ibmsp6709.txt",
  col_types=cols(.default=col_double(),
                  date=col_date(format="%Y%m%d")))
xts.ibm <- xts(log(1 + da[["ibm"]]), da[["date"]])
tclass(xts.ibm) <- "yearmon"
ts.ibm <- ts(log(1 + da[["ibm"]]), start=c(1967,1), frequency=12)

对数收益率的时间序列图:

plot(xts.ibm, 
  format.labels="%Y",
  xlab="Year", ylab="Log return")
IBM股票月对数收益率

图19.1: IBM股票月对数收益率

ACF图:

forecast::Acf(xts.ibm, lag.max = 36, main="")
IBM股票月对数收益率ACF

图19.2: IBM股票月对数收益率ACF

从对数收益率的ACF来看基本是白噪声。作Ljung-Box白噪声检验:

Box.test(xts.ibm, lag=12, type="Ljung")
## 
##  Box-Ljung test
## 
## data:  xts.ibm
## X-squared = 7.4042, df = 12, p-value = 0.8298

基本可以认为是(宽)白噪声。

\(a_t^2\)的ACF:

forecast::Acf(xts.ibm^2, lag.max = 36, main="")
IBM股票月对数收益率平方ACF

图19.3: IBM股票月对数收益率平方ACF

这表明有轻微的ARCH效应。

用rugarch扩展包可以估计EGARCH模型,程序如:

library(rugarch)
spec1 <- ugarchspec(
  mean.model = list(
    armaOrder=c(0,0),
    include.mean=TRUE  ),
  variance.model = list(
    model = "eGARCH", # exponential GARCH model
    garchOrder = c(1,1) ) )
mod1ru <- ugarchfit(spec = spec1, data = xts.ibm)
show(mod1ru)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : eGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.006649    0.002963   2.2442 0.024820
## omega  -0.423208    0.223673  -1.8921 0.058480
## alpha1 -0.094813    0.039373  -2.4081 0.016037
## beta1   0.920485    0.041729  22.0586 0.000000
## gamma1  0.218711    0.060802   3.5971 0.000322
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.006649    0.003073   2.1635 0.030502
## omega  -0.423208    0.308230  -1.3730 0.169743
## alpha1 -0.094813    0.051835  -1.8291 0.067382
## beta1   0.920485    0.057270  16.0728 0.000000
## gamma1  0.218711    0.061770   3.5407 0.000399
## 
## LogLikelihood : 651.634 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -2.5063
## Bayes        -2.4652
## Shibata      -2.5065
## Hannan-Quinn -2.4902
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      1.237  0.2661
## Lag[2*(p+q)+(p+q)-1][2]     1.344  0.3989
## Lag[4*(p+q)+(p+q)-1][5]     1.867  0.6501
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                   0.009633  0.9218
## Lag[2*(p+q)+(p+q)-1][5]  1.087446  0.8396
## Lag[4*(p+q)+(p+q)-1][9]  2.511467  0.8360
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]   0.09128 0.500 2.000  0.7626
## ARCH Lag[5]   1.10479 1.440 1.667  0.7021
## ARCH Lag[7]   2.20197 2.315 1.543  0.6746
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.1719
## Individual Statistics:              
## mu     0.21948
## omega  0.61756
## alpha1 0.15868
## beta1  0.61824
## gamma1 0.06386
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.28 1.47 1.88
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.1014 0.9192    
## Negative Sign Bias  0.2560 0.7980    
## Positive Sign Bias  0.1888 0.8503    
## Joint Effect        0.3726 0.9458    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     13.07       0.8350
## 2    30     22.26       0.8094
## 3    40     28.03       0.9040
## 4    50     42.53       0.7314
## 
## 
## Elapsed time : 0.07496881

估计的模型可以写成:

\[\begin{aligned} r_t =& 0.0066 + a_t, \quad a_t = \sigma_t \varepsilon_t, \quad \varepsilon_t \text{ iid N}(0,1) \\ \ln\sigma_t^2 =& -0.4232 - 0.0948 \varepsilon_{t-1} + 0.2187(|\varepsilon_{t-1}| - E|\varepsilon_{t-1}|) + 0.9205 \ln\sigma_{t-1}^2 . \end{aligned}\] 这种模型是关于\(\{ \ln\sigma_t^2 \}\)的一个AR(1)模型。 输出中参数alpha1代表了正负收益率对波动率不对称影响的杠杆效应, 在0.05水平下这一项显著, 说明杠杆效存在。

\(\varepsilon_{t-1}=\pm 2\)时, 不同正负号的\(\varepsilon_{t-1}\)引起的\(\sigma_t^2\)的变化比值为 \[ e^{4 \times 0.0948} = 1.46, \] 所以这时幅度相等的负的扰动比正的扰动额外增加46%的波动率平方。 幅度越大,额外增加的波动率也越大。

计算标准化残差\(\tilde a_t = a_t/\sigma_t\)和波动率\(\sigma_t\)

stdresi <- residuals(mod1ru, standardize=TRUE)
vol <- sigma(mod1ru)

标准化残差图:

plot(stdresi, format.labels="%Y",
  xlab="Year", ylab="", 
  main="IBM股票收益率模型标准化残差")
IBM股票收益率模型标准化残差

图19.4: IBM股票收益率模型标准化残差

估计的波动率:

plot(vol, format.labels="%Y",
  xlab="Year", ylab="", 
  main="IBM股票收益率模型波动率")
IBM股票收益率模型波动率

图19.5: IBM股票收益率模型波动率

标准化残差的ACF:

forecast::Acf(stdresi, lag.max = 36, main="")
IBM股票收益率模型标准化残差ACF

图19.6: IBM股票收益率模型标准化残差ACF

## 或:
## plot(mod1ru, which=10)

标准化残差的Ljung-Box白噪声检验:

Box.test(stdresi, lag=12, type="Ljung")
## 
##  Box-Ljung test
## 
## data:  stdresi
## X-squared = 6.894, df = 12, p-value = 0.8645

标准化残差平方的ACF:

forecast::Acf(stdresi^2, lag.max = 36, main="")
IBM股票收益率模型标准化残差平方的ACF

图19.7: IBM股票收益率模型标准化残差平方的ACF

## 或:
## plot(mod1ru, which=11)

标准化残差平方的Ljung-Box白噪声检验:

Box.test(stdresi^2, lag=12, type="Ljung")
## 
##  Box-Ljung test
## 
## data:  stdresi^2
## X-squared = 5.537, df = 12, p-value = 0.9376

从以上残差诊断看模型比较充分。 从ugarchfit()的结果看, 标准化残差及其残差的白噪声检验通过, 分布的拟合优度检验通过, 模型参数稳定性检验通过, 杠杆效应的检验通过。

19.1.5 预测

以EGARCH(1,1)为例说明EGARCH模型的超前多步预测。 设模型参数已知, \(\varepsilon_t\)服从标准正态分布。 模型为(见(19.5)(19.1)(19.2)\[\begin{aligned} \ln\sigma_t^2 =& \alpha_0 + \beta_1 \ln\sigma_{t-1}^2 + \alpha_1 \varepsilon_{t-1} + \gamma_1(|\varepsilon_{t-1}| - E|\varepsilon_{t-1}| ) . \end{aligned}\] 其中\(E|\varepsilon_{t-1}| = \sqrt{\frac{2}{\pi}}\)。 记 \[ g(x) = \alpha_1 x + \gamma_1(|x| - \sqrt{\frac{2}{\pi}}), \] 则模型可以写成 \[ \ln\sigma_t^2 = \alpha_0 + \beta_1 \ln\sigma_{t-1}^2 + g(\varepsilon_{t-1}), \] 将模型写成指数形式为 \[\begin{aligned} \sigma_t^2 = e^{\alpha_0} \sigma_{t-1}^{2\beta_1} e^{g(\varepsilon_{t-1})} . \end{aligned}\]

\(h\)为预测原点, 超前一步预测是对\(\sigma_{h+1}^2\)的预测,而 \[\begin{aligned} \sigma_{h+1}^2 = e^{\alpha_0} \sigma_{h}^{2\beta_1} e^{g(\varepsilon_{h})} \in F_h , \end{aligned}\] 所以超前一步预测为 \[\begin{align} \hat\sigma_h^2(1) = E(\sigma_{h+1}^2 | F_h) = \sigma_{h+1}^2 . \tag{19.6} \end{align}\]

\(t=h+2\)\[\begin{aligned} \sigma_{h+2}^2 = e^{\alpha_0} \sigma_{h+1}^{2\beta_1} e^{g(\varepsilon_{h+1})} , \end{aligned}\] 注意到\(\sigma_{h+1}^2 \in F_h\)\(\varepsilon_{h+1}\)\(F_h\)独立,所以有 \[\begin{aligned} \hat\sigma_h^2(2) =& E(\sigma_{h+2}^2 | F_h) = e^{\alpha_0} \sigma_{h+1}^{2\beta_1} E[e^{g(\varepsilon_{h+1})} | F_h] \\ =& e^{\alpha_0} \sigma_{h+1}^{2\beta_1} E[e^{g(\varepsilon_{h+1})}] . \end{aligned}\]

只需要计算\(E[e^{g(\varepsilon_{h+1})}]\)。 设\(\varepsilon \sim \text{N}(0,1)\),则 \[\begin{aligned} E e^{g(\varepsilon)} =& \int_{-\infty}^\infty \exp\left\{ \alpha_1 u + \gamma_1 |u| - \gamma_1 \sqrt{\frac{2}{\pi}} \right\} \varphi(u) \,du \\ =& e^{- \gamma_1 \sqrt{\frac{2}{\pi}}} \left\{ \int_0^\infty \frac{1}{\sqrt{2\pi}} \exp\left[ -\frac12 u^2 + (\alpha_1 + \gamma_1)u \right] \,du \right. \\ & \left. + \int_{-\infty}^0 \frac{1}{\sqrt{2\pi}} \exp\left[ -\frac12 u^2 + (\alpha_1 - \gamma_1)u \right] \,du \right\} \end{aligned}\] 其中 \[\begin{aligned} -\frac12 u^2 + (\alpha_1 \pm \gamma_1)u =& -\frac12 [ u - (\alpha_1 \pm \gamma_1)]^2 + \frac12 (\alpha_1 \pm \gamma_1)^2 , \end{aligned}\] 所以 \[\begin{aligned} E e^{g(\varepsilon)} =& e^{- \gamma_1 \sqrt{\frac{2}{\pi}}} \left\{ e^{\frac12 (\alpha_1 + \gamma_1)^2} \int_0^\infty \frac{1}{\sqrt{2\pi}} \exp\left[ -\frac12 (u - (\alpha_1 + \gamma_1))^2 \right] \,du \right. \\ & \left. + e^{\frac12 (\alpha_1 - \gamma_1)^2} \int_{-\infty}^0 \frac{1}{\sqrt{2\pi}} \exp\left[ -\frac12 (u - (\alpha_1 - \gamma_1))^2 \right] \,du \right\} \\ =& e^{- \gamma_1 \sqrt{\frac{2}{\pi}}} \left\{ e^{\frac12 (\alpha_1 + \gamma_1)^2} \int_{-(\alpha_1 + \gamma_1)}^\infty \varphi(v) \,dv + e^{\frac12 (\alpha_1 - \gamma_1)^2} \int_{-\infty}^{-(\alpha_1 - \gamma_1)} \varphi(v) \,dv \right\} \\ =& e^{- \gamma_1 \sqrt{\frac{2}{\pi}}} \left\{ e^{\frac12 (\alpha_1 + \gamma_1)^2} [ 1 - \Phi(-(\alpha_1 + \gamma_1))] + e^{\frac12 (\alpha_1 - \gamma_1)^2} \Phi(-(\alpha_1 - \gamma_1)) \right\} \\ =& e^{- \gamma_1 \sqrt{\frac{2}{\pi}}} \left\{ e^{\frac12 (\gamma_1 + \alpha_1)^2} \Phi(\gamma_1 + \alpha_1) + e^{\frac12 (\gamma_1 - \alpha_1)^2} \Phi(\gamma_1 - \alpha_1) \right\} . \end{aligned}\] 其中\(\varphi(\cdot)\)\(\Phi(\cdot)\)是标准正态分布的密度和分布函数。

于是,超前二步预测为 \[\begin{align} \hat\sigma_h^2(2) =& e^{\alpha_0} E e^{g(\varepsilon)} [\hat\sigma_h^2(1)]^{\beta_1} . \tag{19.7} \end{align}\]

超前三步: \[\begin{aligned} \sigma_{h+3}^2 = e^{\alpha_0} (\sigma_{h+2}^{2})^{\beta_1} e^{g(\varepsilon_{h+2})}, \end{aligned}\] 注意到\(F_h \subset F_{h+1}\), \(\sigma_{h+2}^2 \in F_{h+1}\), \(\varepsilon_{h+2}\)\(F_{h+1}\)独立, 可得 \[\begin{align} \hat\sigma_h^2(3) =& e^{\alpha_0} E \left\{ (\sigma_{h+2}^2)^{\beta_1} e^{g(\varepsilon_{h+2})} | F_h \right\} \\ =& e^{\alpha_0} E \left\{ E\left[ (\sigma_{h+2}^2)^{\beta_1} e^{g(\varepsilon_{h+2})} | F_{h+1} \right] | F_h \right\} \\ =& e^{\alpha_0} E \left\{ (\sigma_{h+2}^2)^{\beta_1} E\left[ e^{g(\varepsilon_{h+2})} | F_{h+1} \right] | F_h \right\} \\ =& e^{\alpha_0} E[e^{g(\varepsilon)}] E \left\{ (\sigma_{h+2}^2)^{\beta_1} | F_h \right\} \\ =& e^{\alpha_0} E[e^{g(\varepsilon)}] E \left\{ e^{\alpha_0 \beta_1} (\sigma_{h+1}^2)^{\beta_1^2} e^{\beta_1 g(\varepsilon_{h+1})} | F_h \right\} \\ =& e^{\alpha_0(1 + \beta_1)} E[e^{g(\varepsilon)}] E[e^{\beta_1 g(\varepsilon)}] (\hat\sigma_{h}^2(1))^{\beta_1^2} . \tag{19.8} \end{align}\]

\(\hat\sigma^2_h(j)\)可类似计算。

19.2 GJR-GARCH模型

GJR-GARCH模型是另一个能够反映杠杆效应的波动率模型, 参见(Glosten, Jagannathan, and Runkle 1993)(Zakoian 1994)。 或称为TGARCH。 GJR-GARCH(\(m,s\))模型的形式为 \[\begin{align} \sigma_t^2 = \alpha_0 + \sum_{i=1}^m (\alpha_i + \gamma_i N_{t-i}) a_{t-i}^2 + \sum_{j=1}^s \beta_j \sigma_{t-j}^2 \tag{19.9} \end{align}\] 其中\(N_{t-i}\)是表示\(a_{t-i}<0\)的示性函数,即 \[\begin{aligned} N_{t-i} = \begin{cases} 1, & a_{t-i} < 0 \\ 0, & a_{t-i} \geq 0 \end{cases} \end{aligned}\] \(\alpha_i, \gamma_i, \beta_j\)是非负参数, 满足与GARCH模型类似的参数条件。

正的\(a_{t-i}\)\(\sigma_t^2\)的影响为\(\alpha_i a_{t-i}^2\), 负的\(a_{t-i}\)\(\sigma_t^2\)的影响为\((\alpha_i + \gamma_i) a_{t-i}^2\), 当\(\gamma_i>0\)时负的\(a_{t-i}\)影响较大。

模型用0作为\(a_{t-i}\)的门限, 还可以用其它实数值作为门限。 关于门限模型参见(Tsay 2010)第4章。

作为例子, 考虑从1999年1月4日到2010年8月20日的欧元对美元汇率的日对数收益率序列。

读入数据:

d.useu <- read_table(
  "d-useu9910.txt", 
  col_types=cols(.default=col_double())
)
xts.useu <- with(d.useu, xts(rate, make_date(year, mon, day)))
xts.useu.lnrtn <- diff(log(xts.useu))[-1,]*100

用rugarch包建模:

library(rugarch)
spec2 <- ugarchspec(
  mean.model = list(
    armaOrder=c(0,0),
    include.mean=TRUE  ),
  variance.model = list(
    model = "gjrGARCH", # GJR-GARCH or TGARCH model
    garchOrder = c(1,1) ) )
mod2ru <- ugarchfit(spec = spec2, data = xts.useu.lnrtn)
show(mod2ru)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : gjrGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.012262    0.010726   1.1432 0.252951
## omega   0.001274    0.000554   2.3013 0.021373
## alpha1  0.022402    0.004277   5.2373 0.000000
## beta1   0.968711    0.001779 544.4570 0.000000
## gamma1  0.012435    0.007027   1.7696 0.076789
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.012262    0.011704    1.0476 0.294800
## omega   0.001274    0.000620    2.0535 0.040021
## alpha1  0.022402    0.004777    4.6891 0.000003
## beta1   0.968711    0.000795 1219.2085 0.000000
## gamma1  0.012435    0.008047    1.5454 0.122258
## 
## LogLikelihood : -2731.85 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       1.8688
## Bayes        1.8790
## Shibata      1.8688
## Hannan-Quinn 1.8725
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.2619  0.6088
## Lag[2*(p+q)+(p+q)-1][2]    0.3038  0.7936
## Lag[4*(p+q)+(p+q)-1][5]    3.3494  0.3469
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      4.935 0.02632
## Lag[2*(p+q)+(p+q)-1][5]     7.136 0.04806
## Lag[4*(p+q)+(p+q)-1][9]     9.101 0.07748
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     2.979 0.500 2.000 0.08435
## ARCH Lag[5]     3.239 1.440 1.667 0.25696
## ARCH Lag[7]     4.681 2.315 1.543 0.25882
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.4726
## Individual Statistics:              
## mu     0.20009
## omega  0.20231
## alpha1 0.06618
## beta1  0.08927
## gamma1 0.08316
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.28 1.47 1.88
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value     prob sig
## Sign Bias            1.663 0.096368   *
## Negative Sign Bias   2.386 0.017101  **
## Positive Sign Bias   1.740 0.082025   *
## Joint Effect        12.898 0.004863 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     59.51    4.628e-06
## 2    30     70.40    2.676e-05
## 3    40     69.44    1.928e-03
## 4    50     81.46    2.454e-03
## 
## 
## Elapsed time : 0.206059

估计的模型为 \[\begin{aligned} r_t =& 0.0123 + a_t, \quad a_t = \sigma_t \varepsilon_t, \quad \varepsilon_t \text{ i.i.d. N}(0,1) \\ \sigma_t^2 =& 0.0013 + (0.0224 + 0.0124 N_{t-1}) a_{t-1}^2 + 0.9687 \sigma_{t-1}^2, \end{aligned}\] 模型结果中mu不显著, 在0.05水平下gamma1也不显著, 说明杠杆效应不显著。 模型标准化残差的白噪声检验通过, 但标准化残差平方的白噪声检验没有通过, 所以模型效果不够好。 模型参数稳定性检验没有通过。 拟合优度检验没有通过。

拟合的波动率\(\sigma_t\)序列图形(单位是百分之一):

plot(sigma(mod2ru), 
  format.labels="%Y",
  main="欧元汇率日对数收益率波动率TGARCH估计(%)", 
  major.ticks="years", minor.ticks=NULL, 
  grid.ticks.on="years")
欧元汇率日对数收益率波动率TGARCH估计

图19.8: 欧元汇率日对数收益率波动率TGARCH估计

波动率的最高值出现在2008年末、2009年初, 次贷危机爆发时期。

为了进行模型诊断,计算标准化残差:

stdresi2 <- residuals(mod2ru, standardize=TRUE)

标准化残差的ACF:

forecast::Acf(stdresi2, lag.max=36, main="")
欧元汇率日对数收益率TGARCH标准化残差ACF

图19.9: 欧元汇率日对数收益率TGARCH标准化残差ACF

## 或:
## plot(mod2ru, which=10)

标准化残差平方的ACF:

forecast::Acf(stdresi2^2, lag.max=36, main="")
欧元汇率日对数收益率TGARCH标准化残差平方ACF

图19.10: 欧元汇率日对数收益率TGARCH标准化残差平方ACF

## 或:
## plot(mod2ru, which=11)

标准化残差的Ljung-Box白噪声检验:

Box.test(stdresi2, lag=12, type="Ljung")
## 
##  Box-Ljung test
## 
## data:  stdresi2
## X-squared = 15.029, df = 12, p-value = 0.2399

标准化残差平方的Ljung-Box白噪声检验:

Box.test(stdresi2^2, lag=12, type="Ljung")
## 
##  Box-Ljung test
## 
## data:  stdresi2^2
## X-squared = 15.094, df = 12, p-value = 0.2363

这些残差诊断的结果说明拟合的模型是基本符合的。

作标准化残差的盒形图:

boxplot(c(coredata(stdresi2)), 
  main="", xlab="标准化残差")
欧元汇率日对数收益率TGARCH标准化残差盒形图

图19.11: 欧元汇率日对数收益率TGARCH标准化残差盒形图

标准化残差的正态QQ图:

plot(mod2ru, which=9)
欧元汇率日对数收益率TGARCH标准化残差正态QQ图

图19.12: 欧元汇率日对数收益率TGARCH标准化残差正态QQ图

从标准化残差的盒形图看有厚尾现象。

可以将GJR-GARCH得到的波动率与标准GARCH得到的波动率进行比较:

library(fGarch, quietly = TRUE)
mod2sg <- garchFit( ~ 1 + garch(1,1), 
  data=xts.useu.lnrtn, trace=FALSE)
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
##   Consider formula(paste(x, collapse = " ")) instead.

将两个波动率作图, 黑色为GARCH结果, 红色虚线为GJR-GARCH结果:

plot(cbind(volatility(mod2sg), sigma(mod2ru)), 
  main="欧元汇率日对数收益率波动率GARCH和NGARCH估计", 
  lty=c(1,2), 
  col=c("black", "red"),
  format.labels="%Y",
  major.ticks="years", minor.ticks=NULL, 
  grid.ticks.on="years")
欧元汇率GARCH和NGARCH估计的波动率

图19.13: 欧元汇率GARCH和NGARCH估计的波动率

放大显示其中2008年以后:

plot(cbind(volatility(mod2sg), sigma(mod2ru)), 
  subset="2008/",
  main="欧元汇率日对数收益率波动率GARCH和NGARCH估计", 
  lty=c(1,2), 
  col=c("black", "red"),
  format.labels="%Y",
  major.ticks="years", minor.ticks=NULL, 
  grid.ticks.on="years")
欧元汇率GARCH和NGARCH估计的波动率

图19.14: 欧元汇率GARCH和NGARCH估计的波动率

19.3 APARCH模型

(Ding, Granger, and Engle 1993)提出了非对称幂(asymmetric power) ARCH模型(APARCH模型), 模型形式为 \[\begin{align} r_t =& \mu_t + a_t, \quad a_t = \sigma_t \varepsilon_t, \quad \varepsilon_t \sim D(0,1) \\ \sigma_t^\delta =& \omega + \sum_{i=1}^m \alpha_i (|a_{t-i}| - \gamma_i a_{t-i})^\delta + \sum_{j=1}^s \beta_j \sigma_{t-j}^\delta \tag{19.10} \end{align}\] 其中\(\mu_t\)是条件均值, \(D(0,1)\)表示某个零均值单位方差分布, \(\delta\)为正实数, 系数\(\omega, \alpha_i, \gamma_i, \beta_j\)满足某些正则性条件使得波动率为正。

最常用的是最简单的APARCH(1,1)模型。 这个模型中包含了许多其它模型。

\(\delta=2\)\(\gamma_j = 0\)时即普通的GARCH模型。

\(\delta=2\)时即TGARCH模型(形式略有不同)。

\(\delta=1\)时波动率方程直接使用波动率\(\sigma_{t}\)和新息\(a_t\)而非其平方。

APARCH中的幂变换旨在提高拟合程度, 但幂次\(\delta\)没有很好解释。

R的fGarch::garchFit()函数中可以使用aparch(m,s)作为模型设定。 作为例子, 拟合欧元对美元汇率数据:

library(fGarch, quietly = TRUE)
modres3 <- garchFit( ~ 1 + aparch(1,1), 
  data=xts.useu.lnrtn, trace=FALSE)
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
##   Consider formula(paste(x, collapse = " ")) instead.
summary(modres3)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~1 + aparch(1, 1), data = xts.useu.lnrtn, 
##     trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ 1 + aparch(1, 1)
## <environment: 0x0000000027bf8bb0>
##  [data = xts.useu.lnrtn]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##        mu      omega     alpha1     gamma1      beta1      delta  
## 0.0127647  0.0015919  0.0313680  0.1135338  0.9689155  1.6743169  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu     0.0127647   0.0107626    1.186   0.2356    
## omega  0.0015919   0.0007226    2.203   0.0276 *  
## alpha1 0.0313680   0.0053350    5.880 4.11e-09 ***
## gamma1 0.1135338   0.0711908    1.595   0.1108    
## beta1  0.9689155   0.0038405  252.291  < 2e-16 ***
## delta  1.6743169   0.4057086    4.127 3.68e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -2731.172    normalized:  -0.9324587 
## 
## Description:
##  Wed May 11 22:04:40 2022 by user: Lenovo 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  50.20527  1.253331e-11
##  Shapiro-Wilk Test  R    W      0.9956711 1.608392e-07
##  Ljung-Box Test     R    Q(10)  13.37689  0.2033563   
##  Ljung-Box Test     R    Q(15)  20.19634  0.1645294   
##  Ljung-Box Test     R    Q(20)  22.84736  0.2963513   
##  Ljung-Box Test     R^2  Q(10)  13.15611  0.2150739   
##  Ljung-Box Test     R^2  Q(15)  16.58009  0.3445799   
##  Ljung-Box Test     R^2  Q(20)  27.44885  0.1231015   
##  LM Arch Test       R    TR^2   14.35739  0.2784705   
## 
## Information Criterion Statistics:
##      AIC      BIC      SIC     HQIC 
## 1.869014 1.881269 1.869006 1.873428

拟合的白噪声检验说明模型是充分的, 但是\(\delta\)的估计为1.67, 标准误差为0.41, 其95%置信区间为\((0.85, 2.49)\)包含2, 所以与\(\delta=2\)没有显著差异, 可以用TGARCH模型。 \(\gamma_1\)在0.05水平下不显著, 可以不需要反映杠杆效应的模型。

拟合的APARCH模型为: \[\begin{aligned} r_t =& 0.0128 + a_t, \quad a_t = \sigma_t \varepsilon_t, \quad \varepsilon_t \sim \text{N}(0,1) \\ \sigma_t^{1.67} =& 0.0016 + 0.0313 (|a_{t-1}| - 0.1135 a_{t-1})^{1.67} + 0.9689 \sigma_{t-1}^{1.67} \end{aligned}\]

固定\(\delta=2\)估计APARCH(1,1)模型:

modres4 <- garchFit( ~ 1 + aparch(1,1), 
  data=xts.useu.lnrtn, delta=2, 
  include.delta=FALSE, trace=FALSE)
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
##   Consider formula(paste(x, collapse = " ")) instead.
summary(modres4)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~1 + aparch(1, 1), data = xts.useu.lnrtn, 
##     delta = 2, include.delta = FALSE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ 1 + aparch(1, 1)
## <environment: 0x00000000287623c0>
##  [data = xts.useu.lnrtn]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##        mu      omega     alpha1     gamma1      beta1  
## 0.0122646  0.0012745  0.0282723  0.1100241  0.9687115  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu     0.0122646   0.0107289    1.143   0.2530    
## omega  0.0012745   0.0005752    2.216   0.0267 *  
## alpha1 0.0282723   0.0038637    7.317 2.53e-13 ***
## gamma1 0.1100241   0.0649051    1.695   0.0900 .  
## beta1  0.9687115   0.0039421  245.736  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -2731.85    normalized:  -0.9326902 
## 
## Description:
##  Wed May 11 22:04:40 2022 by user: Lenovo 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  49.97675  1.405032e-11
##  Shapiro-Wilk Test  R    W      0.9956803 1.655882e-07
##  Ljung-Box Test     R    Q(10)  13.38285  0.203047    
##  Ljung-Box Test     R    Q(15)  20.29833  0.1607845   
##  Ljung-Box Test     R    Q(20)  22.87265  0.2950909   
##  Ljung-Box Test     R^2  Q(10)  12.89585  0.2295534   
##  Ljung-Box Test     R^2  Q(15)  16.55288  0.3462878   
##  Ljung-Box Test     R^2  Q(20)  27.24036  0.1286361   
##  LM Arch Test       R    TR^2   14.29661  0.2821698   
## 
## Information Criterion Statistics:
##      AIC      BIC      SIC     HQIC 
## 1.868795 1.879007 1.868789 1.872472

拟合的固定\(\delta=2\)的APARCH模型为: \[\begin{aligned} r_t =& 0.0123 + a_t, \quad a_t = \sigma_t \varepsilon_t, \quad \varepsilon_t \sim \text{N}(0,1) \\ \sigma_t^2 =& 0.00127 + 0.0283 (|a_{t-1}| - 0.1100 a_{t-1})^2 + 0.9687 \sigma_{t-1}^2 \end{aligned}\] 其中的波动率方程可以借助\(N_{t-1}\)\(a_{t-1}<0\)的示性函数写成 \[ \sigma_t^2 = 0.00127 + (0.0224 + 0.0125 N_{t-1}) a_{t-1}^2 + 0.9687 \sigma_{t-1}^2 . \]

与直接估计TGARCH(1,1)的结果对比: \[\begin{aligned} r_t =& 0.0123 + a_t, \quad a_t = \sigma_t \varepsilon_t, \quad \varepsilon_t \text{ i.i.d. N}(0,1) \\ \sigma_t^2 =& 0.0013 + (0.0224 + 0.0124 N_{t-1}) a_{t-1}^2 + 0.9687 \sigma_{t-1}^2 . \end{aligned}\] 两个估计结果基本相同。

也可以用rugarch包估计APARCH模型:

library(rugarch)
spec3 <- ugarchspec(
  mean.model = list(
    armaOrder=c(0,0),
    include.mean=TRUE  ),
  variance.model = list(
    model = "apARCH", # APARCH model
    garchOrder = c(1,1) ) )
mod3ru <- ugarchfit(spec = spec3, 
  data = xts.useu.lnrtn)
show(mod3ru)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : apARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.012516    0.010740    1.1654 0.243858
## omega   0.001683    0.000595    2.8299 0.004656
## alpha1  0.031692    0.000877   36.1495 0.000000
## beta1   0.969388    0.000679 1427.5230 0.000000
## gamma1  0.121817    0.072545    1.6792 0.093115
## delta   1.580101    0.129609   12.1913 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.012516    0.011800    1.0607 0.288841
## omega   0.001683    0.000772    2.1804 0.029229
## alpha1  0.031692    0.002833   11.1859 0.000000
## beta1   0.969388    0.000680 1425.6414 0.000000
## gamma1  0.121817    0.084847    1.4357 0.151084
## delta   1.580101    0.082927   19.0542 0.000000
## 
## LogLikelihood : -2731.113 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       1.8690
## Bayes        1.8812
## Shibata      1.8690
## Hannan-Quinn 1.8734
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.2823  0.5952
## Lag[2*(p+q)+(p+q)-1][2]    0.3207  0.7843
## Lag[4*(p+q)+(p+q)-1][5]    3.3970  0.3393
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                      5.165 0.02304
## Lag[2*(p+q)+(p+q)-1][5]     7.349 0.04264
## Lag[4*(p+q)+(p+q)-1][9]     9.341 0.06909
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     2.896 0.500 2.000 0.08882
## ARCH Lag[5]     3.238 1.440 1.667 0.25703
## ARCH Lag[7]     4.723 2.315 1.543 0.25404
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.5207
## Individual Statistics:              
## mu     0.18789
## omega  0.18056
## alpha1 0.07682
## beta1  0.09516
## gamma1 0.38035
## delta  0.08651
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.49 1.68 2.12
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value     prob sig
## Sign Bias            1.641 0.100814    
## Negative Sign Bias   2.398 0.016537  **
## Positive Sign Bias   1.771 0.076666   *
## Joint Effect        12.991 0.004657 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     62.92    1.323e-06
## 2    30     71.22    2.064e-05
## 3    40     73.65    6.628e-04
## 4    50     79.79    3.561e-03
## 
## 
## Elapsed time : 0.4970009

与fGarch包的估计结果不完全相同, 参数值比较接近。

B 参考文献

Ding, Z., C. W. J. Granger, and R. F. Engle. 1993. “A Long Memory Property of Stock Returns and a New Model.” Journal of Empirical Finance 1: 83–106.
Glosten, L. R., R. Jagannathan, and D. E. Runkle. 1993. “On the Relation Between the Expected Value and the Volatility of Nominal Excess Return on Stocks.” J. Finance 48: 1779–1801.
Nelson, D. B. 1991. “Conditional Heteroskedasticity in Asset Returns: A New Approach.” Econometrica 59: 347–70.
Tsay, Ruey S. 2010. Analysis of Financial Time Series. 3rd Ed. John Wiley & Sons, Inc.
———. 2013. 金融数据分析导论:基于R语言. 机械工业出版社.
Zakoian, J. M. 1994. “Threshold Heteroscedastic Models.” J Econ Dyn Control 18: 931–55.