R語言用多元ARMA,GARCH ,EWMA, ETS,隨機波動率SV模型對金融時間序列數(shù)據(jù)建模
原文鏈接:http://tecdat.cn/?p=20015?
?
本文將說明單變量和多變量金融時間序列的不同模型,特別是條件均值和條件協(xié)方差矩陣、波動率的模型。
均值模型
本節(jié)探討條件均值模型。
iid模型
我們從簡單的iid模型開始。iid模型假定對數(shù)收益率xt為N維高斯時間序列:

均值和協(xié)方差矩陣的樣本估計量分別是樣本均值

和樣本協(xié)方差矩陣

我們從生成數(shù)據(jù)開始,熟悉該過程并確保估計過程給出正確的結果(即完整性檢查)。然后使用真實的市場數(shù)據(jù)并擬合不同的模型。
讓我們生成合成iid數(shù)據(jù)并估算均值和協(xié)方差矩陣:
# 生成綜合收益數(shù)據(jù)
X <- rmvnorm(n = T, mean = mu, sigma = Sigma)
# 樣本估計(樣本均值和樣本協(xié)方差矩陣)
mu_sm <- colMeans(X)
Sigma_scm <- cov(X)
# 誤差
norm(mu_sm ? ? - mu, "2")
#> [1] 2.44
norm(Sigma_scm - Sigma, "F")
#> [1] 70.79
現(xiàn)在,讓我們針對不同數(shù)量的觀測值T再做一次:
# 首先生成所有數(shù)據(jù)
X <- rmvnorm(n = T_max, mean = mu, sigma = Sigma)
# 現(xiàn)在遍歷樣本的子集
for (T_ in T_sweep) {
?# 樣本估算
?mu_sm <- colMeans(X_)
?Sigma_scm <- cov(X_)
?# 計算誤差
?error_mu_vs_T ? ?<- c(error_mu_vs_T, ? ?norm(mu_sm ? ? - mu, "2"))
?error_Sigma_vs_T <- c(error_Sigma_vs_T, norm(Sigma_scm - Sigma, "F"))
# 繪圖
plot(T_sweep, error_mu_vs_T,
? ? main = "mu估計誤差",

plot(T_sweep, error_Sigma_vs_T
? ? main = "Sigma估計中的誤差", ylab = "誤差"

單變量ARMA模型
對數(shù)收益率xt上的ARMA(p,q)模型是

其中wt是均值為零且方差為σ2的白噪聲序列。模型的參數(shù)是系數(shù)?i,θi和噪聲方差σ2。
請注意,ARIMA(p,d,q)模型是時間差分為d階的ARMA(p,q)模型。因此,如果我們用xt代替對數(shù)價格,那么先前的對數(shù)收益模型實際上就是ARIMA(p,1,q)模型,因為一旦對數(shù)價格差分,我們就獲得對數(shù)收益。
rugarch生成數(shù)據(jù)?
我們將使用rugarch包??生成單變量ARMA數(shù)據(jù),估計參數(shù)并進行預測。
首先,我們需要定義模型:
# 指定具有給定系數(shù)和參數(shù)的AR(1)模型
#>
#> *----------------------------------*
#> * ? ? ? ARFIMA Model Spec ? ? ? ? ?*
#> *----------------------------------*
#> Conditional Mean Dynamics
#> ------------------------------------
#> Mean Model ? ? ? ? ? : ARFIMA(1,0,0)
#> Include Mean ? ? : TRUE
#>
#> Conditional Distribution
#> ------------------------------------
#> Distribution : ?norm
#> Includes Skew ? ?: ?FALSE
#> Includes Shape ? : ?FALSE
#> Includes Lambda ?: ?FALSE
#> ? ? ? ? ?Level Fixed Include Estimate LB UB
#> mu ? ? ? ?0.01 ? ? 1 ? ? ? 1 ? ? ? ?0 NA NA
#> ar1 ? ? ?-0.90 ? ? 1 ? ? ? 1 ? ? ? ?0 NA NA
#> ma ? ? ? ?0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> arfima ? ?0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> archm ? ? 0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> mxreg ? ? 0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> sigma ? ? 0.20 ? ? 1 ? ? ? 1 ? ? ? ?0 NA NA
#> alpha ? ? 0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> beta ? ? ?0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> gamma ? ? 0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> eta1 ? ? ?0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> eta2 ? ? ?0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> delta ? ? 0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> lambda ? ?0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> vxreg ? ? 0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> skew ? ? ?0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> shape ? ? 0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> ghlambda ?0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> xi ? ? ? ?0.00 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
fixed.pars
#> $mu
#> [1] 0.01
#>
#> $ar1
#> [1] -0.9
#>
#> $sigma
#> [1] 0.2
true_params
#> ? ?mu ? ar1 sigma
#> ?0.01 -0.90 ?0.20
然后,我們可以生成時間序列:
# 模擬一條路徑
apath(spec, n.sim = T)
# 轉(zhuǎn)換為xts并繪圖
plot(synth_log_returns, main = "ARMA模型的對數(shù)收益率"
plot(synth_log_prices, main = "ARMA模型的對數(shù)價格"


ARMA模型
現(xiàn)在,我們可以估計參數(shù)(我們已經(jīng)知道):
# 指定AR(1)模型
arfimaspec(mean.model = list(armaOrder = c(1,0), include.mean = TRUE))
# 估計模型
#> ? ? ? ? ? mu ? ? ? ? ?ar1 ? ? ? ?sigma
#> ? ? ? 0.0083 ? ? ?-0.8887 ? ? ? 0.1987
#> ? ?mu ? ar1 sigma
#> ?0.01 -0.90 ?0.20
我們還可以研究樣本數(shù)量T對參數(shù)估計誤差的影響:
# 循環(huán)
for (T_ in T_sweep) {
?estim_coeffs_vs_T <- rbind(estim_coeffs_vs_T, coef(arma_fit))
?error_coeffs_vs_T <- rbind(error_coeffs_vs_T, abs(coef(arma_fit) - true_params)/true_params)
# 繪圖
matplot(T_sweep, estim_coeffs_vs_T,
? ? ? ?main = "估計的ARMA系數(shù)", xlab = "T", ylab = "值",

matplot(T_sweep, 100*error_coeffs_vs_T,
? ? ? ?main = "估計ARMA系數(shù)的相對誤差", xlab = "T", ylab = "誤差 (%)",

首先,真正的μ幾乎為零,因此相對誤差可能顯得不穩(wěn)定。在T = 800個樣本之后,其他系數(shù)得到了很好的估計。
ARMA預測
為了進行健全性檢查,我們現(xiàn)在將比較兩個程序包?Forecast?和?rugarch的結果:
# 指定具有給定系數(shù)和參數(shù)的AR(1)模型
spec(mean.model = list(armaOrder = c(1,0), include.mean = TRUE),
? ? ? ? ? ? ? ? ? ? ? ? ? ? fixed.pars = list(mu = 0.005, ar1 = -0.9, sigma = 0.1))
# 生成長度為1000的序列
arfima(arma_fixed_spec, n.sim = 1000)@path$seriesSim
# 使用 rugarch包指定和擬合模型
spec(mean.model = list(armaOrder = c(1,0), include.mean = TRUE))
# 使用包“ forecast”擬合模型
#> ARIMA(1,0,0) with non-zero mean
#>
#> Coefficients:
#> ? ? ? ? ? ar1 ? ?mean
#> ? ? ? -0.8982 ?0.0036
#> s.e. ? 0.0139 ?0.0017
#>
#> sigma^2 estimated as 0.01004: ?log likelihood=881.6
#> AIC=-1757.2 ? AICc=-1757.17 ? BIC=-1742.47
# 比較模型系數(shù)
#> ? ? ? ? ?ar1 ? ?intercept ? ? ? ?sigma
#> -0.898181148 ?0.003574781 ?0.100222964
#> ? ? ? ? ? mu ? ? ? ? ?ar1 ? ? ? ?sigma
#> ?0.003605805 -0.898750138 ?0.100199956
確實,這兩個軟件包給出了相同的結果。
ARMA模型選擇?
在先前的實驗中,我們假設我們知道ARMA模型的階數(shù),即p = 1和q = 0。實際上,階數(shù)是未知的,因此必須嘗試不同的階數(shù)組合。階數(shù)越高,擬合越好,但這將不可避免地導致過度擬合。已經(jīng)開發(fā)出許多方法來懲罰復雜性的增加以避免過度擬合,例如AIC,BIC,SIC,HQIC等。
# 嘗試不同的組合
# 查看排名
#> ? ?AR MA Mean ARFIMA ? ? ? ? BIC converged
#> 1 ? 1 ?0 ? ?1 ? ? ?0 -0.38249098 ? ? ? ? 1
#> 2 ? 1 ?1 ? ?1 ? ? ?0 -0.37883157 ? ? ? ? 1
#> 3 ? 2 ?0 ? ?1 ? ? ?0 -0.37736340 ? ? ? ? 1
#> 4 ? 1 ?2 ? ?1 ? ? ?0 -0.37503980 ? ? ? ? 1
#> 5 ? 2 ?1 ? ?1 ? ? ?0 -0.37459177 ? ? ? ? 1
#> 6 ? 3 ?0 ? ?1 ? ? ?0 -0.37164609 ? ? ? ? 1
#> 7 ? 1 ?3 ? ?1 ? ? ?0 -0.37143480 ? ? ? ? 1
#> 8 ? 2 ?2 ? ?1 ? ? ?0 -0.37107841 ? ? ? ? 1
#> 9 ? 3 ?1 ? ?1 ? ? ?0 -0.36795491 ? ? ? ? 1
#> 10 ?2 ?3 ? ?1 ? ? ?0 -0.36732669 ? ? ? ? 1
#> 11 ?3 ?2 ? ?1 ? ? ?0 -0.36379209 ? ? ? ? 1
#> 12 ?3 ?3 ? ?1 ? ? ?0 -0.36058264 ? ? ? ? 1
#> 13 ?0 ?3 ? ?1 ? ? ?0 -0.11875575 ? ? ? ? 1
#> 14 ?0 ?2 ? ?1 ? ? ?0 ?0.02957266 ? ? ? ? 1
#> 15 ?0 ?1 ? ?1 ? ? ?0 ?0.39326050 ? ? ? ? 1
#> 16 ?0 ?0 ? ?1 ? ? ?0 ?1.17294875 ? ? ? ? 1
#選最好的
armaOrder
#> AR MA
#> ?1 ?0
在這種情況下,由于觀察次數(shù)T = 1000足夠大,因此階數(shù)被正確地檢測到。相反,如果嘗試使用T = 200,則檢測到的階數(shù)為p = 1,q = 3。
ARMA預測?
一旦估計了ARMA模型參數(shù)?i? ^ i和θ^j,就可以使用該模型預測未來的值。例如,根據(jù)過去的信息對xt的預測是

并且預測誤差將為xt-x ^ t = wt(假設參數(shù)已被估計),其方差為σ2。軟件包?rugarch?使對樣本外數(shù)據(jù)的預測變得簡單:
# 估計模型(不包括樣本外)
coef(arma_fit)
#> ? ? ? ? ? mu ? ? ? ? ?ar1 ? ? ? ?sigma
#> ?0.007212069 -0.898745183 ?0.200400119
# 整個樣本外的預測對數(shù)收益
forecast_log_returns <- xts(arma_fore@forecast$seriesFor[1, ], dates_out_of_sample)
# 恢復對數(shù)價格
prev_log_price <- head(tail(synth_log_prices, out_of_sample+1), out_of_sample)
# 對數(shù)收益圖
plot(cbind("fitted" ? = fitted(arma_fit),
# 對數(shù)價格圖
plot(cbind("forecast" = forecast_log_prices,
? ? main = "對數(shù)價格預測", legend.loc = "topleft")


多元VARMA模型
對數(shù)收益率xt上的VARMA(p,q)模型是

其中wt是具有零均值和協(xié)方差矩陣Σw的白噪聲序列。該模型的參數(shù)是矢量/矩陣系數(shù)?0,Φi,Θj和噪聲協(xié)方差矩陣Σw。
?
比較
讓我們首先加載S&P500:
# 加載標普500數(shù)據(jù)
head(SP500_index_prices)
#> ? ? ? ? ? ? ?SP500
#> 2012-01-03 1277.06
#> 2012-01-04 1277.30
#> 2012-01-05 1281.06
#> 2012-01-06 1277.81
#> 2012-01-09 1280.70
#> 2012-01-10 1292.08
# 準備訓練和測試數(shù)據(jù)
logreturns_trn <- logreturns[1:T_trn]
logreturns_tst <- logreturns[-c(1:T_trn)]
# 繪圖
{ plot(logreturns,
?addEventLines(xts("訓練"

現(xiàn)在,我們使用訓練數(shù)據(jù)(即,對于t = 1,…,Ttrnt = 1,…,Ttrn)來擬合不同的模型(請注意,通過指示排除了樣本外數(shù)據(jù)?out.sample = T_tst
)。特別是,我們將考慮iid模型,AR模型,ARMA模型以及一些ARCH和GARCH模型(稍后將對方差建模進行更詳細的研究)。
# 擬合i.i.d.模型
coef(iid_fit)
#> ? ? ? ? ? mu ? ? ? ?sigma
#> 0.0005712982 0.0073516993
mean(logreturns_trn)
#> [1] 0.0005681388
sd(logreturns_trn)
#> [1] 0.007360208
# 擬合AR(1)模型
coef(ar_fit)
#> ? ? ? ? ? ?mu ? ? ? ? ? ar1 ? ? ? ? sigma
#> ?0.0005678014 -0.0220185181 ?0.0073532716
# 擬合ARMA(2,2)模型
coef(arma_fit)
#> ? ? ? ? ? ?mu ? ? ? ? ? ar1 ? ? ? ? ? ar2 ? ? ? ? ? ma1 ? ? ? ? ? ma2 ? ? ? ? sigma
#> ?0.0007223304 ?0.0268612636 ?0.9095552008 -0.0832923604 -0.9328475211 ?0.0072573570
# 擬合ARMA(1,1)+ ARCH(1)模型
coef(arch_fit)
#> ? ? ? ? ? ?mu ? ? ? ? ? ar1 ? ? ? ? ? ma1 ? ? ? ? omega ? ? ? ?alpha1
#> ?6.321441e-04 ?8.720929e-02 -9.391019e-02 ?4.898885e-05 ?9.986975e-02
#擬合ARMA(0,0)+ARCH(10)模型
coef(long_arch_fit)
#> ? ? ? ? ? mu ? ? ? ?omega ? ? ? alpha1 ? ? ? alpha2 ? ? ? alpha3 ? ? ? alpha4 ? ? ? alpha5
#> 7.490786e-04 2.452099e-05 6.888561e-02 7.207551e-02 1.419938e-01 1.909541e-02 3.082806e-02
#> ? ? ? alpha6 ? ? ? alpha7 ? ? ? alpha8 ? ? ? alpha9 ? ? ?alpha10
#> 4.026539e-02 3.050040e-07 9.260183e-02 1.150128e-01 1.068426e-06
# 擬合ARMA(1,1)+GARCH(1,1)模型
coef(garch_fit)
#> ? ? ? ? ? ?mu ? ? ? ? ? ar1 ? ? ? ? ? ma1 ? ? ? ? omega ? ? ? ?alpha1 ? ? ? ? beta1
#> ?6.660346e-04 ?9.664597e-01 -1.000000e+00 ?7.066506e-06 ?1.257786e-01 ?7.470725e-01
我們使用不同的模型來預測對數(shù)收益率:
# 準備預測樣本外周期的對數(shù)收益
# i.i.d.模型預測
forecast(iid_fit, n.ahead = 1, n.roll = T_tst - 1)
? ? ? ? ? ? ? ? ? ? ? ? ?dates_out_of_sample)
# AR(1)模型進行預測
forecast(ar_fit, n.ahead = 1, n.roll = T_tst - 1)
? ? ? ? ? ? ? ? ? ? ? ? dates_out_of_sample)
# ARMA(2,2)模型進行預測
forecast(arma_fit, n.ahead = 1, n.roll = T_tst - 1)
? ? ? ? ? ? ? ? ? ? ? ? ? dates_out_of_sample)
# 使用ARMA(1,1)+ ARCH(1)模型進行預測
forecast(arch_fit, n.ahead = 1, n.roll = T_tst - 1)
? ? ? ? ? ? ? ? ? ? ? ? ? dates_out_of_sample)
# ARMA(0,0)+ARCH(10)模型預測
forecast(long_arch_fit, n.ahead = 1, n.roll = T_tst - 1)
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?dates_out_of_sample)
# ARMA(1,1)+GARCH(1,1)模型預測
forecast(garch_fit, n.ahead = 1, n.roll = T_tst - 1)
? ? ? ? ? ? ? ? ? ? ? ? ? ?dates_out_of_sample)
我們可以計算不同模型的預測誤差(樣本內(nèi)和樣本外):
print(error_var)
#> ? ? ? ? ? ? ? ? ? ? ? ? ? in-sample out-of-sample
#> iid ? ? ? ? ? ? ? ? ? ?5.417266e-05 ?8.975710e-05
#> AR(1) ? ? ? ? ? ? ? ? ?5.414645e-05 ?9.006139e-05
#> ARMA(2,2) ? ? ? ? ? ? ?5.265204e-05 ?1.353213e-04
#> ARMA(1,1) + ARCH(1) ? ?5.415836e-05 ?8.983266e-05
#> ARCH(10) ? ? ? ? ? ? ? 5.417266e-05 ?8.975710e-05
#> ARMA(1,1) + GARCH(1,1) 5.339071e-05 ?9.244012e-05
我們可以觀察到,隨著模型復雜度的增加,樣本內(nèi)誤差趨于變?。ㄓ捎跀M合數(shù)據(jù)的自由度更高),盡管差異可以忽略不計。重要的實際上是樣本外誤差:我們可以看到,增加模型復雜度可能會得出較差的結果。就預測收益的誤差而言,似乎最簡單的iid模型已經(jīng)足夠了。
最后,讓我們展示一些樣本外誤差的圖表:
plot(error,
? ? main = "不同模型收益預測的樣本外誤差",

請注意,由于我們沒有重新擬合模型,因此隨著時間的發(fā)展,誤差越大(對于ARCH建模尤其明顯)。
滾動窗口比較
讓我們首先通過一個簡單的示例比較靜態(tài)預測與滾動預測的概念:
#ARMA(2,2)模型
spec <- spec(mean.model = list(armaOrder = c(2,2), include.mean = TRUE))
# 靜態(tài)擬合和預測
ar_static_fit <- fit(spec = spec, data = logreturns, out.sample = T_tst)
#滾動擬合和預測
modelroll <- aroll(spec = spec, data = logreturns, n.ahead = 1,
# 預測圖
plot(cbind("static forecast" ?= ar_static_fore_logreturns,
? ? main = "使用ARMA(2,2)模型進行預測", legend.loc = "topleft")
# 預測誤差圖
plot(error_logreturns, col = c("black", "red"), lwd = 2,
? ? main = "ARMA(2,2)模型的預測誤差", legend.loc = "topleft")


我們可以清楚地觀察到滾動窗口過程對時間序列的影響。
現(xiàn)在,我們可以在滾動窗口的基礎上重做所有模型的所有預測:
# 基于i.i.d.模型的滾動預測
roll(iid_spec, data = logreturns, n.ahead = 1, forecast.length = T_t
# AR(1)模型的滾動預測
roll(ar_spec, data = logreturns, n.ahead = 1, forecast.length = T_tst,
# ARMA(2,2)模型的滾動預測
roll(arma_spec, data = logreturns, n.ahead = 1, forecast.length = T_tst,
# ARMA(1,1)+ ARCH(1)模型的滾動預測
roll(arch_spec, data = logreturns, n.ahead = 1, forecast.length = T_tst,
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? refit.every = 50, refit.win
# ARMA(0,0)+ ARCH(10)模型的滾動預測
roll(long_arch_spec, data = logreturns, n.ahead = 1, forecast.length = T_tst,
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?refit.every = 50,
# ARMA(1,1)+ GARCH(1,1)模型的滾動預測
roll(garch_spec, data = logreturns, n.ahead = 1, forecast.length = T_tst,
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?refit.every = 50, refit.window
讓我們看看滾動基準情況下的預測誤差:
print(rolling_error_var)
#> ? ? ? ? ? ? ? ? ? ? ? ? ? in-sample out-of-sample
#> iid ? ? ? ? ? ? ? ? ? ?5.417266e-05 ?8.974166e-05
#> AR(1) ? ? ? ? ? ? ? ? ?5.414645e-05 ?9.038057e-05
#> ARMA(2,2) ? ? ? ? ? ? ?5.265204e-05 ?8.924223e-05
#> ARMA(1,1) + ARCH(1) ? ?5.415836e-05 ?8.991902e-05
#> ARCH(10) ? ? ? ? ? ? ? 5.417266e-05 ?8.976736e-05
#> ARMA(1,1) + GARCH(1,1) 5.339071e-05 ?8.895682e-05
和一些圖表:
plot(error_logreturns,
? ? main = "不同模型的滾動預測誤差", legend.loc = "topleft"

我們看到,現(xiàn)在所有模型都擬合了時間序列。此外,我們在模型之間沒有發(fā)現(xiàn)任何顯著差異。
我們最終可以比較靜態(tài)誤差和滾動誤差:
barplot(rbind(error_var[, "out-of-sample"], rolling_error_var[, "out-of-sample"])
? ? ? ?col = c("darkblue", "darkgoldenrod"),
? ? ? ?legend = c("靜態(tài)預測", "滾動預測"),

我們可以看到,滾動預測在某些情況下是必須的。因此,實際上,我們需要定期進行滾動預測改進。
方差模型
ARCH和GARCH模型
對數(shù)收益率殘差wt的ARCH(m)模型為
其中zt是具有零均值和恒定方差的白噪聲序列,而條件方差σ2t建模為
其中,m為模型階數(shù),ω> 0,αi≥0為參數(shù)。
GARCH(m,s)模型使用σ2t上的遞歸項擴展了ARCH模型:
其中參數(shù)ω> 0,αi≥0,βj≥0需要滿足∑mi =1αi+ ∑sj = 1βj≤1的穩(wěn)定性。
rugarch生成數(shù)據(jù)?
首先,我們需要定義模型:
# 指定具有給定系數(shù)和參數(shù)的GARCH模型
#>
#> *---------------------------------*
#> * ? ? ? GARCH Model Spec ? ? ? ? ?*
#> *---------------------------------*
#>
#> Conditional Variance Dynamics ? ?
#> ------------------------------------
#> GARCH Model ? ? ?: sGARCH(1,1)
#> Variance Targeting ? : FALSE
#>
#> Conditional Mean Dynamics
#> ------------------------------------
#> Mean Model ? ? ? : ARFIMA(1,0,0)
#> Include Mean ? ? : TRUE
#> GARCH-in-Mean ? ? ? ?: FALSE
#>
#> Conditional Distribution
#> ------------------------------------
#> Distribution : ?norm
#> Includes Skew ? ?: ?FALSE
#> Includes Shape ? : ?FALSE
#> Includes Lambda ?: ?FALSE
#> ? ? ? ? ? Level Fixed Include Estimate LB UB
#> mu ? ? ? ?0.005 ? ? 1 ? ? ? 1 ? ? ? ?0 NA NA
#> ar1 ? ? ?-0.900 ? ? 1 ? ? ? 1 ? ? ? ?0 NA NA
#> ma ? ? ? ?0.000 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> arfima ? ?0.000 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> archm ? ? 0.000 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> mxreg ? ? 0.000 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> omega ? ? 0.001 ? ? 1 ? ? ? 1 ? ? ? ?0 NA NA
#> alpha1 ? ?0.300 ? ? 1 ? ? ? 1 ? ? ? ?0 NA NA
#> beta1 ? ? 0.650 ? ? 1 ? ? ? 1 ? ? ? ?0 NA NA
#> gamma ? ? 0.000 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> eta1 ? ? ?0.000 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> eta2 ? ? ?0.000 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> delta ? ? 0.000 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> lambda ? ?0.000 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> vxreg ? ? 0.000 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> skew ? ? ?0.000 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> shape ? ? 0.000 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> ghlambda ?0.000 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> xi ? ? ? ?0.000 ? ? 0 ? ? ? 0 ? ? ? ?0 NA NA
#> $mu
#> [1] 0.005
#>
#> $ar1
#> [1] -0.9
#>
#> $omega
#> [1] 0.001
#>
#> $alpha1
#> [1] 0.3
#>
#> $beta1
#> [1] 0.65
true_params
#> ? ? mu ? ?ar1 ?omega alpha1 ?beta1
#> ?0.005 -0.900 ?0.001 ?0.300 ?0.650
然后,我們可以生成收益率時間序列:
# 模擬一條路徑
hpath(garch_spec, n.sim = T)
#> ?num [1:2000, 1] 0.167 -0.217
# 繪圖對數(shù)收益
{ plot(synth_log_returns, main = "GARCH模型的對數(shù)收益", lwd = 1.5)
?lines(synth_volatility
GARCH
現(xiàn)在,我們可以估計參數(shù):
# 指定一個GARCH模型
ugarchspec(mean.model = list(armaOrder = c(1,0)
# 估計模型
coef(garch_fit)
#> ? ? ? ? ? ?mu ? ? ? ? ? ar1 ? ? ? ? omega ? ? ? ?alpha1 ? ? ? ? beta1
#> ?0.0036510100 -0.8902333595 ?0.0008811434 ?0.2810460728 ?0.6717486402
#> ? ? mu ? ?ar1 ?omega alpha1 ?beta1
#> ?0.005 -0.900 ?0.001 ?0.300 ?0.650
# 系數(shù)誤差
#> ? ? ? ? ? mu ? ? ? ? ?ar1 ? ? ? ?omega ? ? ? alpha1 ? ? ? ?beta1
#> 0.0013489900 0.0097666405 0.0001188566 0.0189539272 0.0217486402
我們還可以研究樣本數(shù)量T對參數(shù)估計誤差的影響:
# 循環(huán)
for (T_ in T_sweep) {
?garch_fit
?error_coeffs_vs_T <- rbind(error_coeffs_vs_T, abs((coef(garch_fit) - true_params)/true_params))
?estim_coeffs_vs_T <- rbind(estim_coeffs_vs_T, coef(garch_fit))
# 繪圖
matplot(T_sweep, 100*error_coeffs_vs_T,
? ? ? ?main = "估計GARCH系數(shù)的相對誤差", xlab = "T", ylab = "誤差 (%)",
真實的ω幾乎為零,因此誤差非常不穩(wěn)定。至于其他系數(shù),就像在ARMA情況下一樣,μ的估計確實很差(相對誤差超過50%),而其他系數(shù)似乎在T = 800個樣本后得到了很好的估計。
GARCH結果比較?
作為健全性檢查,我們現(xiàn)在將比較兩個軟件包?fGarch?和?rugarch的結果:
# 指定具有特定參數(shù)值的ARMA(0,0)-GARCH(1,1)作為數(shù)據(jù)生成過程
garch_spec
#生成長度為1000的數(shù)據(jù)
path(garch_fixed_spec, n.sim = 1000)@path$
# 使用“ rugarch”包指定和擬合模型
rugarch_fit <- ugarchfit(spec = garch_spec, data = x)
# 使用包“ fGarch”擬合模型
garchFit(formula = ~ garch(1, 1), data = x, trace = FALSE)
# 比較模型系數(shù)
#> ? ? ? ? mu ? ? ?omega ? ? alpha1 ? ? ?beta1
#> 0.09749904 0.01395109 0.13510445 0.73938595
#> ? ? ? ? mu ? ? ?omega ? ? alpha1 ? ? ?beta1
#> 0.09750394 0.01392648 0.13527024 0.73971658
# 比較擬合的標準偏差
print(head(fGarch_fi
#> [1] 0.3513549 0.3254788 0.3037747 0.2869034 0.2735266 0.2708994
print(head(rugar
#> [1] 0.3538569 0.3275037 0.3053974 0.2881853 0.2745264 0.2716555
確實,這兩個軟件包給出了相同的結果。
使用rugarch包進行GARCH預測?
一旦估計出GARCH模型的參數(shù),就可以使用該模型預測未來的值。例如,基于過去的信息對條件方差的單步預測為
給定ω^ /(1-∑mi =1α^ i-∑sj =1β^ j)。軟件包?rugarch?使對樣本外數(shù)據(jù)的預測變得簡單:
# 估計模型,不包括樣本外
garch_fit
coef(garch_fit)
#> ? ? ? ? ? ?mu ? ? ? ? ? ar1 ? ? ? ? omega ? ? ? ?alpha1 ? ? ? ? beta1
#> ?0.0034964331 -0.8996287630 ?0.0006531088 ?0.3058756796 ?0.6815452241
# 預測整個樣本的對數(shù)收益
garch_fore@forecast$sigmaFor[1, ]
# 對數(shù)收益圖
plot(cbind("fitted" ? = fitted(garch_fit),
? ? main = "合成對數(shù)收益預測", legend.loc = "topleft")
#波動率對數(shù)收益圖
plot(cbind("fitted volatility" ? = sigma(garch_fit),
? ? main = "預測合成對數(shù)收益率的波動性", legend.loc = "topleft")
不同方法
讓我們首先加載S&P500:
# 加載標準普爾500指數(shù)數(shù)據(jù)
head(SP500_index_prices)
#> ? ? ? ? ? ? ?SP500
#> 2008-01-02 1447.16
#> 2008-01-03 1447.16
#> 2008-01-04 1411.63
#> 2008-01-07 1416.18
#> 2008-01-08 1390.19
#> 2008-01-09 1409.13
# 準備訓練和測試數(shù)據(jù)
x_trn <- x[1:T_trn]
x_tst <- x[-c(1:T_trn)]
# 繪圖
{ plot(x, main = "收益"
?addEventLines(xts("訓練", in
常數(shù)
讓我們從常數(shù)開始:
plot(cbind(sqrt(var_constant), x_trn)
? ? main = "常數(shù)")
移動平均值
現(xiàn)在,讓我們使用平方收益的移動平均值:
plot(cbind(sqrt(var_t), x_trn),
? ? main = "基于簡單滾動平方均值的包絡線(時間段=20)
EWMA
指數(shù)加權移動平均線(EWMA):
請注意,這也可以建模為ETS(A,N,N)狀態(tài)空間模型:
plot(cbind(std_t, x_trn),
? ? main = "基于平方EWMA的包絡")
乘法ETS
我們還可以嘗試ETS模型的不同變體。例如,具有狀態(tài)空間模型的乘性噪聲版本ETS(M,N,N):

plot(cbind(std_t, x_trn), col = c("red", "black")
? ? main = "基于平方的ETS(M,N,N)的包絡"

ARCH
現(xiàn)在,我們可以使用更復雜的ARCH建模:

plot(cbind(std_t, x_trn), col = c("red", "black")
? ? main = "基于ARCH(5)的包絡")

GARCH
我們可以將模型提升到GARCH:

plot(cbind(std_t, x_trn), col = c("red", "black")
? ? main = "基于GARCH(1,1)的包絡")

SV隨機波動率
最后,我們可以使用隨機波動率建模:

或者,等效地,

plot(cbind(std_t, x_trn), col = c("red", "black"),
? ? main = "基于隨機波動率的包絡分析")

比較
現(xiàn)在,我們可以比較每種方法在樣本外期間的方差估計中的誤差:
#> ? ? ? ? ? MA ? ? ? ? EWMA ? ETS(M,N,N) ? ? ?ARCH(5) ? GARCH(1,1) ? ? ? ? ? SV
#> 2.204965e-05 7.226188e-06 3.284057e-06 7.879039e-05 6.496545e-06 6.705059e-06
barplot(error_all, main = "樣本外方差估計中的誤差"

滾動窗口比較
六種方法的滾動窗口比較:MA,EWMA,ETS(MNN),ARCH(5),GARCH(1,1)和SV。
?
#滾動窗口
lookback <- 200
len_tst <- 40
for (i in seq(lookback, T-len_tst, by = len_tst)) {
?# MA
?var_t <- roll_meanr(x_trn^2, n = 20, fill = NA)
?var_fore <- var(x_trn/sqrt(var_t), na.rm = TRUE) * tail(var_t, 1)
?error_ma <- c(error_ma, abs(var_fore - var_tst))
?
?# EWMA
?error_ewma <- c(error_ewma, abs(var_fore - var_tst))
?
?# ETS(M,N,N)
?error_ets_mnn <- c(error_ets_mnn, abs(var_fore - var_tst))
?
?# ARCH
?error_arch <- c(error_arch, abs(var_fore - var_tst))
?
?# GARCH
?error_garch <- c(error_garch, abs(var_fore - var_tst))
?
?# SV
?error_sv <- c(error_sv, abs(var_fore - var_tst))
}
barplot(error_all, main = "方差估計誤差",

多元GARCH模型
出于說明目的,我們將僅考慮恒定條件相關(CCC)和動態(tài)條件相關(DCC)模型,因為它們是最受歡迎的模型。對數(shù)收益率殘差wt建模為

其中zt是具有零均值和恒定協(xié)方差矩陣II的iid白噪聲序列。條件協(xié)方差矩陣Σt建模為

其中Dt = Diag(σ1,t,...,σN,t)是標準化噪聲向量C,協(xié)方差矩陣ηt=C-1wt(即,它包含等于1的對角線元素)。
基本上,使用此模型,對角矩陣Dt包含一組單變量GARCH模型,然后矩陣C包含序列之間的一些相關性。該模型的主要缺點是矩陣C是恒定的。為了克服這個問題,DCC被提議為

其中Ct包含等于1的對角元素。要強制等于1的對角元素,Engle將其建模為

Qt具有任意對角線元素并遵循模型

我們將生成數(shù)據(jù),估計參數(shù)和預測。
從加載多元ETF數(shù)據(jù)開始:
SPDR S&P 500 ETF
20年以上國債ETF
IEF:7-10年期國債ETF
# 下載數(shù)據(jù)
prices <- xts()
head(prices)
#> ? ? ? ? ? ? ? ? SPY ? ? ?TLT ? ? ?IEF
#> 2013-01-02 127.8779 99.85183 93.65224
#> 2013-01-03 127.5890 98.49886 93.17085
#> 2013-01-04 128.1493 98.88306 93.21463
#> 2013-01-07 127.7991 98.92480 93.26714
#> 2013-01-08 127.4314 99.57622 93.49468
#> 2013-01-09 127.7553 99.48438 93.54719
# 繪制三個對數(shù)價格序列
plot(log(prices)
? ? main = "三個ETF的對數(shù)價格", legend.loc = "topleft")

首先,我們定義模型:
# 指定i.i.d.單變量時間序列模型
ugarch_spec
# 指定DCC模型
spec( multispec(replicate(spec, n = 3))
? ?
接下來,我們擬合模型:
# 估計模型
#>
#> *---------------------------------*
#> * ? ? ? ? ?DCC GARCH Fit ? ? ? ? ?*
#> *---------------------------------*
#>
#> Distribution ? ? ? ? : ?mvnorm
#> Model ? ? ? ? ? ? ? ?: ?DCC(1,1)
#> No. Parameters ? ? ? : ?44
#> [VAR GARCH DCC UncQ] : [30+9+2+3]
#> No. Series ? ? ? ? ? : ?3
#> No. Obs. ? ? ? ? ? ? : ?1007
#> Log-Likelihood ? ? ? : ?12198.4
#> Av.Log-Likelihood ? ?: ?12.11
#>
#> Optimal Parameters
#> -----------------------------------
#> ? ? ? ? ? ? ? Estimate ?Std. Error ? t value Pr(>|t|)
#> [SPY].omega ? 0.000004 ? ?0.000000 ?11.71585 0.000000
#> [SPY].alpha1 ?0.050124 ? ?0.005307 ? 9.44472 0.000000
#> [SPY].beta1 ? 0.870051 ? ?0.011160 ?77.96041 0.000000
#> [TLT].omega ? 0.000001 ? ?0.000001 ? 0.93156 0.351563
#> [TLT].alpha1 ?0.019716 ? ?0.010126 ? 1.94707 0.051527
#> [TLT].beta1 ? 0.963760 ? ?0.006434 149.79210 0.000000
#> [IEF].omega ? 0.000000 ? ?0.000001 ? 0.46913 0.638979
#> [IEF].alpha1 ?0.031741 ? ?0.023152 ? 1.37097 0.170385
#> [IEF].beta1 ? 0.937777 ? ?0.016498 ?56.84336 0.000000
#> [Joint]dcca1 ?0.033573 ? ?0.014918 ? 2.25044 0.024421
#> [Joint]dccb1 ?0.859787 ? ?0.079589 ?10.80278 0.000000
#>
#> Information Criteria
#> ---------------------
#> ? ? ? ? ? ? ? ? ? ?
#> Akaike ? ? ? -24.140
#> Bayes ? ? ? ?-23.925
#> Shibata ? ? ?-24.143
#> Hannan-Quinn -24.058
#>
#>
#> Elapsed time : 0.8804049
我們可以繪制時變相關性:
# 提取時變協(xié)方差和相關矩陣
dim(dcc_cor)
#> [1] ? ?3 ? ?3 1007
#繪圖
plot(corr_t
? ? main = "時變相關", legend.loc = "left")

我們看到兩個收益ETF之間的相關性非常高且相當穩(wěn)定。與SPY的相關性較小,在小于0的區(qū)間波動。
?

最受歡迎的見解
1.HAR-RV-J與遞歸神經(jīng)網(wǎng)絡(RNN)混合模型預測和交易大型股票指數(shù)的高頻波動率
2.R語言中基于混合數(shù)據(jù)抽樣(MIDAS)回歸的HAR-RV模型預測GDP增長
3.波動率的實現(xiàn):ARCH模型與HAR-RV模型
4.R語言ARMA-EGARCH模型、集成預測算法對SPX實際波動率進行預測
5.GARCH(1,1),MA以及歷史模擬法的VaR比較
6.R語言多元COPULA GARCH 模型時間序列預測
7.R語言基于ARMA-GARCH過程的VAR擬合和預測
8.matlab預測ARMA-GARCH 條件均值和方差模型
9.R語言對S&P500股票指數(shù)進行ARIMA + GARCH交易策略