在 R 中估計 GARCH 參數存在的問題
本文翻譯自《Problems In Estimating GARCH Parameters in R 》
原文鏈接：https://ntguardian.wordpress.com/2017/11/02/problems-estimating-garch-parameters-r/
更新（11/2/17 3:00 PM MDT）：我從 R 的金融板塊郵件列表收到一位知名金融工具包貢獻者——Brian Peterson 的郵件：
我強烈建議關注 rugarch 和 rmgarch。RMetrics 序列包的主要維護者 Diethelm Wuertz 在 2016 年死於車禍，目前的代碼基本處於無維護狀態。
我會看看這是否解決了這個問題。謝謝 Brian！我把這篇帖子留給其他人，作為將來避免使用 gGarch 的告誡。這對我來說是個新聞，因為書籍經常引用 fGarch，所以這可能是那些尋求在 R 中使用 GARCH 模型的人的資源——為什么不要使用 fGarch。
更新（11/2/17 11:30 PM MDT）：我用 rugarch 進行了一次快速實驗，看起來它同樣被這個問題困擾。下面是我運行的代碼，我會盡快在明天貼出一份全面的研究。
library(rugarch)
spec = ugarchspec(
variance.model = list(garchOrder = c(1, 1)),
mean.model = list(
armaOrder = c(0, 0), include.mean = FALSE),
fixed.pars = list(
alpha1 = 0.2, beta1 = 0.2, omega = 0.2))
ugarchpath(
spec = spec, n.sim = 1000, n.start = 1000) -> x
srs = x@path$seriesSim
spec1 = ugarchspec(
variance.model = list(garchOrder = c(1, 1)),
mean.model = list(
armaOrder = c(0, 0), include.mean = FALSE))
ugarchfit(spec = spec1, data = srs)
ugarchfit(spec = spec1, data = srs[1:100])
這些天我的研究重點是變點檢測方法。這些是用於檢測數據序列中出現結構性變化的統計檢驗和過程。來自質量控制的早期示例是在生產小部件時檢測機器是否未校准。可能存在一些感興趣的測量值，例如我們觀察到的滾珠軸承的直徑。機器按順序生成這些小部件。在原假設下，滾珠軸承的平均直徑不會改變，而在備擇假設中，在制造過程中的某些未知點處，機器變得未校准並且滾珠軸承的平均直徑發生變化。然后，檢驗在這兩個假設之間做出決定。
這些類型的檢驗對經濟學家和金融業工作者也很重要，特別是對於預測。我們再次使用按時間索引的序列數據，我首選的例子是股票的價格，在給出了股票的常見時間序列圖，人們可以立即將其識別為時間序列，但還有更多的數據集，如州的 GDP 或失業率。經濟學家希望使用過去的數據和統計量預測這些數量。統計方法的一個假設是要預測的序列是平穩的：數據的產生過程具有單一的均值、自相關性、分布等。這個假設並不總是經過檢驗，但對成功的預測至關重要。結構性變化的檢驗檢查了這個假設，如果結果是假的，預測者可能需要在訓練他們的模型時分割他們的數據集。
我之前寫過關於這些檢驗的文章，介紹了 CUSUM 統計量，這是檢測結構性變化的最流行的統計量。我的導師和他以前的博士生（Greg Rice，目前是滑鐵盧大學的教授）開發了一個新的檢驗統計量，可以更好地檢測數據集中早期或晚期發生的結構性變化（想象一下，制作小部件的機器幾乎沒有被校准，一百個小部件中只有最后一打樣本受到影響）。我們提交論文的期刊正在要求我們進行修訂，其中一個修訂是更好的示例應用（我們最初使用上述博客文章中討論的工資 / 生產率數據，審稿人抱怨這些變量是被相同因素決定的（codetermined），所以使用一個對另一個做回歸沒意義，我不同意這樣的抱怨，但我不會在這里爭論）。
我們希望將我們的檢驗應用於檢測 GARCH 模型中的結構性變化，這是金融時間序列中的常見模型。據我所知，用於 GARCH 模型估計和推斷（以及其他工作）的“最新技術” R 包是 fGarch。特別是，函數 garchFit() 用於從數據中估計 GARCH 模型。但是，當我們嘗試在我們的檢驗中使用此函數時，我們得到了明顯病態的數值（我們已經完成了模擬研究以了解預期的行為）。沒有變化的原假設在模擬序列上被完全拒絕。我從未看到檢驗未能拒絕原假設，即使原假設總是對的。即使樣本量為 10000，幾乎不是小樣本也是如此。
我們認為問題可能在於參數估計的協方差矩陣的估計，並且我煞費苦心地推導和編寫函數以使該矩陣不使用數值微分，但這並沒有阻止不良行為。最后我的導師和我上周三戲弄了 garchFit() 一把，並發現該函數應該受到責備。當我估計參數（不一定是我們最初認為的協方差矩陣，盡管它可能也被污染）時，函數對模擬數據的行為是如此不穩定，依我來看，該函數基本上是無用的。
這個函數應該是眾所周知的，問題當然可能在於我，而不是 fGarch（或者可能有更好的包）。這個函數是如此重要，讓我感到我應該分享我的發現。在這篇文章中，我展示了一序列數值實驗，證明了 garchFit() 的病態行為。
GARCH 模型基礎
\(\text{GARCH}(1,1)\) 是一種時間序列模型，用於為金融工具（例如股票）回報的波動率建模。用 \(\epsilon_t\) 表示 \(\text{GARCH}(1,1)\) 過程，這可以表示諸如股票回報的方差。\(\text{GARCH}(1,1)\) 模型（無均值參數）以遞歸的形式表示：
\(\sigma_t\) 是隨機過程的條件標准差，也就是條件波動率，\(\eta_t\) 是一個隨機過程。
關注金融的人[1]注意到金融工具（如股票或共同基金）的回報表現出稱為波動率聚集的行為。有些時期金融工具相對溫和，沒有劇烈的市場走勢。在其他情況下，金融工具的價格可能會大幅波動，這些時期不是一次性的單日行為，而是可以持續一段時間。GARCH 模型來建模波動率聚集。
一些人認為，即使股票的日常變動基本上是不可預測的(股票在任何一天都等概率地被高估或低估)，但波動率是可預測的。即使對那些不狂熱相信未來的回報可以預測的人來說，這些模型依然很重要。例如，如果使用模型 \(R_t = \alpha + \beta R_{M,t} + \epsilon_t\) 來估計股票的 beta 統計量（其中 \(R_t\) 是股票在時間 \(t\) 上的回報；\(R_{M,t}\) 是市場回報；\(\epsilon_t\) 是“隨機噪音”），\(\epsilon_t\) 很可能不是 i.i.d 的隨機數序列（通常在其他上下文中假設），而實際上是一個 GARCH 序列。建模者想知道在這種情況下她的估計值的行為。因此，GARCH 模型被認為是重要的。實際上，我剛剛描述的波動率聚集行為有時被描述為“GARCH行為”，因為它經常出現，GARCH 模型是解決它們的常用工具。（首字母縮略詞 GARCH 代表廣義自回歸條件異方差性，它是對波動率時變性的統計學語言描述。）
\(\eta_t\) 可以是任何隨機過程，但經常選擇使用 i.i.d 的標准正態隨機變量序列。這里 \(\eta_t\) 是模型中唯一的隨機源。為了讓 \(\text{GARCH}(1,1)\) 過程有一個穩定解，我們必須要求 \(\alpha + \beta > 0\)。在這種情況下，該過程具有 \(\frac{\omega} {1 - \alpha - \beta}\) 的長期方差。
估計 GARCH 參數
我上面寫的過程是一個無限過程。索引 \(t\) 可以一直擴展到負數。顯然在實踐中我們沒有觀察到無限序列，因此如果我們想在實踐中使用 \(\text{GARCH}(1,1)\) 模型，我們需要考慮類似的序列：
下面是新序列的獨家秘方：
我們選擇這個序列的初始值（前面描述的理論 \(\text{GARCH}(1,1)\) 序列沒有初始值）！這個序列非常類似於理論序列，但它的整體上是可觀察的，並且可以證明使用該序列估計的參數非常接近理論上無限 \(\text{GARCH}(1,1)\) 過程的參數。
當然，這些過程最重要的任務之一就是估算它們的參數。對於 \(\text{GARCH}(1,1)\) 過程，這些是 \(\omega\)、\(\alpha\) 和 \(\beta\)。一種基本方法是找到擬最大似然估計（QMLE）的估計量。假設我們有 \(t-1\) 時的觀察值。在 QMLE 中，當假設 \(\eta_t\) 遵循標准正態分布（即 \(\eta_t \sim N(0,1)\)）時，我們使用 \(\epsilon_t\) 的條件分布。我們假設一直到 \(t-1\)，整個過程的歷史記錄已知，這意味着 \(\tilde{\sigma}_t\) 也是已知的(實際上我們需要知道的是 \(t-1\) 時間的值，但我離題了)。在這種情況下，我們有 \(\epsilon_t \sim N(0，\tilde{\sigma}_t^2)\)。設 \((x | \mathscr{f}_t) = f(x | \tilde{\sigma}_t^2)\) 是 \(\epsilon_t\) 的條件分布（所以 \(f(x | \tilde{\sigma}_t^2) = \frac{1}{\sqrt{2 \pi \tilde{\sigma}_t^2}} \exp \left(-\frac{x^2}{2 \tilde{\sigma}_t^2} \right)\)）。然后是擬似然方程
像大多數似然方法一樣，統計學家不是直接優化擬似然函數，而是嘗試優化對數似然函數，\(\log \left(\mathscr{L}_T(\epsilon_1，...，\epsilon_T)\right)\)，經過一些工作后，不難發現這等同於最小化
請注意，\(\omega\)、\(\alpha\) 和 \(\beta\) 通過 \(\tilde{\sigma}_t^2\) 包含在其中。對於最小化此數量的參數，沒有閉式解。這意味着必須應用數值優化技術來找到參數。
可以證明，當以這種方式計算時，參數 \(\omega\)、\(\alpha\) 和 \(\beta\) 的估計量是一致的（意味着它們依概率收斂到它們的真實值），並且漸近地服從隨高斯分布。[2]這些屬性我們可以與樣本均值相聯系，但是我們可能樂觀地認為這些估計量的收斂速度與樣本均值的收斂速度一樣好，我們可以預期可比較的漸近行為。
理想情況下，參數應該像下面說明的過程一樣。
library(ggplot2)
x <- rnorm(1000, sd = 1/3)
df <- t(
sapply(
50:1000,
function(t)
{
return(c("mean" = mean(x[1:t]), "mean.se" = sd(x[1:t])/sqrt(t)))
}))
df <- as.data.frame(df)
df$t <- 50:1000
ggplot(
df,
aes(x = t, y = mean)) +
geom_line() +
geom_ribbon(
aes(x = t, ymin = mean - 2 * mean.se, ymax = mean + 2 * mean.se),
color = "grey", alpha = 0.5) +
geom_hline(color = "blue", yintercept = 0) +
coord_cartesian(ylim = c(-0.5, 0.5))
fGarch 參數估計的行為
在繼續之前，讓我們生成 \(\text{GARCH}(1,1)\) 序列。在本文中，我使用了所有參數都等於 0.2 的過程。注意，對於 \(\text{GARCH}(1,1)\) 過程，長期方差將為 \(\frac{0.2}{1 - 0.2 - 0.2} = \frac{1}{3}\)。
set.seed(110117)
library(fGarch)
x <- garchSim(
garchSpec(
model = list(
"alpha" = 0.2, "beta" = 0.2, "omega" = 0.2)),
n.start = 1000,
n = 1000)
plot(x)
讓我們看看 fGarch 的函數 garchFit() 所使用的參數。
args(garchFit)
## function (formula = ~garch(1, 1), data = dem2gbp, init.rec = c("mci",
## "uev"), delta = 2, skew = 1, shape = 4, cond.dist = c("norm",
## "snorm", "ged", "sged", "std", "sstd", "snig", "QMLE"), include.mean = TRUE,
## include.delta = NULL, include.skew = NULL, include.shape = NULL,
## leverage = NULL, trace = TRUE, algorithm = c("nlminb", "lbfgsb",
## "nlminb+nm", "lbfgsb+nm"), hessian = c("ropt", "rcd"),
## control = list(), title = NULL, description = NULL, ...)
## NULL
該函數提供了一些選項，要最大化的分布（cond.dist）和用於優化的算法（algorithm）。除非另有說明，否則我將始終選擇 cond.dist = QMLE 來指示函數使用 QMLE 估算器。
這是一次調用。
garchFit(
data = x, cond.dist = "QMLE", include.mean = FALSE)
##
## Series Initialization:
## ARMA Model: arma
## Formula Mean: ~ arma(0, 0)
## GARCH Model: garch
## Formula Variance: ~ garch(1, 1)
## ARMA Order: 0 0
## Max ARMA Order: 0
## GARCH Order: 1 1
## Max GARCH Order: 1
## Maximum Order: 1
## Conditional Dist: QMLE
## h.start: 2
## llh.start: 1
## Length of Series: 1000
## Recursion Init: mci
## Series Scale: 0.5320977
##
## Parameter Initialization:
## Initial Parameters: $params
## Limits of Transformations: $U, $V
## Which Parameters are Fixed? $includes
## Parameter Matrix:
## U V params includes
## mu -0.15640604 0.156406 0.0 FALSE
## omega 0.00000100 100.000000 0.1 TRUE
## alpha1 0.00000001 1.000000 0.1 TRUE
## gamma1 -0.99999999 1.000000 0.1 FALSE
## beta1 0.00000001 1.000000 0.8 TRUE
## delta 0.00000000 2.000000 2.0 FALSE
## skew 0.10000000 10.000000 1.0 FALSE
## shape 1.00000000 10.000000 4.0 FALSE
## Index List of Parameters to be Optimized:
## omega alpha1 beta1
## 2 3 5
## Persistence: 0.9
##
##
## --- START OF TRACE ---
## Selected Algorithm: nlminb
##
## R coded nlminb Solver:
##
## 0: 1419.0152: 0.100000 0.100000 0.800000
## 1: 1418.6616: 0.108486 0.0998447 0.804683
## 2: 1417.7139: 0.109746 0.0909961 0.800931
## 3: 1416.7807: 0.124977 0.0795152 0.804400
## 4: 1416.7215: 0.141355 0.0446605 0.799891
## 5: 1415.5139: 0.158059 0.0527601 0.794304
## 6: 1415.2330: 0.166344 0.0561552 0.777108
## 7: 1415.0415: 0.195230 0.0637737 0.743465
## 8: 1415.0031: 0.200862 0.0576220 0.740088
## 9: 1414.9585: 0.205990 0.0671331 0.724721
## 10: 1414.9298: 0.219985 0.0713468 0.712919
## 11: 1414.8226: 0.230628 0.0728325 0.697511
## 12: 1414.4689: 0.325750 0.0940514 0.583114
## 13: 1413.4560: 0.581449 0.143094 0.281070
## 14: 1413.2804: 0.659173 0.157127 0.189282
## 15: 1413.2136: 0.697840 0.155964 0.150319
## 16: 1413.1467: 0.720870 0.142550 0.137645
## 17: 1413.1416: 0.726527 0.138146 0.135966
## 18: 1413.1407: 0.728384 0.137960 0.134768
## 19: 1413.1392: 0.731725 0.138321 0.132991
## 20: 1413.1392: 0.731146 0.138558 0.133590
## 21: 1413.1392: 0.730849 0.138621 0.133850
## 22: 1413.1392: 0.730826 0.138622 0.133869
##
## Final Estimate of the Negative LLH:
## LLH: 782.211 norm LLH: 0.782211
## omega alpha1 beta1
## 0.2069173 0.1386221 0.1338686
##
## R-optimhess Difference Approximated Hessian Matrix:
## omega alpha1 beta1
## omega -8858.897 -1839.6144 -2491.9827
## alpha1 -1839.614 -782.8005 -531.7393
## beta1 -2491.983 -531.7393 -729.7246
## attr(,"time")
## Time difference of 0.04132652 secs
##
## --- END OF TRACE ---
##
##
## Time to Estimate Parameters:
## Time difference of 0.3866439 secs
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(data = x, cond.dist = "QMLE", include.mean = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 1)
## <environment: 0xa636ba4>
## [data = x]
##
## Conditional Distribution:
## QMLE
##
## Coefficient(s):
## omega alpha1 beta1
## 0.20692 0.13862 0.13387
##
## Std. Errors:
## robust
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## omega 0.20692 0.05102 4.056 5e-05 ***
## alpha1 0.13862 0.04928 2.813 0.00491 **
## beta1 0.13387 0.18170 0.737 0.46128
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## -782.211 normalized: -0.782211
##
## Description:
## Thu Nov 2 13:01:14 2017 by user:
參數不接近真實參數。人們可能最初將其歸因於隨機性，但事實似乎並非如此。
例如，當我在前 500 個數據點上擬合模型時，我能得到什么呢？
garchFit(
data = x[1:500], cond.dist = "QMLE", include.mean = FALSE)
##
## Series Initialization:
## ARMA Model: arma
## Formula Mean: ~ arma(0, 0)
## GARCH Model: garch
## Formula Variance: ~ garch(1, 1)
## ARMA Order: 0 0
## Max ARMA Order: 0
## GARCH Order: 1 1
## Max GARCH Order: 1
## Maximum Order: 1
## Conditional Dist: QMLE
## h.start: 2
## llh.start: 1
## Length of Series: 500
## Recursion Init: mci
## Series Scale: 0.5498649
##
## Parameter Initialization:
## Initial Parameters: $params
## Limits of Transformations: $U, $V
## Which Parameters are Fixed? $includes
## Parameter Matrix:
## U V params includes
## mu -0.33278068 0.3327807 0.0 FALSE
## omega 0.00000100 100.0000000 0.1 TRUE
## alpha1 0.00000001 1.0000000 0.1 TRUE
## gamma1 -0.99999999 1.0000000 0.1 FALSE
## beta1 0.00000001 1.0000000 0.8 TRUE
## delta 0.00000000 2.0000000 2.0 FALSE
## skew 0.10000000 10.0000000 1.0 FALSE
## shape 1.00000000 10.0000000 4.0 FALSE
## Index List of Parameters to be Optimized:
## omega alpha1 beta1
## 2 3 5
## Persistence: 0.9
##
##
## --- START OF TRACE ---
## Selected Algorithm: nlminb
##
## R coded nlminb Solver:
##
## 0: 706.37230: 0.100000 0.100000 0.800000
## 1: 706.27437: 0.103977 0.100309 0.801115
## 2: 706.19091: 0.104824 0.0972295 0.798477
## 3: 706.03116: 0.112782 0.0950253 0.797812
## 4: 705.77389: 0.122615 0.0858136 0.788169
## 5: 705.57316: 0.134608 0.0913105 0.778144
## 6: 705.43424: 0.140011 0.0967118 0.763442
## 7: 705.19541: 0.162471 0.102711 0.739827
## 8: 705.16325: 0.166236 0.0931680 0.737563
## 9: 705.09943: 0.168962 0.100977 0.731085
## 10: 704.94924: 0.203874 0.0958205 0.702986
## 11: 704.78210: 0.223975 0.108606 0.664678
## 12: 704.67414: 0.250189 0.122959 0.630886
## 13: 704.60673: 0.276532 0.131788 0.595346
## 14: 704.52185: 0.335952 0.146435 0.520961
## 15: 704.47725: 0.396737 0.157920 0.448557
## 16: 704.46540: 0.442499 0.164111 0.396543
## 17: 704.46319: 0.440935 0.161566 0.400606
## 18: 704.46231: 0.442951 0.159225 0.400940
## 19: 704.46231: 0.443022 0.159284 0.400863
## 20: 704.46230: 0.443072 0.159363 0.400851
## 21: 704.46230: 0.443112 0.159367 0.400807
##
## Final Estimate of the Negative LLH:
## LLH: 405.421 norm LLH: 0.810842
## omega alpha1 beta1
## 0.1339755 0.1593669 0.4008074
##
## R-optimhess Difference Approximated Hessian Matrix:
## omega alpha1 beta1
## omega -8491.005 -1863.4127 -2488.5700
## alpha1 -1863.413 -685.6071 -585.4327
## beta1 -2488.570 -585.4327 -744.1593
## attr(,"time")
## Time difference of 0.02322888 secs
##
## --- END OF TRACE ---
##
##
## Time to Estimate Parameters:
## Time difference of 0.1387401 secs
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(data = x[1:500], cond.dist = "QMLE", include.mean = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 1)
## <environment: 0xa85f084>
## [data = x[1:500]]
##
## Conditional Distribution:
## QMLE
##
## Coefficient(s):
## omega alpha1 beta1
## 0.13398 0.15937 0.40081
##
## Std. Errors:
## robust
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## omega 0.13398 0.11795 1.136 0.2560
## alpha1 0.15937 0.07849 2.030 0.0423 *
## beta1 0.40081 0.44228 0.906 0.3648
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## -405.421 normalized: -0.810842
##
## Description:
## Thu Nov 2 13:01:15 2017 by user:
請注意，參數 \(\beta\)（列為 beta1）發生了巨大變化。不同的截止點又會怎么樣？
garchFit(
data = x[1:200], cond.dist = "QMLE", include.mean = FALSE)
##
## Series Initialization:
## ARMA Model: arma
## Formula Mean: ~ arma(0, 0)
## GARCH Model: garch
## Formula Variance: ~ garch(1, 1)
## ARMA Order: 0 0
## Max ARMA Order: 0
## GARCH Order: 1 1
## Max GARCH Order: 1
## Maximum Order: 1
## Conditional Dist: QMLE
## h.start: 2
## llh.start: 1
## Length of Series: 200
## Recursion Init: mci
## Series Scale: 0.5746839
##
## Parameter Initialization:
## Initial Parameters: $params
## Limits of Transformations: $U, $V
## Which Parameters are Fixed? $includes
## Parameter Matrix:
## U V params includes
## mu -0.61993813 0.6199381 0.0 FALSE
## omega 0.00000100 100.0000000 0.1 TRUE
## alpha1 0.00000001 1.0000000 0.1 TRUE
## gamma1 -0.99999999 1.0000000 0.1 FALSE
## beta1 0.00000001 1.0000000 0.8 TRUE
## delta 0.00000000 2.0000000 2.0 FALSE
## skew 0.10000000 10.0000000 1.0 FALSE
## shape 1.00000000 10.0000000 4.0 FALSE
## Index List of Parameters to be Optimized:
## omega alpha1 beta1
## 2 3 5
## Persistence: 0.9
##
##
## --- START OF TRACE ---
## Selected Algorithm: nlminb
##
## R coded nlminb Solver:
##
## 0: 280.63354: 0.100000 0.100000 0.800000
## 1: 280.63302: 0.100315 0.100088 0.800223
## 2: 280.63262: 0.100695 0.0992822 0.800059
## 3: 280.63258: 0.102205 0.0983397 0.800404
## 4: 280.63213: 0.102411 0.0978709 0.799656
## 5: 280.63200: 0.102368 0.0986702 0.799230
## 6: 280.63200: 0.101930 0.0984977 0.800005
## 7: 280.63200: 0.101795 0.0983937 0.799987
## 8: 280.63197: 0.101876 0.0984197 0.799999
## 9: 280.63197: 0.102003 0.0983101 0.799965
## 10: 280.63197: 0.102069 0.0983780 0.799823
## 11: 280.63197: 0.102097 0.0983703 0.799827
## 12: 280.63197: 0.102073 0.0983592 0.799850
## 13: 280.63197: 0.102075 0.0983616 0.799846
##
## Final Estimate of the Negative LLH:
## LLH: 169.8449 norm LLH: 0.8492246
## omega alpha1 beta1
## 0.03371154 0.09836156 0.79984610
##
## R-optimhess Difference Approximated Hessian Matrix:
## omega alpha1 beta1
## omega -26914.901 -6696.498 -8183.925
## alpha1 -6696.498 -2239.695 -2271.547
## beta1 -8183.925 -2271.547 -2733.098
## attr(,"time")
## Time difference of 0.02161336 secs
##
## --- END OF TRACE ---
##
##
## Time to Estimate Parameters:
## Time difference of 0.09229803 secs
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(data = x[1:200], cond.dist = "QMLE", include.mean = FALSE)
##
## Mean and Variance Equation:
## data ~ garch(1, 1)
## <environment: 0xad38a84>
## [data = x[1:200]]
##
## Conditional Distribution:
## QMLE
##
## Coefficient(s):
## omega alpha1 beta1
## 0.033712 0.098362 0.799846
##
## Std. Errors:
## robust
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## omega 0.03371 0.01470 2.293 0.0218 *
## alpha1 0.09836 0.04560 2.157 0.0310 *
## beta1 0.79985 0.03470 23.052 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## -169.8449 normalized: -0.8492246
##
## Description:
## Thu Nov 2 13:01:15 2017 by user:
對於 200 個觀察，\(\beta\) 的估計是巨大的，標准差相對較小！
讓我們深入探討這一點。我在猶他大學數學系的超級計算機上進行了一些數值實驗（譯注：實際上，普通家用電腦也能應付）。下面是一個輔助函數，用於通過 garchFit()（在計算過程中屏蔽所有 garchFit() 的輸出）來提取特定擬合的系數和標准差。
getFitData <- function(x,
cond.dist = "QMLE",
include.mean = FALSE,
...)
{
args <- list(...)
args$data = x
args$cond.dist = cond.dist
args$include.mean = include.mean
log <- capture.output(
{
fit <- do.call(garchFit, args = args)
})
res <- coef(fit)
res[paste0(names(fit@fit$se.coef), ".se")] <- fit@fit$se.coef
return(res)
}
第一個實驗是在每個可能的截止點計算該特定序列的系數。
(在編寫此文檔時，不會評估以下代碼塊。我已將結果保存在 Rda 文件中。對於涉及並行計算的每個代碼塊都是如此。我在猶他大學數學系的超級計算機上執行了這些計算，在這里保存結果。)
library(doParallel)
set.seed(110117)
cl <- makeCluster(detectCores() - 1)
registerDoParallel(cl)
x <- garchSim(
garchSpec(
model = list(
alpha = 0.2, beta = 0.2, omega = 0.2)),
n.start = 1000, n = 1000)
params <- foreach(
t = 50:1000,
.combine = rbind,
.packages = c("fGarch")) %dopar%
{
getFitData(x[1:t])
}
rownames(params) <- 50:1000
下面我繪制這些系數，以及對應於兩倍標准差的區域。該區域應大致對應 95％ 的置信區間。
params_df <- as.data.frame(params)
params_df$t <- as.numeric(rownames(params))
ggplot(params_df) +
geom_line(
aes(x = t, y = beta1)) +
geom_hline(
yintercept = 0.2, color = "blue") +
geom_ribbon(
aes(x = t,
ymin = beta1 - 2 * beta1.se,
ymax = beta1 + 2 * beta1.se),
color = "grey", alpha = 0.5) +
ylab(expression(hat(beta))) +
scale_y_continuous(
breaks = c(0, 0.2, 0.25, 0.5, 1)) +
coord_cartesian(ylim = c(0, 1))
這是一幅令人震驚的畫面（但不是我見過的最令人震驚的畫面，這是最好的案例之一）。請注意，置信區間無法覆蓋 \(\beta = 0.2\) 的真實值，直到大約 375 個數據點。這些間隔本應該在大約 95％ 的時間內包含真實值！除此之外，置信區間相當大。
讓我們看看其他參數的行為。
library(reshape2)
library(plyr)
library(dplyr)
param_reshape <- function(p)
{
p <- as.data.frame(p)
p$t <- as.integer(rownames(p))
pnew <- melt(
p, id.vars = "t", variable.name = "parameter")
pnew$parameter <- as.character(pnew$parameter)
pnew.se <- pnew[grepl("*.se", pnew$parameter), ]
pnew.se$parameter <- sub(".se", "", pnew.se$parameter)
names(pnew.se)[3] <- "se"
pnew <- pnew[!grepl("*.se", pnew$parameter), ]
return(
join(
pnew, pnew.se,
by = c("t", "parameter"),
type = "inner"))
}
ggp <- ggplot(
param_reshape(params),
aes(x = t, y = value)) +
geom_line() +
geom_ribbon(
aes(ymin = value - 2 * se,
ymax = value + 2 * se),
color = "grey",
alpha = 0.5) +
geom_hline(yintercept = 0.2, color = "blue") +
scale_y_continuous(
breaks = c(0, 0.2, 0.25, 0.5, 0.75, 1)) +
coord_cartesian(ylim = c(0, 1)) +
facet_grid(. ~ parameter)
print(ggp + ggtitle("NLMINB Optimization"))
這種現象不僅限於 \(\beta\)，\(\omega\) 也表現出不良行為。(\(\alpha\) 也不是很好，但要好得多。)
這種行為並不罕見，這是典型的。下面是使用不同種子生成的類似序列的圖。
seeds <- c(
103117, 123456, 987654, 101010,
8675309, 81891, 222222, 999999, 110011)
experiments1 <- foreach(
s = seeds) %do%
{
set.seed(s)
x <- garchSim(
garchSpec(
model = list(
alpha = 0.2, beta = 0.2, omega = 0.2)),
n.start = 1000, n = 1000)
params <- foreach(
t = 50:1000,
.combine = rbind,
.packages = c("fGarch")) %dopar%
{
getFitData(x[1:t])
}
rownames(params) <- 50:1000
params
}
names(experiments1) <- seeds
experiments1 <- lapply(experiments1, param_reshape)
names(experiments1) <- c(
103117, 123456, 987654, 101010,
8675309, 81891, 222222, 999999, 110011)
experiments1_df <- ldply(experiments1, .id = "seed")
head(experiments1_df)
## seed t parameter value se
## 1 103117 50 omega 0.1043139 0.9830089
## 2 103117 51 omega 0.1037479 4.8441246
## 3 103117 52 omega 0.1032197 4.6421147
## 4 103117 53 omega 0.1026722 1.3041128
## 5 103117 54 omega 0.1020266 0.5334988
## 6 103117 55 omega 0.2725939 0.6089607
ggplot(
experiments1_df,
aes(x = t, y = value)) +
geom_line() +
geom_ribbon(
aes(ymin = value - 2 * se,
ymax = value + 2 * se),
color = "grey", alpha = 0.5) +
geom_hline(yintercept = 0.2, color = "blue") +
scale_y_continuous(
breaks = c(0, 0.2, 0.25, 0.5, 0.75, 1)) +
coord_cartesian(ylim = c(0, 1)) +
facet_grid(seed ~ parameter) +
ggtitle(
"Successive parameter estimates using NLMINB optimization")
在這個圖中我們看到了 \(\beta\) 的其他類型的病態行為，特別是種子是 222222 和 999999 的情況，其中 \(\beta\) 長期遠低於正確值。對於所有這些模擬，\(\beta\) 開始比正確的值大得多，接近 1，對於前面提到的兩個種子，\(\beta\) 從非常高的水平突然跳到非常低的水平。(此處未顯示種子 110131 和 110137 的結果，它們甚至更糟！)
其他參數也存在自己的病態行為，但情況似乎並不那么嚴峻。我們看到的病態行為可能與 \(\beta\) 的估計有關。實際上，如果我們看一下 \(\text{ARCH}(1)\) 過程的類似實驗(這是一個 \(\text{GARCH}(1,0)\) 過程，相當於設置 \(\beta = 0\))，我們會看到更好的行為。
set.seed(110117)
x <- garchSim(
garchSpec(
model = list(alpha = 0.2, beta = 0.2, omega = 0.2)),
n.start = 1000, n = 1000)
xarch <- garchSim(
garchSpec(
model = list(omega = 0.2, alpha = 0.2, beta = 0)),
n.start = 1000, n = 1000)
params_arch <- foreach(
t = 50:1000,
.combine = rbind,
.packages = c("fGarch")) %dopar%
{
getFitData(
xarch[1:t], formula = ~ garch(1, 0))
}
rownames(params_arch) <- 50:1000
print(ggp %+% param_reshape(params_arch) + ggtitle("ARCH(1) Model"))
病態行為似乎是數值性的，並且與 \(\beta\) 密切相關。默認情況下，garchFit() 使用 nlminb()（帶約束的擬牛頓方法）來解決優化問題，使用數值計算出的梯度。不過，我們可以選擇其他方法。我們可以使用 L-BFGS-B 方法，以及 Nelder-Mead 方法。
不幸的是，這些替代優化算法沒有做得更好，他們甚至可能會做得更糟。
# lbfgsb algorithm
params_lbfgsb <- foreach(
t = 50:1000,
.combine = rbind,
.packages = c("fGarch")) %dopar%
{
getFitData(x[1:t], algorithm = "lbfgsb")
}
rownames(params_lbfgsb) <- 50:1000
# nlminb+nm algorithm
params_nlminbnm <- foreach(
t = 50:1000,
.combine = rbind,
.packages = c("fGarch")) %dopar%
{
getFitData(x[1:t], algorithm = "nlminb+nm")
}
rownames(params_nlminbnm) <- 50:1000
# lbfgsb+nm algorithm
params_lbfgsbnm <- foreach(
t = 50:1000,
.combine = rbind,
.packages = c("fGarch")) %dopar%
{
getFitData(x[1:t], algorithm = "lbfgsb+nm")
}
rownames(params_lbfgsbnm) <- 50:1000
# cond.dist is norm (default)
params_norm <- foreach(
t = 50:1000,
.combine = rbind,
.packages = c("fGarch")) %dopar%
{
getFitData(x[1:t], cond.dist = "norm")
}
rownames(params_norm) <- 50:1000
print(ggp %+% param_reshape(params_lbfgsb) + ggtitle("L-BFGS-B Optimization"))
print(ggp %+% param_reshape(params_nlminbnm) + ggtitle("nlminb Optimization with Nelder-Mead"))
print(ggp %+% param_reshape(params_lbfgsbnm) + ggtitle("L-BFGS-B Optimization with Nelder-Mead"))
誠然，QMLE 並非 garchFit() 的默認估計方法，默認的是正態分布。不幸的是，這沒有帶來更好的結果。
print(ggp %+% param_reshape(params_norm) + ggtitle("cond.dist = 'norm'"))
在 CRAN，fGarch 自 2013 年以來沒有看到更新！有可能 fGarch 開始顯現出它的落伍老邁，新的包裝已經解決了我在這里強調的一些問題。包 tseries 提供了一個函數 garch()，它也通過 QMLE 擬合 \(\text{GARCH}(1,1)\) 模型，並且比 fGarch 更新。它是我所知道的唯一可以擬合 \(\text{GARCH}(1,1)\) 模型的其他包。
不幸的是，garch() 沒有做得更好。事實上，它似乎更糟糕。問題再一次出現在 \(\beta\) 身上。
library(tseries)
getFitDatagarch <- function(x)
{
garch(x)$coef
}
params_tseries <- foreach(
t = 50:1000,
.combine = rbind,
.packages = c("tseries")) %dopar%
{
getFitDatagarch(x[1:t])
}
rownames(params_tseries) <- 50:1000
param_reshape_tseries <- function(p)
{
p <- as.data.frame(p)
p$t <- as.integer(rownames(p))
pnew <- melt(
p, id.vars = "t", variable.name = "parameter")
pnew$parameter <- as.character(pnew$parameter)
return(pnew)
}
ggplot(
param_reshape_tseries(params_tseries),
aes(x = t, y = value)) +
geom_line() +
geom_hline(
yintercept = 0.2, color = "blue") +
scale_y_continuous(
breaks = c(0, 0.2, 0.25, 0.5, 0.75, 1)) +
coord_cartesian(ylim = c(0, 1)) +
facet_grid(. ~ parameter)
所有這些實驗均在固定（但隨機選擇）的序列上進行。實驗顯示，對於樣本量小於 300（可能更大的數字）的情況，\(\text{GARCH}(1,1)\) 參數估計的分布是可疑的。當我們模擬許多過程並查看參數的分布時會發生什么？
我模擬了 10000 個樣本大小為 100、500 和 1000 的 \(\text{GARCH}(1,1)\) 過程（使用與之前相同的參數）。以下是參數估計的經驗分布。
experiments2 <- foreach(
n = c(100, 500, 1000)) %do%
{
mat <- foreach(
i = 1:10000,
.combine = rbind,
.packages = c("fGarch")) %dopar%
{
x <- garchSim(
garchSpec(
model = list(
omega = 0.2, alpha = 0.2, beta = 0.2)),
n.start = 1000,
n = n)
getFitData(x)
}
rownames(mat) <- NULL
mat
}
names(experiments2) <- c(100, 500, 1000)
save(params, x, experiments1,
xarch, params_arch, params_lbfgsb,
params_nlminbnm, params_lbfgsbnm,
params_norm, params_tseries,
experiments2,
file = "garchfitexperiments.Rda")
param_sim <- lapply(
experiments2,
function(mat)
{
df <- as.data.frame(mat)
df <- df[c("omega", "alpha1", "beta1")]
return(df)
}) %>%
ldply(.id = "n")
param_sim <- param_sim %>%
melt(id.vars = "n", variable.name = "parameter")
head(param_sim)
## n parameter value
## 1 100 omega 8.015968e-02
## 2 100 omega 2.493595e-01
## 3 100 omega 2.300699e-01
## 4 100 omega 3.674244e-07
## 5 100 omega 2.697577e-03
## 6 100 omega 2.071737e-01
ggplot(
param_sim,
aes(x = value)) +
geom_density(
fill = "grey", alpha = 0.7) +
geom_vline(
xintercept = 0.2, color = "blue") +
facet_grid(
n ~ parameter)
當樣本量為 100 時，這些估計遠非可靠。\(\omega\) 和 \(\alpha\) 以一種令人不安的傾向趨近於 0，而 \(\beta\) 幾乎可以說是任何東西。如上所述，garchFit() 報告的標准差不會捕獲這種行為。對於較大的樣本量，\(\omega\) 和 \(\alpha\) 表現得更好，但 \(\beta\) 仍顯示出令人不安的行為。它的變化幅度幾乎沒有變化，並且它仍然有過小的傾向。
最讓我困擾的是樣本量為 1000 讓我感覺很大。如果一個人正在查看股票價格的每日數據，那么此樣本大小大致相當於 4 年的數據。這告訴我，這種病態行為正在影響人們現在試圖估計並在模型中使用的 GARCH 模型。
結論
由 John C. Nash 撰寫的題為《On best practice optimization methods in R》的文章，發表於 2014 年 9 月的 Journal of Statistical Software，討論了 R 需要更好的優化計算實踐。特別是，他強調了 garchFit() 使用了過時的方法（或至少它們的 R 實現）。他主張在社區中提高對優化問題的認識，並提高包的靈活性，而不僅僅是使用 optim() 提供的不同算法。
我在本文中強調的問題讓我更加意識到選擇在優化方法中的重要性。我最初的目標是編寫一個函數，用於根據 GARCH 模型中的結構性變化執行統計檢驗。正如我在此演示的那樣，這些檢驗嚴重依賴於對模型參數的連續估計。至少我的實驗表明，參數的變化沒有被標准差充分捕獲，同時也存在參數估計中不可接受的高度不穩定性。它們是如此不穩定，它將使檢驗不拒絕無變化的原假設成為一項奇跡（譯注：幾乎必然拒絕原假設）。畢竟，只要看一下模擬數據的圖片就可以得出結論，結構性變化的備擇假設是正確的。因此，每當我嘗試對原假設被認為是真的數據進行檢驗時，檢驗就明確拒絕它，其中 \(p\)-值基本上為 0。
我聽說有人們正對 GARCH 模型中的結構性變化進行假設檢驗研究，所以如果我在這里寫到的數值不穩定性可以避免，我不會對此感到驚訝。這是一個我自認知之甚少的主題，如果 R 社區中的某個人已經觀察到了這種行為並且知道如何解決它，我希望他們會在評論或電子郵件中告訴我。我可以寫一個回復並展示如何使用 garchFit() 生成參數的穩定估計。也許關鍵在於函數 garchFitControl()。
我也考慮過根據我的測試編寫自己的優化程序。Nash 教授在他的論文中強調了針對特定問題的需求定制優化程序的重要性。我已經寫下了要優化的量，我有一個 \(\text{GARCH}(1,1)\) 的梯度和 Hessian 矩陣的公式。也許我們的檢驗所要求的連續優化可以使用先前迭代中的參數作為初始值，從而有助於防止優化計算找到離群的、局部最優而全局次優的解。
雖然這使得問題比我最初想找一個我們檢驗的例子更難。我現在正在計划檢測 GARCH 模型中的結構性變化，但是僅涉及使用線性回歸的示例（一個更易處理的問題）。但我希望聽到別人對我在這里寫的內容的意見。
sessionInfo()
## R version 3.3.3 (2017-03-06)
## Platform: i686-pc-linux-gnu (32-bit)
## Running under: Ubuntu 16.04.2 LTS
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_0.7.2 plyr_1.8.4 reshape2_1.4.2
## [4] fGarch_3010.82.1 fBasics_3011.87 timeSeries_3022.101.2
## [7] timeDate_3012.100 ggplot2_2.2.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.11 bindr_0.1 knitr_1.17 magrittr_1.5
## [5] munsell_0.4.3 colorspace_1.3-2 R6_2.2.0 rlang_0.1.2
## [9] stringr_1.2.0 tools_3.3.3 grid_3.3.3 gtable_0.2.0
## [13] htmltools_0.3.6 assertthat_0.1 yaml_2.1.14 lazyeval_0.2.0
## [17] rprojroot_1.2 digest_0.6.12 tibble_1.3.4 bindrcpp_0.2
## [21] glue_1.1.1 evaluate_0.10 rmarkdown_1.6 labeling_0.3
## [25] stringi_1.1.5 scales_0.4.1 backports_1.0.5 pkgconfig_2.0.1
譯后記
原作者所指出的 GARCH 模型參數估計的不穩定性不僅使其本人困惑，也同樣令我震驚。我之前從未懷疑或質疑過統計軟件的計算結果，甚至沒有考慮過這個問題。今后在處理其他統計模型的參數估計問題時，務必首先用模擬數據檢驗一下相關軟件的結果穩健性。
回到 GARCH 模型參數估計的話題，我猜測 \(\beta\) 的不穩定性可能來自以下原因：
- GARCH 序列的統計性質對 \(\alpha\) 和 \(\beta\) 不敏感，特別是 \(\beta\)；
- \(\omega\)、\(\alpha\)、\(\beta\) 以及長期方差之間存在一個硬性的等式約束，但是在優化計算中沒有體現出這種等式約束。
\(\omega\)、\(\alpha\)、\(\beta\) 以及長期方差之間存在的硬性等式約束也許提供了技術性的補救手段：先估計 \(\omega\) 和 \(\alpha\)，再用等式約束和經驗長期方差推算出 \(\beta\)。
延續原作者的思路，下面是后續研究工作和博客寫作的安排：
- 驗證上述補救手段；
- 測試
rugarch的數值穩定性； - 測試
rmgarch和MTS等數值穩定性； - 測試
pyflux、arch和statsmodels等 python 包的數值穩定性。
GARCH 模型參數估計的不穩定性也引出了另一個問題，對於不可觀測的波動率的建模，參數估計以及校准的結果都是值得懷疑的。所以，某些 SDE 參數的估計和校准的穩定性實驗應該提上日程。
Bollerslev 在他 1986 年的文章（《General autoregressive conditional heteroscedasticity》）中介紹了 GARCH 模型。 ↩︎
《GARCH Models: Structure, Statistical Inference and Financial Applications》，Christian Francq 和 Jean-Michel Zakoian 著。 ↩︎