R语言股票市场指数：ARMA-GARCH模型和对数收益率数据探索性分析（中）：https://developer.aliyun.com/article/1491668
峰度
有正峰度的年份是：
## \[1\] "2007" "2008" "2009" "2010" "2011" "2012" "2013" "2014" "2015" "2016" ## \[11\] "2017" "2018"
均按升序排列。
## 2010 2009 2008 2017 2018 2016 2013 ## Kurtosis 1.353797 1.500979 2.616615 3.411036 4.335748 4.647785 4.68112 ## 2015 2007 2012 2014 2011 ## Kurtosis 4.754926 5.345212 8.115847 9.850061 14.55464
箱形图
可以在2011、2014和2016年发现正的极端值。在2007、2011、2012、2014年可以发现负的极端值。
密度图
shapiro检验
## result ## 2007 3.695053e-09 ## 2008 6.160136e-07 ## 2009 2.083475e-04 ## 2010 1.500060e-03 ## 2011 3.434415e-18 ## 2012 8.417627e-12 ## 2013 1.165184e-10 ## 2014 1.954662e-16 ## 2015 5.261037e-11 ## 2016 7.144940e-11 ## 2017 1.551041e-08 ## 2018 3.069196e-09
基于报告的p值，我们可以拒绝所有正态分布的零假设。
QQ图
在所有报告的年份都可以发现偏离正态状态。
对数收益率GARCH模型
我将为工业平均指数（DJIA）的每日对数收益率建立一个ARMA-GARCH模型。
这是工业平均指数每日对数收益的图。
plot(ret)
离群值检测
Performance Analytics程序包中的Return.clean函数能够清除异常值。在下面，我们将原始时间序列与调整离群值后的进行比较。
clean(ret, "boudt")
作为对波动率评估的更为保守的方法，本文将以原始时间序列进行分析。
相关图
以下是自相关和偏相关图。
acf(ret)
pacf(dj_ret)
上面的相关图表明p和q> 0的一些ARMA（p，q）模型。将在本分析的该范围内对此进行验证。
单位根检验
我们运行Augmented Dickey-Fuller检验。
## ## ############################################### ## # Augmented Dickey-Fuller Test Unit Root Test # ## ############################################### ## ## Test regression none ## ## ## Call: ## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag) ## ## Residuals: ## Min 1Q Median 3Q Max ## -0.081477 -0.004141 0.000762 0.005426 0.098777 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## z.lag.1 -1.16233 0.02699 -43.058 < 2e-16 *** ## z.diff.lag 0.06325 0.01826 3.464 0.000539 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.01157 on 2988 degrees of freedom ## Multiple R-squared: 0.5484, Adjusted R-squared: 0.5481 ## F-statistic: 1814 on 2 and 2988 DF, p-value: < 2.2e-16 ## ## ## Value of test-statistic is: -43.0578 ## ## Critical values for test statistics: ## 1pct 5pct 10pct ## tau1 -2.58 -1.95 -1.62
基于报告的检验统计数据与临界值的比较，我们拒绝单位根存在的零假设。
ARMA模型
现在，我们确定时间序列的ARMA结构，以便对结果残差进行ARCH效应检验。ACF和PACF系数拖尾表明存在ARMA（2,2）。我们利用auto.arima（）函数开始构建。
## Series: ret ## ARIMA(2,0,4) with zero mean ## ## Coefficients: ## ar1 ar2 ma1 ma2 ma3 ma4 ## 0.4250 -0.8784 -0.5202 0.8705 -0.0335 -0.0769 ## s.e. 0.0376 0.0628 0.0412 0.0672 0.0246 0.0203 ## ## sigma^2 estimated as 0.0001322: log likelihood=9201.19 ## AIC=-18388.38 AICc=-18388.34 BIC=-18346.29 ## ## Training set error measures: ## ME RMSE MAE MPE MAPE MASE ## Training set 0.0002416895 0.01148496 0.007505056 NaN Inf 0.6687536 ## ACF1 ## Training set -0.002537238
建议使用ARMA（2,4）模型。但是，ma3系数在统计上并不显着，进一步通过以下方法验证：
## z test of coefficients: ## ## Estimate Std. Error z value Pr(>|z|) ## ar1 0.425015 0.037610 11.3007 < 2.2e-16 *** ## ar2 -0.878356 0.062839 -13.9779 < 2.2e-16 *** ## ma1 -0.520173 0.041217 -12.6204 < 2.2e-16 *** ## ma2 0.870457 0.067211 12.9511 < 2.2e-16 *** ## ma3 -0.033527 0.024641 -1.3606 0.1736335 ## ma4 -0.076882 0.020273 -3.7923 0.0001492 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
因此，我们将MA阶q <= 2作为约束。
## Series: dj_ret ## ARIMA(2,0,2) with zero mean ## ## Coefficients: ## ar1 ar2 ma1 ma2 ## -0.5143 -0.4364 0.4212 0.3441 ## s.e. 0.1461 0.1439 0.1512 0.1532 ## ## sigma^2 estimated as 0.0001325: log likelihood=9196.33 ## AIC=-18382.66 AICc=-18382.64 BIC=-18352.6 ## ## Training set error measures: ## ME RMSE MAE MPE MAPE MASE ## Training set 0.0002287171 0.01150361 0.007501925 Inf Inf 0.6684746 ## ACF1 ## Training set -0.002414944
现在，所有系数都具有统计意义。
## z test of coefficients: ## ## Estimate Std. Error z value Pr(>|z|) ## ar1 -0.51428 0.14613 -3.5192 0.0004328 *** ## ar2 -0.43640 0.14392 -3.0322 0.0024276 ** ## ma1 0.42116 0.15121 2.7853 0.0053485 ** ## ma2 0.34414 0.15323 2.2458 0.0247139 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
使用ARMA（2,1）和ARMA（1,2）进行的进一步验证得出的AIC值高于ARMA（2,2）。因此，ARMA（2,2）是更可取的。这是结果。
## Series: dj_ret ## ARIMA(2,0,1) with zero mean ## ## Coefficients: ## ar1 ar2 ma1 ## -0.4619 -0.1020 0.3646 ## s.e. 0.1439 0.0204 0.1438 ## ## sigma^2 estimated as 0.0001327: log likelihood=9194.1 ## AIC=-18380.2 AICc=-18380.19 BIC=-18356.15 ## ## Training set error measures: ## ME RMSE MAE MPE MAPE MASE ## Training set 0.0002370597 0.01151213 0.007522059 Inf Inf 0.6702687 ## ACF1 ## Training set 0.0009366271 coeftest(auto_model3) ## ## z test of coefficients: ## ## Estimate Std. Error z value Pr(>|z|) ## ar1 -0.461916 0.143880 -3.2104 0.001325 ** ## ar2 -0.102012 0.020377 -5.0062 5.552e-07 *** ## ma1 0.364628 0.143818 2.5353 0.011234 * ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
所有系数均具有统计学意义。
## ARIMA(1,0,2) with zero mean ## ## Coefficients: ## ar1 ma1 ma2 ## -0.4207 0.3259 -0.0954 ## s.e. 0.1488 0.1481 0.0198 ## ## sigma^2 estimated as 0.0001328: log likelihood=9193.01 ## AIC=-18378.02 AICc=-18378 BIC=-18353.96 ## ## Training set error measures: ## ME RMSE MAE MPE MAPE MASE ## Training set 0.0002387398 0.0115163 0.007522913 Inf Inf 0.6703448 ## ACF1 ## Training set -0.001958194 coeftest(auto_model4) ## ## z test of coefficients: ## ## Estimate Std. Error z value Pr(>|z|) ## ar1 -0.420678 0.148818 -2.8268 0.004702 ** ## ma1 0.325918 0.148115 2.2004 0.027776 * ## ma2 -0.095407 0.019848 -4.8070 1.532e-06 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
所有系数均具有统计学意义。此外，我们使用TSA软件包报告中的eacf（）函数。
## AR/MA ## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 ## 0 x x x o x o o o o o o o o x ## 1 x x o o x o o o o o o o o o ## 2 x o o x x o o o o o o o o o ## 3 x o x o x o o o o o o o o o ## 4 x x x x x o o o o o o o o o ## 5 x x x x x o o x o o o o o o ## 6 x x x x x x o o o o o o o o ## 7 x x x x x o o o o o o o o o
以“ O”为顶点的左上三角形位于（p，q）= {（1,2 ,,（2,2），（1,3）}}内，它表示一组潜在候选对象（p，q）值。ARMA（1,2）模型已经过验证。ARMA（2,2）已经是候选模型。让我们验证ARMA（1,3）。
## Call: ## ## Coefficients: ## ar1 ma1 ma2 ma3 ## -0.2057 0.1106 -0.0681 0.0338 ## s.e. 0.2012 0.2005 0.0263 0.0215 ## ## sigma^2 estimated as 0.0001325: log likelihood = 9193.97, aic = -18379.94 coeftest(arima_model5) ## ## z test of coefficients: ## ## Estimate Std. Error z value Pr(>|z|) ## ar1 -0.205742 0.201180 -1.0227 0.306461 ## ma1 0.110599 0.200475 0.5517 0.581167 ## ma2 -0.068124 0.026321 -2.5882 0.009647 ** ## ma3 0.033832 0.021495 1.5739 0.115501 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
只有一个系数具有统计意义。
结论是，我们选择ARMA（2,2）作为均值模型。现在，我们可以继续进行ARCH效果检验。
ARCH效应检验
现在，我们可以检验模型残差上是否存在ARCH效应。如果ARCH效应对于我们的时间序列的残差在统计上显着，则需要GARCH模型。
## ARCH LM-test; Null hypothesis: no ARCH effects ## ## data: model\_residuals - mean(model\_residuals) ## Chi-squared = 986.82, df = 12, p-value < 2.2e-16
基于报告的p值，我们拒绝没有ARCH效应的原假设。
让我们看一下残差相关图。
条件波动率
条件均值和方差定义为：
μt：= E（rt | Ft-1）σt2：= Var（rt | Ft-1）= E [（rt-μt）2 | Ft-1]
条件波动率可以计算为条件方差的平方根。
eGARCH模型
将sGARCH作为方差模型的尝试未获得具有统计显着性系数的结果。而指数GARCH（eGARCH）方差模型能够捕获波动率内的不对称性。要检查DJIA对数收益率内的不对称性，显示汇总统计数据和密度图。
## DAdjusted ## nobs 3019.000000 ## NAs 0.000000 ## Minimum -0.082005 ## Maximum 0.105083 ## 1. Quartile -0.003991 ## 3. Quartile 0.005232 ## Mean 0.000207 ## Median 0.000551 ## Sum 0.625943 ## SE Mean 0.000211 ## LCL Mean -0.000206 ## UCL Mean 0.000621 ## Variance 0.000134 ## Stdev 0.011593 ## Skewness -0.141370 ## Kurtosis 10.200492
负偏度值确认分布内不对称性的存在。
这给出了密度图。
我们继续提出eGARCH模型作为方差模型（针对条件方差）。更准确地说，我们将使用ARMA（2,2）作为均值模型，指数GARCH（1,1）作为方差模型对ARMA-GARCH进行建模。
在此之前，我们进一步强调ARMA（0,0）在这种情况下不令人满意。ARMA-GARCH：ARMA（0,0）+ eGARCH（1,1）
## ## *---------------------------------* ## * GARCH Model Fit * ## *---------------------------------* ## ## Conditional Variance Dynamics ## ----------------------------------- ## GARCH Model : eGARCH(1,1) ## Mean Model : ARFIMA(0,0,0) ## Distribution : sstd ## ## Optimal Parameters ## ------------------------------------ ## Estimate Std. Error t value Pr(>|t|) ## mu 0.000303 0.000117 2.5933 0.009506 ## omega -0.291302 0.016580 -17.5699 0.000000 ## alpha1 -0.174456 0.013913 -12.5387 0.000000 ## beta1 0.969255 0.001770 547.6539 0.000000 ## gamma1 0.188918 0.021771 8.6773 0.000000 ## skew 0.870191 0.021763 39.9848 0.000000 ## shape 6.118380 0.750114 8.1566 0.000000 ## ## Robust Standard Errors: ## Estimate Std. Error t value Pr(>|t|) ## mu 0.000303 0.000130 2.3253 0.020055 ## omega -0.291302 0.014819 -19.6569 0.000000 ## alpha1 -0.174456 0.016852 -10.3524 0.000000 ## beta1 0.969255 0.001629 595.0143 0.000000 ## gamma1 0.188918 0.031453 6.0063 0.000000 ## skew 0.870191 0.022733 38.2783 0.000000 ## shape 6.118380 0.834724 7.3298 0.000000 ## ## LogLikelihood : 10138.63 ## ## Information Criteria ## ------------------------------------ ## ## Akaike -6.7119 ## Bayes -6.6980 ## Shibata -6.7119 ## Hannan-Quinn -6.7069 ## ## Weighted Ljung-Box Test on Standardized Residuals ## ------------------------------------ ## statistic p-value ## Lag\[1\] 5.475 0.01929 ## Lag\[2*(p+q)+(p+q)-1\]\[2\] 6.011 0.02185 ## Lag\[4*(p+q)+(p+q)-1\]\[5\] 7.712 0.03472 ## d.o.f=0 ## H0 : No serial correlation ## ## Weighted Ljung-Box Test on Standardized Squared Residuals ## ------------------------------------ ## statistic p-value ## Lag\[1\] 1.342 0.2467 ## Lag\[2*(p+q)+(p+q)-1\]\[5\] 2.325 0.5438 ## Lag\[4*(p+q)+(p+q)-1\]\[9\] 2.971 0.7638 ## d.o.f=2 ## ## Weighted ARCH LM Tests ## ------------------------------------ ## Statistic Shape Scale P-Value ## ARCH Lag\[3\] 0.3229 0.500 2.000 0.5699 ## ARCH Lag\[5\] 1.4809 1.440 1.667 0.5973 ## ARCH Lag\[7\] 1.6994 2.315 1.543 0.7806 ## ## Nyblom stability test ## ------------------------------------ ## Joint Statistic: 4.0468 ## Individual Statistics: ## mu 0.2156 ## omega 1.0830 ## alpha1 0.5748 ## beta1 0.8663 ## gamma1 0.3994 ## skew 0.1044 ## shape 0.4940 ## ## Asymptotic Critical Values (10% 5% 1%) ## Joint Statistic: 1.69 1.9 2.35 ## Individual Statistic: 0.35 0.47 0.75 ## ## Sign Bias Test ## ------------------------------------ ## t-value prob sig ## Sign Bias 1.183 0.23680 ## Negative Sign Bias 2.180 0.02932 ** ## Positive Sign Bias 1.554 0.12022 ## Joint Effect 8.498 0.03677 ** ## ## ## Adjusted Pearson Goodness-of-Fit Test: ## ------------------------------------ ## group statistic p-value(g-1) ## 1 20 37.24 0.00741 ## 2 30 42.92 0.04633 ## 3 40 52.86 0.06831 ## 4 50 65.55 0.05714 ## ## ## Elapsed time : 0.6527421
所有系数均具有统计学意义。但是，根据以上报告的p值的标准化残差加权Ljung-Box检验，我们确认该模型无法捕获所有ARCH效果（我们拒绝了残差内无相关性的零假设） ）。
作为结论，我们通过在下面所示的GARCH拟合中指定ARMA（2,2）作为均值模型来继续进行。
ARMA-GARCH：ARMA（2,2）+ eGARCH（1,1）
## ## *---------------------------------* ## * GARCH Model Fit * ## *---------------------------------* ## ## Conditional Variance Dynamics ## ----------------------------------- ## GARCH Model : eGARCH(1,1) ## Mean Model : ARFIMA(2,0,2) ## Distribution : sstd ## ## Optimal Parameters ## ------------------------------------ ## Estimate Std. Error t value Pr(>|t|) ## ar1 -0.47642 0.026115 -18.2433 0 ## ar2 -0.57465 0.052469 -10.9523 0 ## ma1 0.42945 0.025846 16.6157 0 ## ma2 0.56258 0.054060 10.4066 0 ## omega -0.31340 0.003497 -89.6286 0 ## alpha1 -0.17372 0.011642 -14.9222 0 ## beta1 0.96598 0.000027 35240.1590 0 ## gamma1 0.18937 0.011893 15.9222 0 ## skew 0.84959 0.020063 42.3469 0 ## shape 5.99161 0.701313 8.5434 0 ## ## Robust Standard Errors: ## Estimate Std. Error t value Pr(>|t|) ## ar1 -0.47642 0.007708 -61.8064 0 ## ar2 -0.57465 0.018561 -30.9608 0 ## ma1 0.42945 0.007927 54.1760 0 ## ma2 0.56258 0.017799 31.6074 0 ## omega -0.31340 0.003263 -96.0543 0 ## alpha1 -0.17372 0.012630 -13.7547 0 ## beta1 0.96598 0.000036 26838.0412 0 ## gamma1 0.18937 0.013003 14.5631 0 ## skew 0.84959 0.020089 42.2911 0 ## shape 5.99161 0.707324 8.4708 0 ## ## LogLikelihood : 10140.27 ## ## Information Criteria ## ------------------------------------ ## ## Akaike -6.7110 ## Bayes -6.6911 ## Shibata -6.7110 ## Hannan-Quinn -6.7039 ## ## Weighted Ljung-Box Test on Standardized Residuals ## ------------------------------------ ## statistic p-value ## Lag\[1\] 0.03028 0.8619 ## Lag\[2*(p+q)+(p+q)-1\]\[11\] 5.69916 0.6822 ## Lag\[4*(p+q)+(p+q)-1\]\[19\] 12.14955 0.1782 ## d.o.f=4 ## H0 : No serial correlation ## ## Weighted Ljung-Box Test on Standardized Squared Residuals ## ------------------------------------ ## statistic p-value ## Lag\[1\] 1.666 0.1967 ## Lag\[2*(p+q)+(p+q)-1\]\[5\] 2.815 0.4418 ## Lag\[4*(p+q)+(p+q)-1\]\[9\] 3.457 0.6818 ## d.o.f=2 ## ## Weighted ARCH LM Tests ## ------------------------------------ ## Statistic Shape Scale P-Value ## ARCH Lag\[3\] 0.1796 0.500 2.000 0.6717 ## ARCH Lag\[5\] 1.5392 1.440 1.667 0.5821 ## ARCH Lag\[7\] 1.6381 2.315 1.543 0.7933 ## ## Nyblom stability test ## ------------------------------------ ## Joint Statistic: 4.4743 ## Individual Statistics: ## ar1 0.07045 ## ar2 0.37070 ## ma1 0.07702 ## ma2 0.39283 ## omega 1.00123 ## alpha1 0.49520 ## beta1 0.79702 ## gamma1 0.51601 ## skew 0.07163 ## shape 0.55625 ## ## Asymptotic Critical Values (10% 5% 1%) ## Joint Statistic: 2.29 2.54 3.05 ## Individual Statistic: 0.35 0.47 0.75 ## ## Sign Bias Test ## ------------------------------------ ## t-value prob sig ## Sign Bias 0.4723 0.63677 ## Negative Sign Bias 1.7969 0.07246 * ## Positive Sign Bias 2.0114 0.04438 ** ## Joint Effect 7.7269 0.05201 * ## ## ## Adjusted Pearson Goodness-of-Fit Test: ## ------------------------------------ ## group statistic p-value(g-1) ## 1 20 46.18 0.0004673 ## 2 30 47.73 0.0156837 ## 3 40 67.07 0.0034331 ## 4 50 65.51 0.0574582 ## ## ## Elapsed time : 0.93679
所有系数均具有统计学意义。在标准化残差或标准化平方残差内未发现相关性。模型正确捕获所有ARCH效果。然而：
*对于某些模型参数，Nyblom稳定性检验无效假设认为模型参数随时间是恒定的
*正偏差为零的假设在5％的显着性水平上被拒绝；这种检验着重于正面冲击的影响
*拒绝了标准化残差的经验和理论分布相同的Pearson拟合优度检验原假设
_注意_：ARMA（1,2）+ eGARCH（1,1）拟合还提供统计上显着的系数，标准化残差内没有相关性，标准化平方残差内没有相关性，并且正确捕获了所有ARCH效应。但是，偏差检验在5％时不如ARMA（2,2）+ eGARCH（1,1）模型令人满意。
进一步显示诊断图。
我们用平均模型拟合（红线）和条件波动率（蓝线）显示了原始的对数收益时间序列。
p <- addSeries(mean\_model\_fit, col = 'red', on = 1) p <- addSeries(cond_volatility, col = 'blue', on = 1) p
模型方程式
结合ARMA（2,2）和eGARCH模型，我们可以：
yt − ϕ1yt−1 − ϕ2yt−2 = ϕ0 + ut + θ1ut−1 +θ2ut-2ut= σtϵt，ϵt = N（0,1）ln（σt2）=ω+ ∑j = 1q（αjϵt−j2 +γ （ϵt−j–E | ϵt−j |））+ ∑i =1pβiln（σt−12）
使用模型结果系数，结果如下。
yt +0.476 yt-1 +0.575 yt-2 = ut +0.429 ut-1 +0.563 ut-2ut = σtϵt，ϵt = N（0,1）ln（σt2）= -0.313 -0.174ϵt-12 +0.189（ ϵt−1–E | ϵt−1 |））+ 0.966 ln（σt−12）
波动率分析
这是由ARMA（2,2）+ eGARCH（1,1）模型得出的条件波动图。
plot(cond_volatility)
显示了年条件波动率的线线图。
pl <- lapply(2007:2018, function(x) { plot(cond_volatility\[as.character(x)\]) pl
显示了按年列出的条件波动率箱图。
2008年之后，日波动率基本趋于下降。在2017年，波动率低于其他任何年。不同的是，与2017年相比，我们在2018年的波动性显着增加。
