18 GARCH模型

本章来自(Tsay 2013)§4.6-4.8内容。

18.1 GARCH模型

ARCH模型用来描述波动率能得到很好的效果, 但实际建模时可能需要较高的阶数, 比如§17.5.3的欧元汇率波动率建模用了11阶的ARCH模型。

18.1.1 模型方程

(Bollerslev 1986)提出了ARCH模型的一种重要推广模型, 称为GARCH模型。 对于一个对数收益率序列\(r_t\), 令\(a_t = r_t - \mu_t = r_t - E(r_t | F_{t-1})\)为其新息序列, 称\(\{ a_t \}\)服从GARCH(\(m,s\))模型, 如果\(a_t\)满足 \[\begin{align} a_t = \sigma_t \varepsilon_t, \quad \sigma_t^2 = \alpha_0 + \sum_{i=1}^m \alpha_i a_{t-i}^2 + \sum_{j=1}^s \beta_j \sigma_{t-j}^2 \tag{18.1} \end{align}\] 其中\(\{ \varepsilon_t \}\)为零均值单位方差的独立同分布白噪声列, \(\alpha_0>0\), \(\alpha_i \geq 0\), \(\beta_j \geq 0\), \(0 < \sum_{i=1}^m \alpha_i + \sum_{j=1}^s \beta_j < 1\), 这最后一个条件用来保证满足模型的\(a_t\)的无条件方差有限且不变, 而条件方差\(\sigma_t^2\)可以随时间\(t\)而变。

18.1.2 与ARMA模型比较

模型(18.1)像是ARMA(\(p,q\))模型, 但是这里的\(\sigma_{t-i}^2\)\(a_{t-i}^2\)是有关系的, \(\sigma_{t-i}^2\)\(a_{t-i}\)的条件方差, ARMA模型中的\(x_{t-i}\)\(\varepsilon_{t-i}\)并没有这样的关系。

为了利用GARCH模型与ARMA模型的相似性, 令\(\alpha_i=0\), 当\(i>m\); 令\(\beta_j = 0\), 当\(j>s\)。 令\(\eta_t = a_t^2 - \sigma_t^2\), 下一小节证明了\(E a_t = 0\), 而\(\sigma_t^2 = \text{Var}(a_t | F_{t-1}) = E(a_t^2 | F_{t-1})\), 当\(a_t\)为严平稳列时\(\eta_t\)是鞅差序列,这是比宽白噪声严一些, 比零均值独立同分布白噪声宽一些的条件。 将\(\sigma_{t-i}^2 = a_{t-i}^2 - \eta_{t-i}\)代入模型(18.1)\[\begin{align} a_t^2 = \alpha_0 + \sum_{i=1}^{\max(m,s)} (\alpha_i + \beta_i) a_{t-i}^2 + \eta_t - \sum_{j=1}^s \beta_j \eta_{t-j} \tag{18.2} \end{align}\] 这就是关于\(\{a_t^2\}\)的ARMA(\(\max(m,s)\), \(s\))模型, 由ARMA模型的无条件期望的公式得 \[ E a_t^2 = \frac{\alpha_0}{1 - \sum_{i=1}^{\max(m,s)} (\alpha_i + \beta_i)} = \frac{\alpha_0}{1 - \sum_{i=1}^m \alpha_i - \sum_{j=1}^s \beta_j} \] 这要求分母为正,即要求\(\sum_{i=1}^m \alpha_i + \sum_{j=1}^s \beta_j < 1\)。 这时\(a_t\)的无条件方差\(\text{Var}(a_t)\)也等于上式。

18.1.3 GARCH模型的性质

下面以最简单的GARCH(1,1)为例研究GARCH模型的性质。 令\(F_{t-1}\)表示截止到\(t-1\)时刻的\(a_{t-i}\)\(\sigma_{t-j}\)所包含的信息。 模型为 \[\begin{align} a_t =& \sigma_t \varepsilon_t, \quad \varepsilon_t \text{ i.i.d. WN} (0,1)\\ \sigma_t^2 =& \alpha_0 + \alpha_1 a_{t-1}^2 + \beta_1 \sigma_{t-1}^2 \tag{18.3} \end{align}\]

为了计算无条件均值\(Ea_t\),先计算条件期望 \[ E(a_t | F_{t-1}) = E(\sigma_t \varepsilon_t | F_{t-1}) = \sigma_t E(\varepsilon_t | F_{t-1}) = 0 \] 这里用了\(\sigma_t \in F_{t-1}\)\(\varepsilon_t\)\(F_{t-1}\)独立。 于是 \[ E a_t = E[ E(a_t | F_{t-1})] = 0 \] 即GARCH模型的新息\(a_t\)的无条件期望为零。

来计算\(a_t\)的无条件方差。 设模型(18.1)\(\{ a_t \}\)序列存在严平稳解,则 \[ \begin{aligned} \text{Var}(a_t) =& E(a_t^2) = E[ E(a_t^2 | F_{t-1})] = E[ E(\sigma_t^2 \varepsilon_t^2 | F_{t-1})]\\ =& E[\sigma_t^2 E(\varepsilon_t^2 | F_{t-1})] = E[\sigma_t^2 E(\varepsilon_t^2)] \\ =& E[\sigma_t^2] = E[\alpha_0 + \alpha_1 a_{t-1}^2 + \beta_1 \sigma_{t-1}^2] \\ =& \alpha_0 + \alpha_1 E(a_{t-1}^2) + \beta_1 E[E(a_{t-1}^2|F_{t-2})]\\ =& \alpha_0 + (\alpha_1 + \beta_1) E(a_{t-1}^2) \end{aligned} \]\(E a_t^2 = E a_{t-1}^2\), 解得 \[ \text{Var}(a_t) = E a_t^2 = \frac{\alpha_0}{1 - \alpha_1 - \beta_1} . \]

GARCH(1,1)模型的性质:

第一,像ARCH模型一样,\(a_t\)存在波动率聚集, 一个较大的\(a_{t-1}\)\(\sigma_{t-1}\)使得\(1\)步以后的条件方差变大, 从而倾向于出现较大的对数收益率。

第二,当\(\varepsilon_t\)为标准正态分布时, 在如下条件下\(a_t\)有无条件四阶矩: \[1 - 2 \alpha_1^2 - (\alpha_1 + \beta_1)^2 > 0\] 这时超额峰度为 \[ \frac{E a_t^4}{(Ea_t^2)^2} - 3 = \frac{2\left[1 - (\alpha_1 + \beta_1)^2 + \alpha_1^2 \right]} {1 - (\alpha_1 + \beta_1)^2 - 2\alpha_1^2} > 0 \]\(a_t\)分布厚尾。

第三,GARCH模型给出了一个比较简单的波动率模型。

18.1.4 预测

