本文是时间序列分析课程的作业，基于R、Rnw和Latex进行编写。
GARCH代码实现主要参考自《经济与金融计量方法：原理、应用案例及R语言实现》和对应包的官方文档，代码进一步整合，但每次执行时可能需要较长的时间，建议执行完后将结果导出成excel。如果本文存在问题，随时欢迎交流~

一、数据来源

沪深300指数，是由沪深证券交易所于 2005 年 4 月 8 日联合发布的反映沪深 300 指数编制目标和运行状况的金融指标，并能够作为投资业绩的评价标准，为指数化投资和指数衍生产品创新提供基础条件。因此，本次数据来源于网易财经，研究的数据集对象是沪深 300 指数（股票代码为000300），此次分析选取了沪深 300 指数2000 年1月到2019年12月的工作日收盘价格数据。

二、数据分析

（一）时序图

为了分析数据的波动情况，对其进行对数化和差分得到对数收益率，下图为沪深 300 指数的收盘价时序图和对数收益率时序图。

（二）平稳性检验

Augmented Dickey-Fuller Test (ADF) 是 DF 检验的拓展形式，可以对存在高阶滞后的序列进行单位根检验，原假设存在单位根，即序列不平稳。本文使用adf.test()进行单位根检验，检验结果如下所示，p 值远小于 0.01 说明拒绝原假设，即序列是平稳的。

三、模型建立

（一）均值模型

1. 均值模型的识别

序列平稳后，使用auto.arima()对序列自动识别均值模型。

识别出来的模型为ARMA(4, 4)。经过模型识别后，对模型ARMA(4,4)ARMA(4,4)ARMA(4,4)进行参数显著性检验。检验结果发现部分参数不显著，采用建立疏系数的均值模型，将不显著的参数强制为0。
采用疏系数的ARMAARMAARMA模型后的参数显著性检验结果来看，非零参数结果都显著。（注：这一步存有疑问，后续依旧是采用ARMA(4,4)ARMA(4,4)ARMA(4,4)的非疏系数模型进行拟合）

（二）方差模型建立

1. ARCH效应检验

上述建立模型后，对残差进行ARCH效应检验。Ljung-Box统计量Q(m)Q(m)Q(m)对残差序列进行自相关检验。原假设是序列不存在自相关，在残差的平方序列中可以检验条件异方差。

使用MTS包中的archTest()进行检验
检验结果显示，滞后10阶和滞后20阶的残差序列存在自相关，因此拒绝原假设，残差序列存在ARCH效应。

## [1] "m = 10"
## Q(m) of squared series(LM test):
## Test statistic: 1043.197 p-value: 0
## Rank-based Test:
## Test statistic: 849.5804 p-value: 0
## [2] "m = 20"
## Q(m) of squared series(LM test):
## Test statistic: 1705.592 p-value: 0
## Rank-based Test:
## Test statistic: 1565.588 p-value: 0

2. 标准GARCH模型建立

上述ARCH效应表明，条件方差是依赖于过去值。因此可以考虑GARCH模型对方差方程进行参数估计。

使用tseries包中的garch()函数进行拟合标准GARCH模型。
从结果上看，拟合出来的参数都显著，Box-Ljung test结果中的P值大于显著性，因此可以认为模型的残差无序列相关，说明该模型拟合效果较好。但实际上，其中Jarque Bera Test用于对回归残差的正态性进行检验，Shapiro - Wilk Normality Test也可以用于正态性检验，原假设都是是残差序列服从正态分布，检验结果表明，残差序列是不服从正态分布，因此可以对模型进行优化，考虑其他GARCH模型。

## Call:
## garch(x = r.data, order = c(1, 1))
##
## Model:
## GARCH(1,1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.91948 -0.53035 0.04107 0.57992 5.60582
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## a0 0.011014 0.002491 4.422 9.78e-06 ***
## a1 0.063407 0.004076 15.556 < 2e-16 ***
## b1 0.934953 0.003839 243.530 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Diagnostic Tests:
## Jarque Bera Test
##
## data: Residuals
## X-squared = 659.37, df = 2, p-value < 2.2e-16
##
##
## Box-Ljung test
##
## data: Squared.Residuals
## X-squared = 1.5224, df = 1, p-value = 0.2173

四、模型优化

（一）模型拟合

