garch模型建模 python garch模型步骤

  本文是时间序列分析课程的作业，基于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)$进行参数显著性检验。检验结果发现部分参数不显著，采用建立疏系数的均值模型，将不显著的参数强制为0。
• 采用疏系数的$ARMA$模型后的参数显著性检验结果来看，非零参数结果都显著。（注：这一步存有疑问，后续依旧是采用$ARMA(4,4)$的非疏系数模型进行拟合）

（二）方差模型建立

1. ARCH效应检验

  上述建立模型后，对残差进行ARCH效应检验。Ljung-Box统计量$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)$的标准$GARCH(1,1)$模型，即均值模型的参数均设置为0。而从均值模型的分析来看，可以拟合$ARIMA(4,0,4)$的均值模型与$GARCH(1,1)$的方差模型。而上述的正态性检验结果表明，残差的分布不适用标准正态分布，应该考虑其他类型的分布。

• 使用$rugarch$包对多个模型进行拟合。
• 对于均值模型，考虑不带截距项的$ARIMA(4,0,4)$。
• 对于方差模型，阶数设定1阶ARCH和1阶GARCH，考虑标准GARCH（$sGARCH$）、指数GARCH（$eGARCH$）、$GJR-GARCH$、渐近幂ARCH（$APARCH$）、门限GARCH（$TGARCH$）、非线性非对称GARCH（$NAGARCH$）六类模型。
• 对于残差分布类型，考虑标准正态分布（$norm$）、标准t分布（$std$）、偏t分布（$sstd$）、广义误差分布（$ged$）和Johnson’SU分布（$jsu$）五类分布。

（二）模型比较

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

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

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

（三）模型选择

$ARMA(4,4)-eGARCH(1,1)-SGED$模型最优。拟合$ARMA(4,4)-EGARCH(1,1)-SGED$模型，理论模型和模型参数估计及显著性如下所示。

$\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}$

五、结论

  本文通过对沪深300指数的波动性分析发现，我国股票市场有两段时间出现较大的波动。第一次波动出现在2008年前后，这段期间为全球金融危机，明显可以看出收盘价时序图出现明显的陡峭的波峰，持续时间较长；第二次波动出现在2015年前后，该阶段是由于杠杆资金的加入和政策收紧，形成短暂的牛市和熊市，波动程度不亚于2008年金融危机，但持续时间比较短。
  本文使用多个GARCH模型进行比对，发现非对称模型$ARIMA-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 <- NULL
  for (i in m_name) {
    sub_model <- NULL
    main_model <- i
    if (i == "TGARCH" | i == "NAGARCH") {
      sub_model <- i
      main_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 <- j
      my_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]. 清华大学.