可以用类似ARMA预测的方法预测波动率。 仍以GARCH(1,1)为例, 由模型(18.3), 基于截止到\(h\)时刻的观测作超前一步预测: \[ \sigma_{h+1}^2 = \alpha_0 + \alpha_1 a_{h}^2 + \beta_1 \sigma_{h}^2 \in F_{h} \] 所以 \[\begin{align} \sigma_h^2(1) = E(\sigma_{h+1}^2 | F_{h}) = \sigma_{h+1}^2 = \alpha_0 + \alpha_1 a_{h}^2 + \beta_1 \sigma_{h}^2 . \tag{18.4} \end{align}\]\(\sigma_{h+2}^2\),利用\(a_t^2 = \sigma_t^2 \varepsilon_t^2\),有 \[ \begin{aligned} \sigma_{h+2}^2 =& \alpha_0 + \alpha_1 a_{h+1}^2 + \beta_1 \sigma_{h+1}^2 \\ =& \alpha_0 + \alpha_1 \sigma_{h+1}^2 \varepsilon_{h+1}^2 + \beta_1 \sigma_{h+1}^2 \\ =& \alpha_0 + (\alpha_1 \varepsilon_{h+1}^2 + \beta_1) \sigma_{h+1}^2 \end{aligned} \] 于是 \[ \sigma_h^2(2) = E(\sigma_{h+2}^2 | F_{h}) = \alpha_0 + E(\alpha_1 \varepsilon_{h+1}^2 + \beta_1 | F_h) \sigma_{h+1}^2 = \alpha_0 + (\alpha_1 + \beta_1) \sigma_h^2(1) \] 类似地,对\(\ell \geq 2\)\[ \sigma_{h+\ell}^2 = \alpha_0 + \alpha_1 \varepsilon_{h+\ell-1}^2 \sigma_{h+\ell-1}^2 + \beta_1 \sigma_{h+\ell-1}^2 = \alpha_0 + (\alpha_1 \varepsilon_{h+\ell-1}^2 + \beta_1) \sigma_{h+\ell-1}^2 \] 于是 \[\begin{align} \sigma_h^2(\ell) =& E\left\{ \sigma_{h+\ell}^2 | F_h \right\} = \alpha_0 + E\left\{ (\alpha_1 \varepsilon_{h+\ell-1}^2 + \beta_1) \sigma_{h+\ell-1}^2 | F_h \right\} \\ =& \alpha_0 + E\left\{ E\left[ (\alpha_1 \varepsilon_{h+\ell-1}^2 + \beta_1) \sigma_{h+\ell-1}^2 | F_{h+\ell-2} \right] | F_h \right\} \\ =& \alpha_0 + E\left\{ \sigma_{h+\ell-1}^2 E\left[ \alpha_1 \varepsilon_{h+\ell-1}^2 + \beta_1 | F_{h+\ell-2} \right] | F_h \right\} \quad(\text{注意}\sigma_{h+\ell-1}^2 \in F_{h+\ell-2}) \\ =& \alpha_0 + \left\{ \sigma_{h+\ell-1}^2 (\alpha_1 + \beta_1) | F_h \right\} \\ =& \alpha_0 + (\alpha_1 + \beta_1) \sigma_h^2(\ell-1) \tag{18.5} \end{align}\]

预测公式与自回归系数为\((\alpha_1 + \beta_1)\)的ARMA(1,1)的超前预测公式相同。

\(\ell=2\)迭代计算得 \[ \sigma_h^2(\ell) = \frac{\alpha_0[1 - (\alpha_1 + \beta_1)^{\ell-1}]}{1 - (\alpha_1 + \beta_1)} + (\alpha_1 + \beta_1)^{(\ell-1)} \sigma_h^2(1) \] 只要\(\alpha_1 + \beta_1 < 1\)就有 \[ \sigma_h^2(\ell) \to \frac{\alpha_0}{1 - \alpha_1 - \beta_1} = \text{Var}(a_t) \] 即超前多步条件方差预测趋于\(a_t\)的无条件方差。

GARCH模型有和ARCH模型类似的弱点。 在高频数据研究发现即使使用t分布, 分布厚尾性也不足; 对于收益率的正负不对称性无法反映。

18.1.5 模型估计

ARCH模型的建模步骤也适用于GARCH模型的建模。 GARCH模型的定阶方法研究不多, 一般用试错法尝试较低阶的GARCH模型, 如GARCH(1,1), GARCH(2,1), GARCH(1,2)等。 许多情况下GARCH(1,1)就能解决问题。

为了估计参数, 可以假定初始的\(\sigma_t^2\)已知, 递推计算后续的\(\sigma_t^2\)并计算条件似然函数, 求条件似然函数的最大值点得到参数估计。 有时用\(a_t\)的样本方差作为初始的\(\sigma_t\)的值。

为了检验模型的充分性, 可以计算标准化残差 \[ \tilde a_t = \frac{a_t}{\sigma_t} \] 通过对\(\tilde a_t\)\(\tilde a_t^2\)的白噪声检验确认模型可以接受。

18.1.6 Intel公司股票收益率的波动率建模实例

继续使用§17.5.1中的Intel公司股票从1973-1到2009-12的月度对数收益率数据, 有444个观测值。 §17.5.1中的ARCH(1)模型(17.7)在模型检验中有一些不足, 比如关于标准化残差平方的Ljung-Box检验的滞后15和滞后20检验是显著的。 尝试用GARCH(1,1)模型来改进。记\(r_t\)为对数收益率序列。

数据读入:

d.intel <- read_table2(
  "m-intcsp7309.txt",
  col_types=cols(
    .default=col_double(),
    date=col_date("%Y%m%d")
  ))
xts.intel <- xts(
  log(1 + d.intel[["intc"]]), d.intel[["date"]]
)
tclass(xts.intel) <- "yearmon"
ts.intel <- ts(c(coredata(xts.intel)), start=c(1973,1), frequency=12)
at <- ts.intel - mean(ts.intel)

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

plot(ts.intel, ylab="log return", main="Intel Stock Price Monthly Log Return")
Intel月对数收益率时间序列

图18.1: Intel月对数收益率时间序列

18.1.6.1 使用条件正态分布

采用正态条件分布建立GARCH(1,1)模型:

library(fGarch, quietly = TRUE)
mod1 <- garchFit(~ 1 + garch(1,1), data=ts.intel, trace=FALSE)
summary(mod1)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~1 + garch(1, 1), data = ts.intel, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ 1 + garch(1, 1)
## <environment: 0x0000000017fd0a68>
##  [data = ts.intel]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##         mu       omega      alpha1       beta1  
## 0.01126568  0.00091902  0.08643831  0.85258554  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu     0.0112657   0.0053931    2.089  0.03672 *  
## omega  0.0009190   0.0003888    2.364  0.01808 *  
## alpha1 0.0864383   0.0265439    3.256  0.00113 ** 
## beta1  0.8525855   0.0394322   21.622  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  312.3307    normalized:  0.7034475 
## 
## Description:
##  Tue Apr 21 11:10:29 2020 by user: user 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  174.904   0           
##  Shapiro-Wilk Test  R    W      0.9709615 1.030282e-07
##  Ljung-Box Test     R    Q(10)  8.016844  0.6271916   
##  Ljung-Box Test     R    Q(15)  15.5006   0.4159946   
##  Ljung-Box Test     R    Q(20)  16.41549  0.6905368   
##  Ljung-Box Test     R^2  Q(10)  0.8746345 0.9999072   
##  Ljung-Box Test     R^2  Q(15)  11.35935  0.7267295   
##  Ljung-Box Test     R^2  Q(20)  12.55994  0.8954573   
##  LM Arch Test       R    TR^2   10.51401  0.5709617   
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -1.388877 -1.351978 -1.389037 -1.374326

对标准化残差及其平方的白噪声检验结果都通过了。 比§17.5.1要好, 而且ACI值\(-1.3889\)(17.7)的AIC值\(-1.3375\)也更好。 条件分布的正态性检验仍通不过。

模型可以写成: \[\begin{align} r_t =& 0.0113 + a_t, \quad a_t = \sigma_t \varepsilon_t, \quad \epsilon_t \text{ i.i.d. } \sim \text{N}(0,1) \\ \sigma_t^2 =& 0.00092 + 0.0864 a_{t-1}^2 + 0.8526 \sigma_{t-1}^2 \tag{18.6} \end{align}\] \(\sigma_t^2\)对过去的依赖主要来源于\(\beta_1 = 0.85\)

拟合的波动率图形:

##plot(mod1, which=2)
vola <- volatility(mod1)
plot(ts(vola, start=start(ts.intel), frequency=frequency(ts.intel)), 
     xlab="年", ylab="波动率")
abline(h=sd(ts.intel), col="green")
Intel股票建模拟合波动率

图18.2: Intel股票建模拟合波动率

可以看出波动率在1973-1974石油危机期间和2000年的互联网泡沫期间最高。 波动率图形中绿色横线为样本标准差,值为:

sd(ts.intel)
## [1] 0.1269101

按模型计算的对数收益率无条件标准差为:

tmpx <- coef(mod1)
unname(sqrt(tmpx["omega"]/(1 - tmpx["alpha1"] - tmpx["beta1"])))
## [1] 0.122767

模型计算的无条件标准差与样本标准差很接近。

标准化残差的时间序列图:

