时间序列分析之二 GARCH模型学习

Jean

《时间序列分析与R软件》教材上GARCH模型的例子，用tseries包只做了建模的部分，但该包的garch模型预测比较麻烦，然后试了fGarch包，ARMA+GARCH建模和预测都O.K.，问题是ARMA(p,q)当p>2时会出错（见stackoverflow帖子https://stackoverflow.com/questions/15475869/error-in-predict-for-arma-garch-model/15476326#15476326），不能运行预测函数，最后用rugarch包解决了问题。之前可用forecast包的auto.arima()自动定阶，用FinTS包的ArchTest()作ARCH集聚效应测试，方便了不少。R语言在时间序列分析上有丰富成熟的工具包，这是它的显著优势。

本例用道琼斯工业平均指数2006-04-20~2016-04-20之间的收盘指数计算每天的收益率，建模拟合及预测，修改了教材的例子，最后7天的数据用作预测测试，其余的用作样本建模。

金融量化分析我完全是小白，本文只是依样画葫芦，了解学习一下，纯学习贴。股指只覆盖了上市的大公司，而税务系统掌握了详实全面的一手经济收据，也许可以参考金融行业，开发一些经济指数反映经济运行的情况，如经济景气指数，行业景气指数等，提供各级政府与社会各界参考，这时候GARCH等金融模型分析就可以派上用场了。

先看结果，使用的是ARMA(3,2)+GARCH(1,1)模型，下图是rugarch包预测函数输出最后7天的预测结果，红线是预测值，黄色区域是标准差覆盖的范围。

把最后7天的实际值画上去，蓝色虚线之间正负标准差的范围基本覆盖了真实收益率的范围，模型的效果还是不错的。 下面详细讲一下建模分析的过程，数据来自astsa包的时间序列djia，先看看数据，包括了时间标签，开盘、高点、低点、收盘、成交量，是个xts时间序列。 用quantmod包的chartSeries()函数画出来看看。

chartSeries(djia2)

1、生成收益率时间序列，以收盘指数djia$Close计算每天的收益率。设X1,X2为相邻两天的收盘指数，收益率r=(X2-X1)/X1,当r不大时，r≈log(1+r)=log(X2/X1)=log(X2)-log(X1)，所以用diff()函数作差分运算即可，[-1]是删除中间的缺失值。共2157项，分成样本与测试数据集，画收益率的时序图观察。

djiar = 100*diff(log(djia$Close))[-1]

djiar1= djiar[2511:length(djiar)]

djiar2= djiar[1:2510]

2、平稳性检验。可以看到，它看起来象个均值为0的平稳序列，只是有些集中的波峰波谷。

adf.test(djiar2) # p-value = 0.01<0.05，平稳序列
pp.test(djiar2) # p-value = 0.01<0.05，平稳序列

kpss.test(djiar2, null = "Level") # p-value = 0.1>0.05，平稳序列

3、序列平稳先尝试ARMA模型建模，画自相关函数及偏自相关函数图，两个函数都不是严格在虚线内，也没有收敛的迹象，可能使用ARMA(p,q)模型。

acf2(djiar2)

4、确定模型阶数。现在有forecast包的auto.arima()函数自动根据AIC等准则定阶，比较方便，结果是ARMA(3,2)模型。

pdq = auto.arima(djiar2)

5、建模，ARMA模型不需要差分，所以中间的参数d=0。

(djiar.arma = arima(djiar2, order=c(3,0,2), include.mean = F))

djiar.arma

6、参数显著性检验，各参数原假设为0的t检验p值<<0.05，显著不为0，上下限区间[lwr,upr]不包含0，参数检验通过。 7、模型诊断，分析残差是否符合白噪声，模型不通过，不能用于预测。

tsdiag(djiar.arma)

A、残差零均值，但非等方差，范围超越了[-2,2]水平线。

B、残差自相关函数小且在2σ的虚线小范围内，接近于0，判定残差序列稳定。

C、Ljung-Box检验，P值有在虚线以内接近于0的，残差之间可能非独立。

8、分析原因。画残差平方的自相关函数，可以看到函数大部分在虚线外，衰减缓慢，说明残差的平方是自相关（非独立）的。

res = djiar.arma$residuals

acf(res^2)

画收益率序列平方的自相关函数与偏自相关函数，也看到了同样的自相关性质，说明存在条件异方差，方差受时间轴上前面方差的影响，不符合稳定序列ARMA模型等方差的要求。反映到收益率时序图上，就是波峰波谷的波动集聚性。

acf2(djiar2^2)

Engle（1982）注意到了收益率曲线中波动的集聚性并建立了ARCH模型，随后在金融业中得到了广泛的应用与验证，他因此获得了2003年的诺贝尔经济学奖。Bollerslev（1986）提出了广义ARCH（GARCH）模型有效减少了高阶下参数估计的数量，预测更准确。

稳定序列的ARMA模型，均值随时间变化产生了ARIMA模型，方差随时间变化产生了GARCH模型，是非平稳时间序列随机性分析的两面。下面看看GARCH建模的过程。

9、ARCH效应测试，FinTS包提供了ArchTest函数。

FinTS::ArchTest(x=djiar2)

p-value < 2.2e-16<<0.05，拒绝没有ARCH效应的原假设，有集聚效应，可以应用GARCH模型。