上述的模型考虑的是ARIMA(0,0,0)ARIMA(0,0,0)ARIMA(0,0,0)的标准GARCH(1,1)GARCH(1,1)GARCH(1,1)模型，即均值模型的参数均设置为0。而从均值模型的分析来看，可以拟合ARIMA(4,0,4)ARIMA(4,0,4)ARIMA(4,0,4)的均值模型与GARCH(1,1)GARCH(1,1)GARCH(1,1)的方差模型。而上述的正态性检验结果表明，残差的分布不适用标准正态分布，应该考虑其他类型的分布。

使用rugarchrugarchrugarch包对多个模型进行拟合。
对于均值模型，考虑不带截距项的ARIMA(4,0,4)ARIMA(4,0,4)ARIMA(4,0,4)。
对于方差模型，阶数设定1阶ARCH和1阶GARCH，考虑标准GARCH（sGARCHsGARCHsGARCH）、指数GARCH（eGARCHeGARCHeGARCH）、GJR−GARCHGJR-GARCHGJR−GARCH、渐近幂ARCH（APARCHAPARCHAPARCH）、门限GARCH（TGARCHTGARCHTGARCH）、非线性非对称GARCH（NAGARCHNAGARCHNAGARCH）六类模型。
对于残差分布类型，考虑标准正态分布（normnormnorm）、标准t分布（stdstdstd）、偏t分布（sstdsstdsstd）、广义误差分布（gedgedged）和Johnson’SU分布（jsujsujsu）五类分布。

（二）模型比较

从拟合的模型结果来看，非对称类型的指数GARCH模型在最大似然估计值达到最大，同时AIC和BIC都能达到最小，说明指数GARCH模型比较好。从拟合的模型残差分布来看，非正态分布的AIC和BIC都明显低于正态分布，说明残差是服从重尾类型的分布。

选取指数GARCH模型，对比不同分布的参数显著性、检验结果以及模型效果，这里分布考虑正态分布、t分布、广义误差分布与三种分布对应的偏态分布以及Johnson’SU分布，结果如下所示。

从模型拟合的系数显著性表可以看出，拟合的系数基本都显著，仅有少数参数不显著。从系数稳定性的个别检验和联合检验可以看出，在5%显著性水平下正态分布、偏正态分布、广义误差分布和偏广义误差分布都接收参数是稳定的原假设。从符号偏误检验的结果来看，指数GARCH模型正负残差的受到冲击的差异不明显，说明该非对称模型有效消除了杠杆效应。从调整皮尔逊拟合优度检验的结果来看，原假设是残差分布与理论分布没有差异，结果表明广义误差分布和偏广义误差分布是没有拒绝原假设，说明这两种分布与模型适配较好。

（三）模型选择

从模型拟合的信息准则表来看，偏广义误差分布的LLH达到最大，而BIC和HQ都达到最小，说明ARMA(4,4)−eGARCH(1,1)−SGEDARMA(4,4)-eGARCH(1,1)-SGEDARMA(4,4)−eGARCH(1,1)−SGED模型最优。拟合ARMA(4,4)−EGARCH(1,1)−SGEDARMA(4,4)-EGARCH(1,1)-SGEDARMA(4,4)−EGARCH(1,1)−SGED模型，理论模型和模型参数估计及显著性如下所示。
Xt=∑i=1pφiXt−i−∑j=1qϑjεt−j+εtzt=εtσt=Xtσtlnσt2=ω+αzt−1+βlnσt−12+γ(∣εt−1∣σt−1−E∣zt−1∣)其中，εt∼SGED(0,1,shape,skew)\begin{aligned} X_t &= \sum \limits _{i=1}^{p}\varphi_i X_{t-i} - \sum \limits _{j=1}^{q} \vartheta_j \varepsilon_{t-j} + \varepsilon_t \\ z_t &= \frac{\varepsilon_t}{\sigma_t} = \frac{X_t}{\sigma_t} \\ ln\sigma_t^2 &= \omega +\alpha z_{t-1} + \beta ln\sigma ^2_{t-1} + \gamma (\frac{|\varepsilon_{t-1}|}{\sigma _{t-1}} - E|z_{t-1}|) \\ 其中，& \varepsilon_t \sim SGED(0,1,shape, skew) \end{aligned} Xtztlnσt2其中，=i=1∑pφiXt−i−j=1∑qϑjεt−j+εt=σtεt=σtXt=ω+αzt−1+βlnσt−12+γ(σt−1∣εt−1∣−E∣zt−1∣)εt∼SGED(0,1,shape,skew)

五、结论