##plot(mod1, which=9)
resi <- residuals(mod1, standardize=TRUE)
plot(ts(resi, start=start(ts.intel), frequency=frequency(ts.intel)), 
     xlab="年", ylab="标准化残差")
Intel股票建模标准化残差

图18.3: Intel股票建模标准化残差

除了个别异常值,标准化残差看不出不符合独立同分布序列的特征。

\(\tilde a_t\)的ACF:

##plot(mod1, which=10)
resi <- residuals(mod1, standardize=TRUE)
acf(resi, lag.max=36, main="")
Intel股票建模标准化残差的ACF

图18.4: Intel股票建模标准化残差的ACF

\(\tilde a_t^2\)的ACF:

##plot(mod1, which=11)
resi <- residuals(mod1, standardize=TRUE)
acf(resi^2, lag.max=36, main="")
Intel股票建模标准化残差平方的ACF

图18.5: Intel股票建模标准化残差平方的ACF

\(\tilde a_t^2\)的ACF在滞后12处略超出。 总的来说模型是满意的。

\(\hat\mu + 2\hat\sigma_t\)可以作为\(r_t\)的近似95%置信区间, 将置信下限和置信下限分别延\(t\)轴方向连成曲线,得到如下图形:

##plot(mod1, which=11)
vola <- volatility(mod1)
hatmu <- coef(mod1)["mu"]
lb.intel <- hatmu - 2*vola
ub.intel <- hatmu + 2*vola
ylim <- range(c(ts.intel, lb.intel, ub.intel))
x.intel <- c(time(ts.intel))
plot(x.intel, c(ts.intel), type="l",
     xlab="年", ylab="对数收益率")
lines(x.intel, c(lb.intel), col="red")
lines(x.intel, c(ub.intel), col="red")
Intel股票对数收益率的逐点95%预测界限

图18.6: Intel股票对数收益率的逐点95%预测界限

可见对数收益率的取值基本都在预测区间之内。 预测是滚动的超前一步预测。

tseries包的garch()函数也可以用来拟合GARCH模型, 只能使用条件正态分布。 函数结果支持一系列的信息提取函数, 如summary()

18.1.6.2 使用条件t分布

应用条件t分布的新息,拟合模型:

library(fGarch, quietly = TRUE)
mod2 <- garchFit(~ 1 + garch(1,1), data=ts.intel, 
                 cond.dist="std", trace=FALSE)
summary(mod2)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~1 + garch(1, 1), data = ts.intel, cond.dist = "std", 
##     trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ 1 + garch(1, 1)
## <environment: 0x00000000191fad70>
##  [data = ts.intel]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##        mu      omega     alpha1      beta1      shape  
## 0.0165075  0.0011576  0.1059030  0.8171313  6.7723503  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu     0.0165075   0.0051031    3.235 0.001217 ** 
## omega  0.0011576   0.0005782    2.002 0.045286 *  
## alpha1 0.1059030   0.0372047    2.846 0.004420 ** 
## beta1  0.8171313   0.0580141   14.085  < 2e-16 ***
## shape  6.7723503   1.8572388    3.646 0.000266 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  326.2264    normalized:  0.734744 
## 
## Description:
##  Tue Apr 21 11:10:34 2020 by user: user 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  203.4933  0           
##  Shapiro-Wilk Test  R    W      0.9687607 3.970603e-08
##  Ljung-Box Test     R    Q(10)  7.877778  0.6407741   
##  Ljung-Box Test     R    Q(15)  15.5522   0.4124197   
##  Ljung-Box Test     R    Q(20)  16.50475  0.6848581   
##  Ljung-Box Test     R^2  Q(10)  1.066054  0.9997694   
##  Ljung-Box Test     R^2  Q(15)  11.49875  0.7165045   
##  Ljung-Box Test     R^2  Q(20)  12.61496  0.8932865   
##  LM Arch Test       R    TR^2   10.80739  0.5454935   
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -1.446966 -1.400841 -1.447215 -1.428776

模型的检验也比较满意。 模型可以写成: \[\begin{align} r_t =& 0.0165 + a_t, \quad a_t = \sigma_t \varepsilon_t, \quad \epsilon_t \text{ i.i.d. } \sim t^*(6.77) \\ \sigma_t^2 =& 0.00116 + 0.1059 a_{t-1}^2 + 0.8171 \sigma_{t-1}^2 \tag{18.7} \end{align}\] 其中\(t^*\)表示标准化t分布。 AIC为\(-1.447\), 正态分布时的AIC为\(-1.389\), t分布结果的样本内比较占优。

模型(18.7)隐含的\(a_t\)的无条件标准差为

tmpx <- coef(mod2)
unname(sqrt(tmpx["omega"]/(1 - tmpx["alpha1"] - tmpx["beta1"])))
## [1] 0.1226399

\(0.1226\),样本标准差为\(0.1269\), 正态时模型隐含的无条件标准差为\(0.1228\), 不同分布隐含的无条件标准差基本相同。

18.1.6.3 使用条件有偏t分布

对数收益率序列\(r_t\)的样本偏度为

skewness.test <- function(x, na.rm=TRUE){
  if(na.rm) x <- x[!is.na(x)]
  z <- (x - mean(x))/sd(x)
  n <- length(x)
  sk <- n/(n-1)/(n-2) * sum(z^3)
  t.sk <- sk / sqrt(6/n)
  p.sk <- 2*(1 - pnorm(abs(t.sk)))
  c(estimate=sk, statistic=t.sk, pvalue=p.sk)
}
skewness.test(c(ts.intel))
##      estimate     statistic        pvalue 
## -5.563719e-01 -4.786092e+00  1.700598e-06

偏度估计为\(-0.56\), 偏度为零的零假设的检验p值很小, 所以收益率分布是显著左偏(负偏)的。 这样, \(\varepsilon_t\)的分布最好采用左偏的分布。

在fGarch的garchFit()中指定cond.dist="sstd", 则条件分布为有偏的标准化t分布:

library(fGarch, quietly = TRUE)
mod3 <- garchFit(~ 1 + garch(1,1), data=ts.intel, 
                 cond.dist="sstd", trace=FALSE)
summary(mod3)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~1 + garch(1, 1), data = ts.intel, cond.dist = "sstd", 
##     trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ 1 + garch(1, 1)
## <environment: 0x00000000197b8720>
##  [data = ts.intel]
## 
## Conditional Distribution:
##  sstd 
## 
## Coefficient(s):
##        mu      omega     alpha1      beta1       skew      shape  
## 0.0133343  0.0011621  0.1049289  0.8177875  0.8717220  7.2344225  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu     0.0133343   0.0053430    2.496 0.012572 *  
## omega  0.0011621   0.0005587    2.080 0.037519 *  
## alpha1 0.1049289   0.0358860    2.924 0.003456 ** 
## beta1  0.8177875   0.0559863   14.607  < 2e-16 ***
## skew   0.8717220   0.0629129   13.856  < 2e-16 ***
## shape  7.2344225   2.1018042    3.442 0.000577 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  328.0995    normalized:  0.7389628 
## 
## Description:
##  Tue Apr 21 11:10:34 2020 by user: user 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  195.2178  0           
##  Shapiro-Wilk Test  R    W      0.9692506 4.892686e-08
##  Ljung-Box Test     R    Q(10)  7.882126  0.6403496   
##  Ljung-Box Test     R    Q(15)  15.62496  0.4074054   
##  Ljung-Box Test     R    Q(20)  16.5774   0.6802193   
##  Ljung-Box Test     R^2  Q(10)  1.078429  0.9997569   
##  Ljung-Box Test     R^2  Q(15)  11.95155  0.6826923   
##  Ljung-Box Test     R^2  Q(20)  13.03792  0.8757513   
##  LM Arch Test       R    TR^2   11.18826  0.5128574   
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -1.450899 -1.395550 -1.451257 -1.429071

模型为: \[\begin{align} r_t =& 0.0133 + a_t, \quad a_t = \sigma_t \varepsilon_t, \quad \epsilon_t \text{ i.i.d. } \sim t^*_{0.87, 7.23} \\ \sigma_t^2 =& 0.00116 + 0.1049 a_{t-1}^2 + 0.8178 \sigma_{t-1}^2 \tag{18.8} \end{align}\] 其中\(t^*_{\xi,d}\)表示标准化的有偏的t分布, 偏度参数为\(\xi\),自由度为\(d\)