10、fGuarch包建模。ARMA建模已经处理了均值变化的部分，GARCH建模主要是处理残差的自相关部分，教材给出的是ARMA(2,3)+ARCH(6)模型，没有最后整合模型及预测验证。虽然预测有风险，但没有预测的建模是不完整的。fGarch及rugarch包都提供了ARMA及GARCH模型族组合建模与预测的支持，先看看fGarch包。

# 建模，ARMA(2,3)+GARCH(1,1)-SN模型

GARCH.model_1 = garchFit(~arma(2,3)+garch(1,1), data=djiar2, cond.dist='snorm',trace=FALSE)

GARCH.model_1

# 模型检验，检验残差的正态性与异方差及独立性，检验通过。

summary(GARCH.model_1)

参数p值全部<<0.05，显著性检验通过。残差Jarque-Bera Test p=0<<0.05，正态性检验不通过（后面详细分析）。Ljung-Box Test各阶p值全部>0.05，独立性通过。

# 预测通过

(pred.model_1 = predict(GARCH.model_1, n.ahead = 7, trace = TRUE, mse = 'cond', plot=FALSE))

meanForecast meanError standardDeviation
1 0.08216599 0.6285496 0.6285496

2 0.04002473 0.6397857 0.6381072

......

但当ARMA(p,q)的p>2时，predict()函数执行出错，原因已在文章开头交代。

11、改用rugarch包建模，教材说大多数情况下GARCH(1,1)即可，拟合不充分再尝试更高阶。

# 模型参数，分为三部分

# variance.model是方差方程，GARCH(1,1) sGARCH模型

# mean.model是均值方程，ARMA(3,2)模型

# distribution.model，方差条件分布函数，默认是norm正态分布

garch11.spec = ugarchspec(variance.model =list(model="sGARCH",garchOrder=c(1,1)), mean.model = list(armaOrder=c(3,2)),distribution.model = "norm")
garch11.fit = ugarchfit(spec=garch11.spec, data=djiar2)

12、参数显著性检验。参数t检验p值全部<0.05，显著不为0，检验通过。模型拟合结果解读可以参考https://bbs.pinggu.org/thread-3069231-1-1.html 。

garch11.fit

13、模型诊断，残差符合白噪声，诊断通过，可以用于预测。

garch11.fit

A、ARCH LM Tests p值>0.05，残差已不存在异方差（无ARCH效应），即符合方差齐性。
B、Ljung-Box Test p值>0.05，残差无自相关，残差独立。

C、标准化残差自相关函数基本在2σ虚线内，残差序列稳定。

plot(garch11.fit)

D、残差平方的自相关函数图基本在在2σ虚线内，残差平方不存在自相关，GARCH(1,1)阶数已够。 14、残差的正态性检验。白噪声不要求残差的正态性，但ARCH模型要求。plot(garch11.fit)有12个图可以提供模型诊断，比fGarch方便多了，画QQ图可知残差不符合正态性，Jarque-Bera测试也拒绝了正态性假设。

res = residuals(garch11.fit)
jarqueberaTest(res) # p Value: < 2.2e-16<<0.05

ARCH模型的重要假设是随机项服从条件正态分布，εt~ N(0,ht)，其中ht是前面时间节点随机项的多项式。教材中ARCH模型诊断也是要符合正态性假设的，通过Jarque-Bera测试进行。教材中用ARMA(2,3)+ARCH(6)模型，残差也不是正态分布，最终通过了上述模型，但没有解释具体的理由。

标准化残差经验分布曲线显示非正态分布，plot(garch11.fit)选第8项画出。

在阅读人大版《应用时间序列分析（第5版）》第7章《条件异方差模型》（7.5.5拟合检验）后，解决了这个问题，人大的金融实战经验还是比较丰富。该教材指出，实践中很多金融时间序列具有尖峰厚尾特征，峰态系数特别大，Jarque-Bera检验通常会拒绝正态分布，此时GARCH模型通常能提取条件异方差信息，但不一定能通过白噪声检验。是否容忍不服从正态分布的假定，取决于2个条件：一是模型拟合的精度是否满意，满意就可以；二是是否有更适当的分布可以替换正态分布，这导致了EGARCH、TGARCH等GARCH模型族的衍生。

从上图标准化残差的经验分布曲线可以看出其尖峰厚尾特征显著，下面看看其拟合的精度，在plot(garch11.fit)命令中选择第1项“Series with 2 Conditional SD Superimposed”，看看拟合的精度。

蓝色线是观察值，水平虚线是无条件方差下95%置信区间（2σ），红色线是条件方差下95%置信区间，可见条件方差模型GARCH(1,1)拟合的效果不错，可以接受残差的非正态分布。

15、预测及查看结果，这里还没有使用滚动预测，n.roll=0。

garch11.fcst = ugarchforecast(garch11.fit, n.ahead=7)

plot(garch11.fcst)

16、画预测值与实际值的比较图，可以看到预测值大部分在正负标准差σ覆盖的范围之内，模型拟合和预测的效果还是不错的。 后续测试了EGARCH（指数GARCH）及TGARCH（门限GARCH）模型，AIC等统计量数据显示TGARCH模型统计量数值要小一点，效果会稍好一点，建模调用时改个参数就可以了。

tgarch11.spec = ugarchspec(variance.model = list(model="fGARCH", submodel="TGARCH", garchOrder=c(1,1)),　
mean.model = list(armaOrder=c(3,2)),
distribution.model = "norm")
tgarch11.fit = ugarchfit(spec=tgarch11.spec, data=djiar2)