本文通过对沪深300指数的波动性分析发现，我国股票市场有两段时间出现较大的波动。第一次波动出现在2008年前后，这段期间为全球金融危机，明显可以看出收盘价时序图出现明显的陡峭的波峰，持续时间较长；第二次波动出现在2015年前后，该阶段是由于杠杆资金的加入和政策收紧，形成短暂的牛市和熊市，波动程度不亚于2008年金融危机，但持续时间比较短。
本文使用多个GARCH模型进行比对，发现非对称模型ARIMA−EGARCHARIMA-EGARCHARIMA−EGARCH模型与沪深300指数的对数收益率有较高的匹配度，同时偏广义误差分布与理论分布较为接近，说明沪深300指数的波动性呈现的是尖峰厚尾、非对称性特征，说明我国股票市场存在杠杆效应。

六、实现代码

# Created by: Enguanei
# Created on: 2021/5/7# 导入包
library(knitr)
library(pedquant)
library(zoo)
library(imputeTS)
library(tseries)
library(forecast)
library(MTS)
library(rugarch)
library(dplyr)# 爬取网易财经的数据
datt <- md_stock(symbol = '000300.ss',from = "2002-01-01",to = "2019-12-31", source = "163")# 查看数据
head(datt$`000300.ss`$close)
head(datt$`000300.ss`$date)# 保存数据
write.csv(datt, "D:\\Rproject\\hs300.csv")# 导入数据
raw_data <- read.csv("D:\\Rproject\\hs300.csv", header = TRUE)# 转成时间序列
my_date <- as.Date(raw_data$X000300.ss.date)
my_data <- raw_data$X000300.ss.close
data.ts <- zoo(my_data, my_date)# 查看是否存在缺失值
sum(is.na(my_data))# 转成时间序列（不包括日期）和对数差分化处理
r.data <- diff(log(ts(my_data))) * 100# 画出时序图
par(mfrow = c(2, 1))
plot(data.ts, xlab = 'Year', ylab = 'Close.Price',main = 'Close of CSI300')
plot(r.data, ylab = 'Log.Return',main = 'Log.Return of CSI300')# 平稳性检验
show(adf.test(r.data))# 拟合arma均值模型
md1 <- auto.arima(r.data)
md1# 系数显著性检验
t <- md1$coef / sqrt(diag(md1$var.coef))
p_1 <- 2 * (1 - pnorm(abs(t)))
kable(rbind(p = p_1), caption = "Result of Coef. Test")# 拟合疏系数arma
mean_md_1 <- Arima(r.data,order = c(4, 0, 4),fixed = c(NA, 0, 0, NA, NA, 0, 0, NA, 0))# 疏系数模型的系数显著性检验
t <- mean_md_1$coef / sqrt(diag(mean_md_1$var.coef))
p_2 <- 2 * (1 - pnorm(abs(t)))
kable(rbind(p = p_2[p_2 != 1]), caption = "Result of Coef. Test")# ARCH效应检验（原假设是没有ARCH效应）
print(paste0('m = 10'))
archTest(mean_md_1$res, lag = 10)
print(paste0('m = 20'))
archTest(mean_md_1$res, lag = 20)# 使用tseries拟合sGARCH（只能拟合方差模型，没有考虑均值模型）
out <- garch(r.data, order = c(1, 1))
summary(out)# 模型优化
# 模型选择
md_name <- c("sGARCH", "eGARCH", "gjrGARCH","apARCH", "TGARCH", "NAGARCH")
md_dist <- c("norm", "std", "sstd", "ged", "jsu")# 设定函数
model_train <- function(data_, m_name, m_dist) {res <- NULLfor (i in m_name) {sub_model <- NULLmain_model <- iif (i == "TGARCH" | i == "NAGARCH") {sub_model <- imain_model <- "fGARCH"}for (j in m_dist) {model_name <- paste0('ARMA(4,4)', '-',i, '-', j,sep = "")if (main_model == "fGARCH") {model_name <- paste0('ARMA(4,4)', '-',sub_model, '-', j,sep = "")}print(model_name)# 模型设定mean.spec <- list(armaOrder = c(4, 4), include.mean = F,archm = F, archpow = 1, arfima = F,external.regressors = NULL)var.spec <- list(model = main_model, garchOrder = c(1, 1),submodel = sub_model,external.regressors = NULL,variance.targeting = F)dist.spec <- jmy_spec <- ugarchspec(mean.model = mean.spec,variance.model = var.spec,distribution.model = dist.spec)# 模型拟合my_fit <- ugarchfit(data = data_, spec = my_spec)res <- rbind(res, c(list(ModelName = model_name),list(LogL = likelihood(my_fit)),list(AIC = infocriteria(my_fit)[1]),list(BIC = infocriteria(my_fit)[2]),list(RMSE = sqrt(mean(residuals(my_fit)^2))),list(md = my_fit)))}}return(res)
}# 模型拟合
my_model <- model_train(r.data, md_name, md_dist)# 保存模型结果
write.table(my_model[, 1:5],"D:\\Rproject\\result.txt",sep = " ")# 读取模型结果
res_table <- read.table("D:\\Rproject\\result.txt",header = TRUE,sep = " ")
kable(res_table, caption = "Result of GARCH Model", align = 'l')# 选择模型再拟合
md_name <- c("eGARCH")
md_dist <- c("norm", "snorm", "std", "sstd", "ged", "sged", "jsu")
new_model <- model_train(r.data, md_name, md_dist)# 信息准则表
my_info_cri <- NULL
for (i in 1:dim(new_model)[1]) {new_aic <- round(infocriteria(new_model[i, 6]$md)[1], digits = 4)new_bic <- round(infocriteria(new_model[i, 6]$md)[2], digits = 4)new_sib <- round(infocriteria(new_model[i, 6]$md)[3], digits = 4)new_hq <- round(infocriteria(new_model[i, 6]$md)[4], digits = 4)new_llk <- round(likelihood(new_model[i, 6]$md), digits = 4)my_info_cri <- rbind(my_info_cri, c(new_model[i, 1],Akaike = new_aic,Bayes = new_bic,Shibata = new_sib,HQ = new_hq,LLH = new_llk))
}# 模型系数显著性表
my_cof_table <- NULL
my_cof_names <- NULL
for (i in 1:dim(new_model)[1]) {my_cof_names <- rbind(my_cof_names, new_model[i, 1]$ModelName)my_df <- data.frame(t(round(new_model[i, 6]$md@fit$robust.matcoef[, 4],digits = 4)))if (i == 1) {my_cof_table <- data.frame(my_df)}else {my_cof_table <- full_join(my_cof_table, my_df)}
}
row.names(my_cof_table) <- my_cof_names[, 1]
my_cof_table <- t(my_cof_table)# 参数稳定性个别检验表
my_nyblom_table <- NULL
my_nyblom_name <- NULL
for (i in 1:dim(new_model)[1]) {my_nyblom_name <- rbind(my_nyblom_name,new_model[i, 1]$ModelName)my_df <- data.frame(t(round(nyblom(new_model[i, 6]$md)$IndividualStat,digits = 4)))if (i == 1) {my_nyblom_table <- data.frame(my_df)}else {my_nyblom_table <- full_join(my_nyblom_table, my_df)}
}
row.names(my_nyblom_table) <- my_nyblom_name[, 1]
my_IC <- NULL
for (i in 1:dim(new_model)[1]) {my_IC <- cbind(my_IC, nyblom(new_model[1, 6]$md)$IndividualCritical)
}
my_nyblom_table <- rbind(t(my_nyblom_table), my_IC)# 参数稳定性联合检验表
my_nyblom_table2 <- NULL
my_nyblom_name2 <- NULL
for (i in 1:dim(new_model)[1]) {my_nyblom_name2 <- rbind(my_nyblom_name2,new_model[i, 1]$ModelName)df1 <- data.frame(JoinStat = round(nyblom(new_model[i,6]$md)$JointStat,digits = 4))df2 <- data.frame(mm = round(nyblom(new_model[i, 6]$md)$JointCritical,digits = 4))temp_table <- cbind(df1, t(df2))if (i == 1) {my_nyblom_table2 <- temp_table}else {my_nyblom_table2 <- full_join(my_nyblom_table2, temp_table)}
}
row.names(my_nyblom_table2) <- my_nyblom_name2# 模型分布拟合度p值检验表
my_gof_table <- NULL
my_gof_name <- NULL
group_name <- c(20, 30, 40, 50)
for (i in 1:dim(new_model)[1]) {my_gof_name <- rbind(my_gof_name,new_model[i, 1]$ModelName)df1 <- round(gof(new_model[i, 6]$md,groups = group_name)[1:4, 3],digits = 4)if (i == 1) {my_gof_table <- df1}else {my_gof_table <- rbind(my_gof_table, df1)}
}
row.names(my_gof_table) <- my_gof_name
colnames(my_gof_table) <- c(20, 30, 40, 50)# 符号偏误显著性检验表
my_sign_table <- NULL
my_sign_name <- NULL
row_sign_name <- row.names(signbias(new_model[1, 6]$md))
for (i in 1:dim(new_model)[1]) {my_sign_name <- rbind(my_sign_name,new_model[i, 1]$ModelName)df1 <- round(signbias(new_model[i, 6]$md)$prob, digits = 4)if (i == 1) {my_sign_table <- df1}else {my_sign_table <- rbind(my_sign_table, df1)}
}
row.names(my_sign_table) <- my_sign_name
colnames(my_sign_table) <- row_sign_name# 展示表格
kable(my_cof_table,caption = "P-value Table of Coefficients")
kable(my_nyblom_table,caption = "P-value Table of Nyblom Individual Stability Test")
kable(my_nyblom_table2,caption = "P-value Table of Nyblom Joint Stability Test")
kable(my_sign_table,caption = "P-value Table of Sign Bias Test")
kable(my_gof_table,caption = "P-value Table of Goodness-of-fit")
kable(my_info_cri,caption = "Table of Information Criteria")# 最终模型系数表
kable(new_model[6, 6]$md@fit$robust.matcoef,caption = "Table of Final Model Coef.")

七、参考资料

[1] Ruey S. Tsay, 李洪成, 尚秀芬,等. 金融数据分析导论[M]. 机械工业出版社, 2013.
[2] 何宗武, 马卫锋. 经济与金融计量方法：原理、应用案例及R语言实现[M]. 机械工业出版社, 2019
[3] 张东旭. 基于ARMA-GARCH模型族的上证指数收益率波动的实证分析[D]. 清华大学.

【R语言】GARCH模型的应用相关推荐

R语言GARCH族模型：正态分布、t、GED分布EGARCH、TGARCH的VaR分析股票指数
全文链接:http://tecdat.cn/?p=31023 如何构建合适的模型以恰当的方法对风险进行测量是当前金融研究领域的一个热门话题(点击文末"阅读原文"获取完整代码数据). ...
R语言GARCH建模常用软件包比较、拟合标准普尔SP 500指数波动率时间序列和预测可视化...
原文链接:http://tecdat.cn/?p=24441 我们研究波动聚集,以及使用单变量 GARCH(1,1) 模型对其进行建模. 波动聚集波动聚集--存在相对平稳时期和高波动时期的现象--是 ...
R语言分类模型：逻辑回归模型LR、决策树DT、推理决策树CDT、随机森林RF、支持向量机SVM、Rattle可视化界面数据挖掘、分类模型评估指标（准确度、敏感度、特异度、PPV、NPV）
R语言分类模型:逻辑回归模型LR.决策树DT.推理决策树CDT.随机森林RF.支持向量机SVM.Rattle可视化界面数据挖掘.分类模型评估指标(准确度.敏感度.特异度.PPV.NPV) 目录
R语言回归模型构建、回归模型基本假设（正态性、线性、独立性、方差齐性）、回归模型诊断、car包诊断回归模型、特殊观察样本分析、数据变换、模型比较、特征筛选、交叉验证、预测变量相对重要度
R语言回归模型构建.回归模型基本假设(正态性.线性.独立性.方差齐性).回归模型诊断.car包诊断回归模型.特殊观察样本分析.数据变换.模型比较.特征筛选.交叉验证.预测变量相对重要度目录
R语言KNN模型数据分类实战
R语言KNN模型数据分类实战目录 R语言KNN模型数据分类实战 #KNN基本原理 #案例分析 #数据集说明
R语言xgboost模型构建:基于prima糖尿病数据集
R语言xgboost模型构建:基于prima糖尿病数据集目录 R语言xgboost模型构建:基于prima糖尿病数据集
R语言回归模型协方差分析（Analysis of Covariance）
R语言回归模型协方差分析(Analysis of Covariance) 目录 R语言回归模型协方差分析(Analysis of Covariance) 输入数据 ANCOVA分析
R语言glm模型预测（predict）过程及Error in eval(predvars, data, env) 错误原因
R语言glm模型预测(predict)过程及Error in eval(predvars, data, env) 错误原因目录 R语言glm模型预测(predict)过程及Error in eval ...
R语言︱机器学习模型评价指标+（转）模型出错的四大原因及如何纠错
笔者寄语:机器学习中交叉验证的方式是主要的模型评价方法,交叉验证中用到了哪些指标呢? 交叉验证将数据分为训练数据集.测试数据集,然后通过训练数据集进行训练,通过测试数据集进行测试,验证集进行验证. 模 ...
R语言结构方程模型（SEM）在生态学领域中的应用
前言:结构方程模型(Sructural Equation Model)是一种建立.估计和检验研究系统中多变量间因果关系的模型方法,它可以替代多元回归.因子分析.协方差分析等方法,利用图形化模型方式清晰 ...

【R语言】GARCH模型的应用