输出的检验结果表明模型充分。 AIC为\(-1.451\), 相比对称的标准化t分布的\(-1.447\), 和正态分布时的AIC为\(-1.389\), 有偏t分布的AIC较优。

输出中偏度参数\(\xi\)的检验是\(H_0: \xi=0\)的检验, 但是我们需要的是\(H_0: \xi=1\)的检验。 计算\(z\)检验如下:

z <- (0.8717220 - 1) / 0.0629129
pv <- 2*(1 - pnorm(abs(z)))
c(statistic=z, pvalue=pv)
##   statistic      pvalue 
## -2.03897770  0.04145225

可以认为偏度参数\(\xi\)是显著小于1的, 模型中的条件分布是左偏的。

作标准化残差相对于\(t^*(0.87, 7.23)\)分布的QQ图:

plot(mod3, which=13)
标准化残差在有偏t分布假设时的QQ图

图18.7: 标准化残差在有偏t分布假设时的QQ图

从QQ图来看标准化残差比理论分布略重尾, 但左偏已经不明显。

18.1.6.4 讨论和比较

对Intel的月对数收益率建模, 同一采用了常数均值作为均值模型, GARCH(1,1)作为波动率模型, 但采用了三种不同的条件分布:

  • 正态分布;
  • t分布;
  • 有偏t分布。

从对\(\tilde a_t\)\(\tilde a_t^2\)的白噪声检验来看三个模型都很好地拟合了数据。 如果采用AIC来选择, 那么有偏t分布较优。 如果采用BIC来选择, 则会选择对称t分布, 实际上BIC倾向于参数较少的模型, 而且有偏t分布拟合得到的参数\(\xi\)与1(对称时)差距并不太大。 这也说明不同的比较准则可能会选择不同的模型。

将三个模型估计的波动率绘制在同一坐标系中进行比较, 可以看到三个模型拟合的波动率基本相同:

x.intel <- as.vector(time(ts.intel))
vola1 <- volatility(mod1)
vola2 <- volatility(mod2)
vola3 <- volatility(mod3)
matplot(x.intel, cbind(vola1, vola2, vola3),
        type="l",
        lty=1, col=c("green", "blue", "red"), 
        xlab="年", ylab="波动率")
legend("top", lty=1, col=c("green", "blue", "red"), 
       legend=c("标准正态分布", "对称标准化t分布", "左偏标准化t分布"))
Intel股票三种模型估计的波动率比较

图18.8: Intel股票三种模型估计的波动率比较

下面比较三个模型的波动率超前1到12步的预测结果, 预测基于截止到2009年12月的数据:

p1 <- predict(mod1, n.ahead=12)[["standardDeviation"]]
p2 <- predict(mod2, n.ahead=12)[["standardDeviation"]]
p3 <- predict(mod3, n.ahead=12)[["standardDeviation"]]
pred.tab <- tibble(
  "预测步数"=1:12,
  "正态"=p1,
  "对称t"=p2,
  "有偏t"=p3
)
knitr::kable(pred.tab, digits=4)
预测步数 正态 对称t 有偏t
1 0.0975 0.0951 0.0955
2 0.0993 0.0975 0.0979
3 0.1009 0.0997 0.1000
4 0.1023 0.1016 0.1019
5 0.1037 0.1034 0.1037
6 0.1050 0.1050 0.1053
7 0.1061 0.1065 0.1067
8 0.1072 0.1078 0.1080
9 0.1082 0.1090 0.1092
10 0.1092 0.1101 0.1103
11 0.1100 0.1111 0.1113
12 0.1109 0.1121 0.1122

三个模型的波动率多步预报基本相同。

18.1.6.5 预测的评估

波动率的三种定义都不是直接可观测的, 基本都来自模型估计。 所以, 比较不同波动率模型的预测结果具有挑战性。 有人用样本外预测进行评估, 比较波动率平方的预测值\(\sigma_h^2(\ell)\)与同期的扰动值\(a^2_{h+\ell}\), 但计算两者的相关系数一般是数值很小的。 这是因为虽然\(a_t^2\)\(\sigma_t^2\)的无偏估计, 但仅用单个点作为估计误差太大, 没有实际意义。

评估的更多信息参见Anderson和Bollerslev(1998)。

18.1.6.6 使用rugarch包

rugarch扩展包提供了更多的GARCH变种(改进)模型的估计功能, 我们用ugarch估计intel股票月度对数收益率的GARCH模型:

library(rugarch)
## Warning: 程辑包'rugarch'是用R版本4.0.4 来建造的
## 载入需要的程辑包:parallel
## 
## 载入程辑包:'rugarch'
## The following object is masked from 'package:purrr':
## 
##     reduce
## The following object is masked from 'package:stats':
## 
##     sigma
spec <- ugarchspec(
  mean.model = list(
    armaOrder=c(0,0),
    include.mean=TRUE  ),
  variance.model = list(
    model = "sGARCH", # standard GARCH model
    garchOrder = c(1,1) ) )
mod1ru <- ugarchfit(spec = spec, data = xts.intel)
show(mod1ru)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.011270    0.005393   2.0898 0.036639
## omega   0.000918    0.000389   2.3624 0.018158
## alpha1  0.085916    0.026449   3.2484 0.001160
## beta1   0.852779    0.039490  21.5949 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.011270    0.006457   1.7454 0.080920
## omega   0.000918    0.000451   2.0372 0.041634
## alpha1  0.085916    0.029592   2.9034 0.003692
## beta1   0.852779    0.042174  20.2204 0.000000
## 
## LogLikelihood : 312.343 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -1.3889
## Bayes        -1.3520
## Shibata      -1.3891
## Hannan-Quinn -1.3744
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.5185  0.4715
## Lag[2*(p+q)+(p+q)-1][2]    0.7006  0.6065
## Lag[4*(p+q)+(p+q)-1][5]    1.9298  0.6351
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.1328  0.7155
## Lag[2*(p+q)+(p+q)-1][5]    0.1423  0.9962
## Lag[4*(p+q)+(p+q)-1][9]    0.2784  0.9998
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]  0.001939 0.500 2.000  0.9649
## ARCH Lag[5]  0.015487 1.440 1.667  0.9991
## ARCH Lag[7]  0.025546 2.315 1.543  1.0000
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.381
## Individual Statistics:              
## mu     0.17241
## omega  0.06464
## alpha1 0.07993
## beta1  0.07156
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.07 1.24 1.6
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.1561 0.8760    
## Negative Sign Bias  0.4376 0.6619    
## Positive Sign Bias  0.1852 0.8531    
## Joint Effect        0.2314 0.9724    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     23.30       0.2245
## 2    30     30.73       0.3783
## 3    40     40.50       0.4038
## 4    50     47.44       0.5365
## 
## 
## Elapsed time : 0.189491

结果与18.1.6中用fGARCH包的garchFit()函数得到的结果基本一致, 输出的诊断统计也比较详细。默认的条件分布是正态分布(distribution.model="norm")。

改用条件t分布:

library(rugarch)
spec2 <- ugarchspec(
  mean.model = list(
    armaOrder=c(0,0),
    include.mean=TRUE  ),
  variance.model = list(
    model = "sGARCH", # standard GARCH model
    garchOrder = c(1,1) ),
  distribution.model="std" ) 
