本文是時間序列分析課程的作業，基于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)
【R語言】GARCH模型的應用一、資料來源二、資料分析三、模型建立四、模型優化五、結論六、實作代碼七、參考資料 - 采用疏系數的 A R M A ARMA ARMA模型後的參數顯著性檢驗結果來看，非零參數結果都顯著。（注：這一步存有疑問，後續依舊是采用 A R M A ( 4 , 4 ) ARMA(4,4) ARMA(4,4)的非疏系數模型進行拟合）
【R語言】GARCH模型的應用一、資料來源二、資料分析三、模型建立四、模型優化五、結論六、實作代碼七、參考資料
（二）方差模型建立
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
四、模型優化
（一）模型拟合
上述的模型考慮的是 A R I M A ( 0 , 0 , 0 ) ARIMA(0,0,0) ARIMA(0,0,0)的标準 G A R C H ( 1 , 1 ) GARCH(1,1) GARCH(1,1)模型，即均值模型的參數均設定為0。而從均值模型的分析來看，可以拟合 A R I M A ( 4 , 0 , 4 ) ARIMA(4,0,4) ARIMA(4,0,4)的均值模型與 G A R C H ( 1 , 1 ) GARCH(1,1) GARCH(1,1)的方差模型。而上述的正态性檢驗結果表明，殘差的分布不适用标準正态分布，應該考慮其他類型的分布。
- 使用 r u g a r c h rugarch rugarch包對多個模型進行拟合。
- 對于均值模型，考慮不帶截距項的 A R I M A ( 4 , 0 , 4 ) ARIMA(4,0,4) ARIMA(4,0,4)。
- 對于方差模型，階數設定1階ARCH和1階GARCH，考慮标準GARCH（ s G A R C H sGARCH sGARCH）、指數GARCH（ e G A R C H eGARCH eGARCH）、 G J R − G A R C H GJR-GARCH GJR−GARCH、漸近幂ARCH（ A P A R C H APARCH APARCH）、門限GARCH（ T G A R C H TGARCH TGARCH）、非線性非對稱GARCH（ N A G A R C H NAGARCH NAGARCH）六類模型。
- 對于殘差分布類型，考慮标準正态分布（ n o r m norm norm）、标準t分布（ s t d std std）、偏t分布（ s s t d sstd sstd）、廣義誤差分布（ g e d ged ged）和Johnson’SU分布（ j s u jsu jsu）五類分布。
【R語言】GARCH模型的應用一、資料來源二、資料分析三、模型建立四、模型優化五、結論六、實作代碼七、參考資料
（二）模型比較
從拟合的模型結果來看，非對稱類型的指數GARCH模型在最大似然估計值達到最大，同時AIC和BIC都能達到最小，說明指數GARCH模型比較好。從拟合的模型殘差分布來看，非正态分布的AIC和BIC都明顯低于正态分布，說明殘差是服從重尾類型的分布。
- 選取指數GARCH模型，對比不同分布的參數顯著性、檢驗結果以及模型效果，這裡分布考慮正态分布、t分布、廣義誤差分布與三種分布對應的偏态分布以及Johnson’SU分布，結果如下所示。
【R語言】GARCH模型的應用一、資料來源二、資料分析三、模型建立四、模型優化五、結論六、實作代碼七、參考資料 【R語言】GARCH模型的應用一、資料來源二、資料分析三、模型建立四、模型優化五、結論六、實作代碼七、參考資料 【R語言】GARCH模型的應用一、資料來源二、資料分析三、模型建立四、模型優化五、結論六、實作代碼七、參考資料 【R語言】GARCH模型的應用一、資料來源二、資料分析三、模型建立四、模型優化五、結論六、實作代碼七、參考資料 【R語言】GARCH模型的應用一、資料來源二、資料分析三、模型建立四、模型優化五、結論六、實作代碼七、參考資料 【R語言】GARCH模型的應用一、資料來源二、資料分析三、模型建立四、模型優化五、結論六、實作代碼七、參考資料
從模型拟合的系數顯著性表可以看出，拟合的系數基本都顯著，僅有少數參數不顯著。從系數穩定性的個别檢驗和聯合檢驗可以看出，在5%顯著性水準下正态分布、偏正态分布、廣義誤差分布和偏廣義誤差分布都接收參數是穩定的原假設。從符号偏誤檢驗的結果來看，指數GARCH模型正負殘差的受到沖擊的差異不明顯，說明該非對稱模型有效消除了杠杆效應。從調整皮爾遜拟合優度檢驗的結果來看，原假設是殘差分布與理論分布沒有差異，結果表明廣義誤差分布和偏廣義誤差分布是沒有拒絕原假設，說明這兩種分布與模型适配較好。
（三）模型選擇
從模型拟合的資訊準則表來看，偏廣義誤差分布的LLH達到最大，而BIC和HQ都達到最小，說明 A R M A ( 4 , 4 ) − e G A R C H ( 1 , 1 ) − S G E D ARMA(4,4)-eGARCH(1,1)-SGED ARMA(4,4)−eGARCH(1,1)−SGED模型最優。拟合 A R M A ( 4 , 4 ) − E G A R C H ( 1 , 1 ) − S G E D ARMA(4,4)-EGARCH(1,1)-SGED ARMA(4,4)−EGARCH(1,1)−SGED模型，理論模型和模型參數估計及顯著性如下所示。
X t = ∑ i = 1 p φ i X t − i − ∑ j = 1 q ϑ j ε t − j + ε t z t = ε t σ t = X t σ t l n σ t 2 = ω + α z t − 1 + β l n σ t − 1 2 + γ ( ∣ ε t − 1 ∣ σ t − 1 − E ∣ z t − 1 ∣ ) 其 中 ， ε t ∼ S G E D ( 0 , 1 , s h a p e , s k e w ) \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模型進行比對，發現非對稱模型 A R I M A − E G A R C H ARIMA-EGARCH 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][email protected]$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][email protected]$robust.matcoef,
caption = "Table of Final Model Coef.")
七、參考資料
[1] Ruey S. Tsay, 李洪成, 尚秀芬,等. 金融資料分析導論[M]. 機械工業出版社, 2013.
[2] 何宗武, 馬衛鋒. 經濟與金融計量方法：原理、應用案例及R語言實作[M]. 機械工業出版社, 2019
[3] 張東旭. 基于ARMA-GARCH模型族的上證指數收益率波動的實證分析[D]. 清華大學.