mod2ru <- ugarchfit(spec = spec2, data = xts.intel)
show(mod2ru)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : std 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.016508    0.005103   3.2349 0.001217
## omega   0.001156    0.000578   2.0002 0.045476
## alpha1  0.105542    0.037176   2.8390 0.004526
## beta1   0.817265    0.058094  14.0680 0.000000
## shape   6.787868    1.861988   3.6455 0.000267
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.016508    0.005690   2.9011 0.003718
## omega   0.001156    0.000571   2.0248 0.042886
## alpha1  0.105542    0.031871   3.3116 0.000928
## beta1   0.817265    0.054038  15.1238 0.000000
## shape   6.787868    1.846240   3.6766 0.000236
## 
## LogLikelihood : 326.2307 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -1.4470
## Bayes        -1.4009
## Shibata      -1.4472
## Hannan-Quinn -1.4288
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.5419  0.4617
## Lag[2*(p+q)+(p+q)-1][2]    0.7139  0.6011
## Lag[4*(p+q)+(p+q)-1][5]    1.8460  0.6552
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.2716  0.6023
## Lag[2*(p+q)+(p+q)-1][5]    0.3241  0.9811
## Lag[4*(p+q)+(p+q)-1][9]    0.4854  0.9986
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]  0.007562 0.500 2.000  0.9307
## ARCH Lag[5]  0.038613 1.440 1.667  0.9966
## ARCH Lag[7]  0.082546 2.315 1.543  0.9996
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.0878
## Individual Statistics:              
## mu     0.13157
## omega  0.10169
## alpha1 0.09333
## beta1  0.08261
## shape  0.05793
## 
## 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.34061 0.7336    
## Negative Sign Bias 0.29880 0.7652    
## Positive Sign Bias 0.02069 0.9835    
## Joint Effect       0.16520 0.9830    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     20.50       0.3648
## 2    30     32.22       0.3105
## 3    40     40.68       0.3961
## 4    50     44.29       0.6643
## 
## 
## Elapsed time : 0.1875019

结果与18.1.6.2中用fGARCH包计算的结果基本一致。

18.1.7 两步估计法

参考(18.2)对ARMA的模仿: \[ a_t^2 = \alpha_0 + \sum_{i=1}^{\max(m,s)} (\alpha_i + \beta_i) a_{t-i}^2 + \eta_t - \sum_{j=1}^s \beta_j \eta_{t-j} \] 提出如下的两步估计方法来估计GARCH模型。

第一,忽略ARCH效应, 用线性时间序列建模方法(如最大似然估计)对收益率序列建立均值方程。 残差序列用\(a_t\)表示。

第二, 将\(\{ a_t^2 \}\)作为观测序列, 可以用最大似然估计方法估计式(18.2)中的参数。 用\(\hat\phi_i\)\(\hat\theta_i\)分别表示AR和MA部分的系数估计值。 则GARCH模型参数估计为 \(\hat\beta_i = -\hat\theta_i\), \(\hat\alpha_i = \hat\phi_i + \hat\theta_i\)。 这种估计方法只是一种近似, 没有理论结果证明其合理性, 但是一些经验显示这样的估计往往能够给出不错的近似结果, 尤其是大中样本情形。

例如,对Intel月对数收益率序列, 先用常数作为均值方程,计算残差\(a_t\):

hatmu <- mean(ts.intel); hatmu
## [1] 0.0143273
at <- ts.intel - hatmu
at2 <- at^2

关于\(a_t^2\)用ARMA(1,1)估计:

mod4 <- arima(at2, order=c(1,0,1)); mod4
## 
## Call:
## arima(x = at2, order = c(1, 0, 1))
## 
## Coefficients:
##          ar1      ma1  intercept
##       0.9119  -0.7915     0.0161
## s.e.  0.0430   0.0635     0.0039
## 
## sigma^2 estimated as 0.001223:  log likelihood = 858.64,  aic = -1709.28

这里\(\hat\phi_1=0.9119\), \(\hat\theta_1 = -0.7915\), 所以\(\hat\beta_1 = 0.7915\), \(\hat\alpha_1 = 0.1204\)\[ \hat\alpha_0 = \hat\mu (1 - \hat\phi_1) = 0.0161(1 - 0.9919) = 0.00013 \] (注意arima()函数输出的intecept是模型均值而不是\(\phi_0\)参数)

这样的两步法得到的模型估计为: \[ r_t = 0.0143 + a_t, \quad a_t = \sigma_t \varepsilon_t, \quad \sigma_t^2 = 0.00013 + 0.1204 a_{t-1}^2 + 0.7915 \sigma_{t-1}^2 \]garchFit()得到的估计结果十分接近。

这样得到的误差波动率的估计为:

vola4 <- sqrt(at2 - mod4$residuals)

与有偏t分布的GARCH拟合的波动率计算相关系数:

cor(vola3, vola4)
## [1] 0.9976242

两者十分相似。

18.2 IGARCH模型

GARCH模型可以写成(18.2)这样的关于\(a_t^2\)的ARMA形式, 其中\(\eta_t = a_t^2 - \sigma_t^2\)是模型的扰动, 是鞅差白噪声:

\[\begin{align} a_t^2 = \alpha_0 + \sum_{i=1}^{\max(m,s)} (\alpha_i + \beta_i) a_{t-i}^2 + \eta_t - \sum_{j=1}^s \beta_j \eta_{t-j} \end{align}\]

如果这个模型中的AR部分有单位根(有特征根等于1), 则对应的模型不再满足GARCH模型条件, 称为IGARCH模型,或单位根GARCH模型。 类似于ARIMA模型, IGARCH模型的扰动\(\eta_t = a_t^2 - \sigma_t^2\)\(a_t^2\)的影响是持久的、不衰减的。

IGARCH(1,1)模型可以写成 \[ a_t = \sigma_t \varepsilon_t, \quad \sigma_t^2 = \alpha_0 + \beta_1 \sigma_{t-1}^2 + (1 - \beta_1) a_{t-1}^2 \]

蔡瑞胸教授的IGARCH(1,1)建模的程序参见§18.4

Intel股票的月对数收益率建立IGARCH(1,1)模型:

mod5 <- Igarch(as.vector(ts.intel))
## Estimates:  0.9217433 
## Maximized log-likehood:  -301.412 
## 
## Coefficient(s):
##       Estimate  Std. Error  t value   Pr(>|t|)    
## beta 0.9217433   0.0155534  59.2633 < 2.22e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

估计的模型为:

\[ \sigma_t^2 = 0.9217 \sigma_{t-1}^2 + 0.0793 a_{t-1}^2 \]

如果要估计\(\mu\)\(\alpha_0\)会导致参数估计最大似然估计计算失败。

IGARCH模型的\(a_t^2\)没有无条件方差。

用rugarch包估计IGARCH(1,1)模型:

library(rugarch)
speci <- ugarchspec(
  mean.model = list(
    armaOrder=c(0,0),
    include.mean=TRUE  ),
  variance.model = list(
    model = "iGARCH", # itegrated GARCH model
    garchOrder = c(1,1) ) )
mod5ru <- ugarchfit(spec = speci, data = xts.intel)
show(mod5ru)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : iGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.009665    0.005340   1.8100 0.070298
## omega   0.000345    0.000182   1.8990 0.057561
## alpha1  0.127533    0.033626   3.7927 0.000149
## beta1   0.872467          NA       NA       NA
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.009665    0.006240   1.5488 0.121425
## omega   0.000345    0.000191   1.8060 0.070925
## alpha1  0.127533    0.031820   4.0079 0.000061
## beta1   0.872467          NA       NA       NA
## 
## LogLikelihood : 308.1861 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -1.3747
## Bayes        -1.3470
## Shibata      -1.3748
## Hannan-Quinn -1.3638
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.6333  0.4261
## Lag[2*(p+q)+(p+q)-1][2]    0.8602  0.5459
## Lag[4*(p+q)+(p+q)-1][5]    1.8181  0.6620
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.2326  0.6296
## Lag[2*(p+q)+(p+q)-1][5]    0.5702  0.9463
## Lag[4*(p+q)+(p+q)-1][9]    0.8666  0.9911
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.1295 0.500 2.000  0.7190
## ARCH Lag[5]    0.2649 1.440 1.667  0.9495
## ARCH Lag[7]    0.3838 2.315 1.543  0.9877
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.1698
## Individual Statistics:              
## mu     0.22485
## omega  0.03416
## alpha1 0.21842
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          0.846 1.01 1.35
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias          0.24275 0.8083    
## Negative Sign Bias 0.06684 0.9467    
## Positive Sign Bias 0.25664 0.7976    
## Joint Effect       0.10895 0.9907    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     27.71      0.08914
## 2    30     37.08      0.14416
## 3    40     57.98      0.02571
## 4    50     59.38      0.14717
## 
## 
## Elapsed time : 0.1176879

rugarch包估计了omega参数, 即\(\alpha_0\)参数, 与蔡教授的做法不同。 这里固定了\(\beta_1 = 1 - \alpha_1\)

18.3 GARCH-M模型

有些金融资产的收益率的条件均值受到其波动率的影响,称为风险溢价。 GARCH-M模型就是用来描述这样的现象, M表示条件均值依赖于GARCH模型。 一种简单的GARCH-M(1,1)模型为 \[\begin{align} r_t =& \mu + c \sigma_t^2 + a_t, \quad a_t = \sigma_t \varepsilon_t,\\ \sigma_t^2 =& \alpha_0 + \alpha_1 a_{t-1}^2 + \beta_1 \sigma_{t-1}^2 \tag{18.9} \end{align}\] 模型中的收益率条件均值为 \(E(r_t | F_{t-1}) = \mu + c \sigma_t^2\) 需要用条件方差\(\sigma_t^2 = \text{Var}(r_t | F_{t-1})\)描述。 参数\(c\)称为风险溢价参数, 如果\(c\)为正值则收益率与波动率正相关。

文献中还有其他的风险溢价模型形式,如 \[\begin{aligned} r_t =& \mu + c \sigma_t + a_t \\ r_t =& \mu + c \ln\sigma_t^2 + a_t \end{aligned}\]

(18.9)中的收益率\(r_t\)不再是不相关列, 而是序列相关的, 相关性来自\(\sigma_t^2\)的序列相关性。 风险溢价的存在是股票收益率具有序列相关性的原因之一。

蔡瑞胸教授用来估计GARCH-M(1,1)模型的R函数参见§18.5

18.3.1 Intel股票月对数收益率的GARCH-M建模

mod6 <- garchM(100*as.vector(ts.intel))
## Warning in nlminb(start = params, objective = glkM, lower = lowBounds, upper =
## uppBounds): NA/NaN function evaluation
## Maximized log-likehood:  -1731.983 
## 
## Coefficient(s):
##         Estimate  Std. Error  t value   Pr(>|t|)    
## mu    0.07760752  1.32041495  0.05878  0.9531312    
## gamma 0.00794321  0.00918582  0.86473  0.3871896    
## omega 9.45891888  3.94051890  2.40042  0.0163761 *  
## alpha 0.08761597  0.02673484  3.27722  0.0010483 ** 
## beta  0.84933812  0.03948620 21.50975 < 2.22e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

估计的模型为:

\[\begin{aligned} 100 r_t =& 0.0776 + 0.0089 \sigma_t^2 + a_t, \quad a_t = \sigma_t \varepsilon_t \\ \sigma_t^2 =& 9.4589 + 0.0876 a_{t-1}^2 + 0.8493 \sigma_{t-1}^2 \end{aligned}\] 与教材结果不同, 可能是估计用的R程序经过升级。 估计中参数mu和参数gamma不显著, gamma不显著意味着没有GARCH-M效应。

用rugarch包估计GARCH-M模型:

library(rugarch)
specm <- ugarchspec(
  mean.model = list(
    armaOrder=c(0,0),
    include.mean=TRUE,
    archm=TRUE, archpow=2),
  variance.model = list(
    model = "sGARCH", # standard GARCH model
    garchOrder = c(1,1) ) )
mod6ru <- ugarchfit(spec = specm, data = xts.intel)
show(mod6ru)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error   t value Pr(>|t|)
## mu     -0.000181    0.013692 -0.013241 0.989436
## archm   0.855449    0.949115  0.901312 0.367422
## omega   0.000944    0.000392  2.409028 0.015995
## alpha1  0.087240    0.026657  3.272698 0.001065
## beta1   0.849721    0.039232 21.658619 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu     -0.000181    0.015260 -0.01188 0.990521
## archm   0.855449    1.097417  0.77951 0.435678
## omega   0.000944    0.000465  2.02972 0.042385
## alpha1  0.087240    0.029543  2.95293 0.003148
## beta1   0.849721    0.041844 20.30706 0.000000
## 
## LogLikelihood : 312.7599 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -1.3863
## Bayes        -1.3402
## Shibata      -1.3866
## Hannan-Quinn -1.3681
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.6479  0.4209
## Lag[2*(p+q)+(p+q)-1][2]    0.8671  0.5434
## Lag[4*(p+q)+(p+q)-1][5]    2.1293  0.5882
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.1373  0.7110
## Lag[2*(p+q)+(p+q)-1][5]    0.1496  0.9958
## Lag[4*(p+q)+(p+q)-1][9]    0.3173  0.9997
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3] 0.0002648 0.500 2.000  0.9870
## ARCH Lag[5] 0.0167638 1.440 1.667  0.9990
## ARCH Lag[7] 0.0530451 2.315 1.543  0.9999
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.5298
## Individual Statistics:              
## mu     0.15304
## archm  0.17085
## omega  0.07592
## alpha1 0.09511
## beta1  0.07486
## 
## 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.2498 0.8028    
## Negative Sign Bias  0.4723 0.6370    
## Positive Sign Bias  0.1415 0.8875    
## Joint Effect        0.2692 0.9657    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     17.35       0.5661
## 2    30     22.76       0.7876
## 3    40     30.77       0.8236
## 4    50     49.47       0.4544
## 
## 
## Elapsed time : 0.6173491

注意输入数据没有乘以100。 结果与蔡教授的程序结果基本一致, 其中的archm项不显著, 所以没有风险溢价效应。

18.3.2 标普500指数月超额收益率的GARCH-M建模

考虑标普500指数从1926年到1991年的月超额收益率。 共792个观测。

x.sp <- scan("sp500.txt", quiet=TRUE)*100
ts.sp <- ts(x.sp, start=c(1926,1), frequency = 12)
plot(ts.sp, main="标普500的月超额收益率(%)", xlab="年", ylab="超额收益率")
标普500的月超额收益率(%)

图18.9: 标普500的月超额收益率(%)

用GARCH(1,1)模型:

library(fGarch, quietly = TRUE)
mod7 <- garchFit(~ 1 + garch(1,1) , data=ts.sp, trace=FALSE)
summary(mod7)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~1 + garch(1, 1), data = ts.sp, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ 1 + garch(1, 1)
## <environment: 0x000000001ff9f878>
##  [data = ts.sp]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##      mu    omega   alpha1    beta1  
## 0.74497  0.80615  0.12198  0.85436  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu       0.74497     0.15377    4.845 1.27e-06 ***
## omega    0.80615     0.28333    2.845  0.00444 ** 
## alpha1   0.12198     0.02202    5.540 3.02e-08 ***
## beta1    0.85436     0.02175   39.276  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -2377.84    normalized:  -3.002323 
## 
## Description:
##  Tue Apr 21 11:10:43 2020 by user: user 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  80.32111  0           
##  Shapiro-Wilk Test  R    W      0.98505   3.136885e-07
##  Ljung-Box Test     R    Q(10)  11.2205   0.340599    
##  Ljung-Box Test     R    Q(15)  17.99703  0.262822    
##  Ljung-Box Test     R    Q(20)  24.29896  0.2295768   
##  Ljung-Box Test     R^2  Q(10)  9.920157  0.4475259   
##  Ljung-Box Test     R^2  Q(15)  14.21124  0.509572    
##  Ljung-Box Test     R^2  Q(20)  16.75081  0.6690903   
##  LM Arch Test       R    TR^2   13.04872  0.3655092   
## 
## Information Criterion Statistics:
##      AIC      BIC      SIC     HQIC 
## 6.014746 6.038355 6.014696 6.023820

模型的残差白噪声检验通过, 正态性假设不通过。得到的模型为 \[\begin{align} 100 r_t =& 0.7450 + a_t, \quad a_t = \sigma_t \varepsilon_t, \\ \sigma_t^2 =& 0.8061 + 0.1220 a_{t-1}^2 + 0.8544 \sigma_{t-1}^2 \end{align}\] 模型很接近IGARCH条件。

尝试对标普500超额收益率拟合GARCH-M(1,1)模型:

library(rugarch)
specm <- ugarchspec(
  mean.model = list(
    armaOrder=c(0,0),
    include.mean=TRUE,
    archm=TRUE, archpow=2),
  variance.model = list(
    model = "sGARCH", # standard GARCH model
    garchOrder = c(1,1) ) )
mod8ru <- ugarchfit(spec = specm, data = ts.sp)
show(mod8ru)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.542048    0.236018   2.2966 0.021639
## archm   0.010081    0.008885   1.1346 0.256536
## omega   0.829648    0.293212   2.8295 0.004662
## alpha1  0.123124    0.022286   5.5247 0.000000
## beta1   0.852261    0.022400  38.0479 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.542048    0.243967   2.2218 0.026296
## archm   0.010081    0.008913   1.1311 0.258024
## omega   0.829648    0.347874   2.3849 0.017083
## alpha1  0.123124    0.027802   4.4286 0.000009
## beta1   0.852261    0.027061  31.4946 0.000000
## 
## LogLikelihood : -2377.192 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       6.0156
## Bayes        6.0451
## Shibata      6.0156
## Hannan-Quinn 6.0270
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.7246  0.3946
## Lag[2*(p+q)+(p+q)-1][2]    0.7501  0.5869
## Lag[4*(p+q)+(p+q)-1][5]    3.0243  0.4026
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.7589  0.3837
## Lag[2*(p+q)+(p+q)-1][5]    1.8690  0.6497
## Lag[4*(p+q)+(p+q)-1][9]    3.4842  0.6771
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.5048 0.500 2.000  0.4774
## ARCH Lag[5]    1.0892 1.440 1.667  0.7066
## ARCH Lag[7]    1.6862 2.315 1.543  0.7833
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.327
## Individual Statistics:              
## mu     0.06463
## archm  0.03109
## omega  0.21419
## alpha1 0.26216
## beta1  0.16482
## 
## 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           3.5833 3.601e-04 ***
## Negative Sign Bias  1.1596 2.466e-01    
## Positive Sign Bias  0.6211 5.347e-01    
## Joint Effect       22.1805 5.983e-05 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     23.86      0.20163
## 2    30     41.11      0.06741
## 3    40     40.63      0.39861
## 4    50     56.86      0.20554
## 
## 
## Elapsed time : 0.481709
#garchM(as.vector(ts.sp))

拟合的模型为: \[\begin{align} \mu_t =& 0.54 + 0.010 \sigma_t^2 + a_t, \quad a_t = \sigma_t \varepsilon_t\\ \sigma_t^2 =& 0.83 + 0.12 a_{t-1}^2 + 0.85 \sigma_{t-1}^2 \end{align}\] 其中\(c=0.010\)不显著, 说明标普500的月超额收益率不存在GARCM-M现象, 没有风险溢价。

18.4 附录:蔡瑞胸教授的IGARCH建模估计R函数

"Igarch" <- function(rtn, include.mean=FALSE, volcnt=FALSE){
  # Estimation of a Gaussian IGARCH(1,1) model.
  # rtn: return series 
  # include.mean: flag for the constant in the mean equation.
  # volcnt: flag for the constant term of the volatility equation.
  #### default is the RiskMetrics model
  #
  Idata <<- rtn
  Flag <<- c(include.mean,volcnt)
  #
  Mean=mean(Idata); Var = var(Idata); S = 1e-6
  if((volcnt)&&(include.mean)){
    params=c(mu = Mean,omega=0.1*Var,beta=0.85)
    lowerBounds = c(mu = -10*abs(Mean), omega= S^2, beta= S)
    upperBounds = c(mu = 10*abs(Mean), omega = 100*Var, beta = 1-S)
  }
  if((volcnt)&&(!include.mean)){
    params=c(omega=0.1*Var, beta=0.85)
    lowerBounds=c(omega=S^2,beta=S)
    upperBounds=c(omega=100*Var,beta=1-S)
  }
  #
  if((!volcnt)&&(include.mean)){
    params=c(mu = Mean, beta= 0.8)
    lowerBounds = c(mu = -10*abs(Mean), beta= S)
    upperBounds = c(mu = 10*abs(Mean), beta = 1-S)
  }
  if((!volcnt)&&(!include.mean)){
    params=c(beta=0.85)
    lowerBounds=c(beta=S)
    upperBounds=c(beta=1-S)
  }
  # Step 3: set conditional distribution function:
  igarchDist = function(z,hh){dnorm(x = z/hh)/hh}
  # Step 4: Compose log-likelihood function:
  igarchLLH = function(parm){
    include.mean=Flag[1]
    volcnt=Flag[2]
    mu=0; omega = 0
    if((include.mean)&&(volcnt)){
      my=parm[1]; omega=parm[2]; beta=parm[3]}
    if((!include.mean)&&(volcnt)){
      omega=parm[1];beta=parm[2]}
    if((!include.mean)&&(!volcnt))beta=parm[1]
    if((include.mean)&&(!volcnt)){mu=parm[1]; beta=parm[2]}
    #
    z = (Idata - mu); Meanz = mean(z^2)
    e= omega + (1-beta)* c(Meanz, z[-length(Idata)]^2)
    h = stats::filter(e, beta, "r", init=Meanz)
    hh = sqrt(abs(h))
    llh = -sum(log(igarchDist(z, hh)))
    llh
  }
  # Step 5: Estimate Parameters and Compute Numerically Hessian:
  fit = nlminb(start = params, objective = igarchLLH,
               lower = lowerBounds, upper = upperBounds)
  ##lower = lowerBounds, upper = upperBounds, control = list(trace=3))
  epsilon = 0.0001 * fit$par
  cat("Estimates: ",fit$par,"\n")
  npar=length(params)
  Hessian = matrix(0, ncol = npar, nrow = npar)
  for (i in 1:npar) {
    for (j in 1:npar) {
      x1 = x2 = x3 = x4  = fit$par
      x1[i] = x1[i] + epsilon[i]; x1[j] = x1[j] + epsilon[j]
      x2[i] = x2[i] + epsilon[i]; x2[j] = x2[j] - epsilon[j]
      x3[i] = x3[i] - epsilon[i]; x3[j] = x3[j] + epsilon[j]
      x4[i] = x4[i] - epsilon[i]; x4[j] = x4[j] - epsilon[j]
      Hessian[i, j] = (igarchLLH(x1)-igarchLLH(x2)-igarchLLH(x3)+igarchLLH(x4))/
        (4*epsilon[i]*epsilon[j])
    }
  }
  cat("Maximized log-likehood: ",igarchLLH(fit$par),"\n")
  # Step 6: Create and Print Summary Report:
  se.coef = sqrt(diag(solve(Hessian)))
  tval = fit$par/se.coef
  matcoef = cbind(fit$par, se.coef, tval, 2*(1-pnorm(abs(tval))))
  dimnames(matcoef) = list(names(tval), c(" Estimate",
                                          " Std. Error", " t value", "Pr(>|t|)"))
  cat("\nCoefficient(s):\n")
  printCoefmat(matcoef, digits = 6, signif.stars = TRUE)
  
  if((include.mean)&&(volcnt)){
    mu=fit$par[1]; omega=fit$par[2]; beta = fit$par[3]
  }
  if((include.mean)&&(!volcnt)){
    mu = fit$par[1]; beta = fit$par[2]; omega = 0
  }
  if((!include.mean)&&(volcnt)){
    mu=0; omega=fit$par[1]; beta=fit$par[2]
  }
  if((!include.mean)&&(!volcnt)){
    mu=0; omega=0; beta=fit$par[1]
  }
  z=Idata-mu; Mz = mean(z^2)
  e= omega + (1-beta)*c(Mz,z[-length(z)]^2)
  h = stats::filter(e,beta,"r",init=Mz)
  vol = sqrt(abs(h))
  
  Igarch <- list(par=fit$par,volatility = vol)
}

18.5 附录:蔡瑞胸教授的GARCH-M建模估计R函数

"garchM" <- function(rtn, type=1){
  # Estimation of a Gaussian GARCH(1,1)-M model.
  ##### The program uses GARCH(1,1) results as initial values.
  # rtn: return series 
  # type = 1 for Variance-in-mean
  #      = 2 for volatility-in-mean
  #      = 3 for log(variance)-in-mean
  #
  if(is.matrix(rtn))rtn=c(rtn[,1])
  garchMdata <<- rtn
  # obtain initial estimates
  m1=garch11FIT(garchMdata)
  est=as.numeric(m1$par); v1=m1$ht  ## v1 is sigma.t-square
  Mean=est[1]; cc=est[2]; ar=est[3]; ma=est[4]; S=1e-6
  if(type==2)v1=sqrt(v1)
  if(type==3)v1=log(v1)
  #### Obtain initial estimate of the parameters for the mean equation
  m2=lm(rtn~v1)
  Cnst=as.numeric(m2$coefficients[1])
  gam=as.numeric(m2$coefficients[2])
  params=c(mu=Cnst,gamma=gam, omega=cc, alpha=ar,beta=ma)
  lowBounds=c(mu=-5*abs(Mean),gamma=-20*abs(gam), omega=S, alpha=S, beta=ma*0.6)
  uppBounds=c(mu=5*abs(Mean),gamma=100*abs(gam), omega=cc*5 ,alpha=3*ar,beta=1-S)
  ### Pass model information via defining global variable
  Vtmp <<- c(type,v1[1])
  #
  fit=nlminb(start = params, objective= glkM, lower=lowBounds, upper=uppBounds)
  ##,control=list(trace=3,rel.tol=1e-5))
  epsilon = 0.0001 * fit$par
  npar=length(params)
  Hessian = matrix(0, ncol = npar, nrow = npar)
  for (i in 1:npar) {
    for (j in 1:npar) {
      x1 = x2 = x3 = x4  = fit$par
      x1[i] = x1[i] + epsilon[i]; x1[j] = x1[j] + epsilon[j]
      x2[i] = x2[i] + epsilon[i]; x2[j] = x2[j] - epsilon[j]
      x3[i] = x3[i] - epsilon[i]; x3[j] = x3[j] + epsilon[j]
      x4[i] = x4[i] - epsilon[i]; x4[j] = x4[j] - epsilon[j]
      Hessian[i, j] = (glkM(x1)-glkM(x2)-glkM(x3)+glkM(x4))/
        (4*epsilon[i]*epsilon[j])
    }
  }
  cat("Maximized log-likehood: ",-glkM(fit$par),"\n")
  # Step 6: Create and Print Summary Report:
  se.coef = sqrt(diag(solve(Hessian)))
  tval = fit$par/se.coef
  matcoef = cbind(fit$par, se.coef, tval, 2*(1-pnorm(abs(tval))))
  dimnames(matcoef) = list(names(tval), c(" Estimate",
                                          " Std. Error", " t value", "Pr(>|t|)"))
  cat("\nCoefficient(s):\n")
  printCoefmat(matcoef, digits = 6, signif.stars = TRUE)
  
  m3=ResiVol(fit$par)
  
  garchM <- list(residuals=m3$residuals,sigma.t=m3$sigma.t)
}

glkM = function(pars){
  rtn <- garchMdata
  mu=pars[1]; gamma=pars[2]; omega=pars[3]; alpha=pars[4]; beta=pars[5]
  type=Vtmp[1]
  nT=length(rtn)
  # use conditional variance
  if(type==1){
    ht=Vtmp[2]
    et=rtn[1]-mu-gamma*ht
    at=c(et)
    for (i in 2:nT){
      sig2t=omega+alpha*at[i-1]^2+beta*ht[i-1]
      ept = rtn[i]-mu-gamma*sig2t
      at=c(at,ept)
      ht=c(ht,sig2t)
    }
  }
  # use volatility
  if(type==2){
    ht=Vtmp[2]^2
    et=rtn[1]-mu-gamma*Vtmp[2]
    at=c(et)
    for (i in 2:nT){
      sig2t=omega+alpha*at[i-1]^2+beta*ht[i-1]
      ept=rtn[i]-mu-gamma*sqrt(sig2t)
      at=c(at,ept)
      ht=c(ht,sig2t)
    }
  }
  # use log(variance)
  if(type==3){
    ht=exp(Vtmp[2])
    et=rtn[1]-mu-gamma*Vtmp[2]
    at=c(et)
    for (i in 2:nT){
      sig2t=omega+alpha*at[i-1]^2+beta*ht[i-1]
      ept=rtn[i]-mu-gamma*log(abs(sig2t))
      at=c(at,ept)
      ht=c(ht,sig2t)
    }
  }
  #
  hh=sqrt(abs(ht))
  glk=-sum(log(dnorm(x=at/hh)/hh))
  
  glk
}


ResiVol = function(pars){
  rtn <- garchMdata
  mu=pars[1]; gamma=pars[2]; omega=pars[3]; alpha=pars[4]; beta=pars[5]
  type=Vtmp[1]
  nT=length(rtn)
  # use conditional variance
  if(type==1){
    ht=Vtmp[2]
    et=rtn[1]-mu-gamma*ht
    at=c(et)
    for (i in 2:nT){
      sig2t=omega+alpha*at[i-1]^2+beta*ht[i-1]
      ept = rtn[i]-mu-gamma*sig2t
      at=c(at,ept)
      ht=c(ht,sig2t)
    }
  }
  # use volatility
  if(type==2){
    ht=Vtmp[2]^2
    et=rtn[1]-mu-gamma*Vtmp[2]
    at=c(et)
    for (i in 2:nT){
      sig2t=omega+alpha*at[i-1]^2+beta*ht[i-1]
      ept=rtn[i]-mu-gamma*sqrt(sig2t)
      at=c(at,ept)
      ht=c(ht,sig2t)
    }
  }
  # use log(variance)
  if(type==3){
    ht=exp(Vtmp[2])
    et=rtn[1]-mu-gamma*Vtmp[2]
    at=c(et)
    for (i in 2:nT){
      sig2t=omega+alpha*at[i-1]^2+beta*ht[i-1]
      ept=rtn[i]-mu-gamma*log(abs(sig2t))
      at=c(at,ept)
      ht=c(ht,sig2t)
    }
  }
  #
  
  ResiVol <- list(residuals=at,sigma.t=sqrt(ht))
}

garch11FIT = function(x){
  # Step 1: Initialize Time Series Globally:
  tx <<- x
  # Step 2: Initialize Model Parameters and Bounds:
  Mean = mean(tx); Var = var(tx); S = 1e-6
  params = c(mu = Mean, omega = 0.1*Var, alpha = 0.1, beta = 0.8)
  lowerBounds = c(mu = -10*abs(Mean), omega = S^2, alpha = S, beta = S)
  upperBounds = c(mu = 10*abs(Mean), omega = 100*Var, alpha = 1-S, beta = 1-S)
  # Step 3: Set Conditional Distribution Function:
  garchDist = function(z, hh) { dnorm(x = z/hh)/hh }
  # Step 4: Compose log-Likelihood Function:
  garchLLH = function(parm) {
    mu = parm[1]; omega = parm[2]; alpha = parm[3]; beta = parm[4]
    z = tx-mu; Mean = mean(z^2)
    # Use Filter Representation:
    e = omega + alpha * c(Mean, z[-length(tx)]^2)
    h = stats::filter(e, beta, "r", init = Mean)
    hh = sqrt(abs(h))
    llh = -sum(log(garchDist(z, hh)))
    llh }
  #####print(garchLLH(params))
  # Step 5: Estimate Parameters and Compute Numerically Hessian:
  fit = nlminb(start = params, objective = garchLLH,
               lower = lowerBounds, upper = upperBounds)
  #
  est=fit$par
  # compute the sigma.t^2 series
  z=tx-est[1]; Mean=mean(z^2)
  e=est[2]+est[3]*c(Mean,z[-length(tx)]^2)
  h=stats::filter(e,est[4],"r",init=Mean)
  
  garch11Fit <- list(par=est,ht=h)
}