在 R 中估計 GARCH 參數存在的問題（基於 rugarch 包）
目錄html
在 R 中估計 GARCH 參數存在的問題（基於 rugarch 包）
本文翻譯自《Problems in Estimating GARCH Parameters in R (Part 2; rugarch)》linux
原文連接：https://ntguardian.wordpress.com/2019/01/28/problems-estimating-garch-parameters-r-part-2-rugarch/git
導論
這是一篇本應早就寫完的博客文章。一年前我寫了一篇文章，關於在 R 中估計 GARCH(1, 1) 模型參數時遇到的問題。我記錄了參數估計的行爲（重點是 \(\beta\)），以及使用 fGarch 計算這些估計值時發現的病態行爲。我在 R 社區呼籲幫助，包括經過 R Finance 郵件列表發送個人博客文章。web
反饋沒有讓我感到失望。你能夠看到一些郵件列表反饋，而且一些來自 Reddit 的評論也頗有幫助，但我認爲我獲得的最佳反饋來自於我本身的電子郵件。算法
Dr. Brian G. Peterson 做爲 R 金融社區的一員，給我發送了一些發人深思的電子郵件。首先，他告訴我 fGarch 再也不是處理 GARCH 模型的首選方案。RMetrics 套件包（包括 fGarch）由 ETH Zürich 的 Diethelm Würtz 教授維護。他在 2016 年的車禍中喪生。markdown
Dr. Peterson 建議我研究另外兩個用於 GARCH 建模的現代軟件包，rugarch（適用於單變量 GARCH 模型）和 rmgarch（適用於多變量 GARCH 模型）。以前我沒有據說過這些包（我之因此知道 fGarch 的緣由是由於它在由 Shumway 和 Stoffer 編寫的時間序列教科書——Time Series Analysis and Its Applications with R Examples 中），因此我很是感謝這個建議。因爲我如今對單變量時間序列感興趣，因此我研究了 rugarch。該軟件包彷佛具備比 fGarch 更多的功能和函數，這能夠解釋爲何它彷佛更難以使用。然而，包的 vignette 頗有幫助，值得打印出來。session
Dr. Peterson 對我提出的應用也有一些有趣的評論。他認爲，日內數據應優於日間數據，而且模擬數據（包括模擬 GARCH 過程）具備在實際數據中看不到的特質。獲取日間數據的便利性（特別是亞洲金融危機期間的 USD/JPY，這是我正在研究的檢驗統計量的預期應用）激發了我對日間數據的興趣。不過，他的評論可能會讓我從新考慮這個應用。1（我也許應該試圖經過 EUR/USD 來檢測 2010 年歐元區金融危機。爲此，我能夠從 HistData.com 得到免費的日內數據。）可是，若是對於小樣本而言不能信任標準差的估計，咱們的檢驗統計量仍然會遇到麻煩，由於它涉及小樣本的參數估計。app
他還警告說，模擬數據表現出在實際數據中看不到的行爲。這多是真的，但模擬數據很重要，由於它能夠被認爲是統計學家的最佳情景。另外，生成模擬數據的過程的屬性是先驗已知的，包括生成參數的值，以及哪些假設（例如序列中是否存在結構變化）是真的。這容許對估計器和檢驗進行健全的檢查。這對現實世界來講是不可能的，由於咱們沒有所需的先驗知識。dom
Prof. André Portela Santos 要求我重複模擬，但使用 \(\alpha = 0.6\)，由於這些值比我選擇的 \(\alpha = \beta = 0.2\) 更常見。這是一個很好的建議，除了 \(\alpha = \beta = 0.2\) 以外，我還會在博文裏考慮此範圍內的參數。然而，個人模擬暗示當 \(\alpha = \beta = 0.2\) 時，估計算法彷佛想要接近較大的 \(\beta\)。我也很驚訝，由於個人導師給個人印象是，\(\alpha\) 或 \(\beta\) 大的 GARCH 過程更難以處理。最後，若是估計量嚴重有偏，咱們可能會看到大多數估計參數位於該範圍內，但這並不意味着「正確」值位於該範圍內。個人模擬顯示 fGarch 很難發現 \(\alpha = \beta = 0.2\)，即便這些參數是「真的」。Prof. Santos 的評論讓我想要作一個在真實世界中 GARCH 參數的估計是什麼樣子的元研究（metastudy）。(可能有也可能沒有，我沒有檢查過。若是有人知道，請分享。)ide
個人導師聯繫了另外一位 GARCH 模型的專家，並得到了一些反饋。據推測，\(\beta\) 的標準差很大，所以參數估計應該有很大的變更範圍。即便對於小樣本，個人一些模擬也認同這種行爲，但同時顯示出對 \(\beta = 0\) 和 \(\beta = 1\) 使人不舒服的偏向。正如我假設的那樣，這多是優化程序的結果。
所以，鑑於此反饋，我將進行更多的模擬實驗。我不會再研究 fGarch 或 tseries 了，我將專門研究 rugarch。我將探討包支持的不一樣優化程序。我不會像我在第一篇文章中那樣畫圖，這些圖只是爲了代表存在的問題及其嚴重性。相反，我將考察由不一樣優化程序生成的估計器的特性。
rugarch 簡介
如上所述，rugarch 是一個用於處理 GARCH 模型的軟件包，一個主要的用例顯然是估計模型的參數。在這裏，我將演示如何指定 GARCH 模型、模擬模型的數據以及估計參數。在此以後，咱們能夠深刻了解模擬研究。
library(rugarch)
## Loading required package: parallel ## ## Attaching package: 'rugarch' ## The following object is masked from 'package:stats': ## ## sigma
指定一個 \(\text{GARCH}(1, 1)\) 模型
要使用 GARCH 模型，咱們須要指定它。執行此操做的函數是 ugarchspec()。我認爲最重要的參數是 variance.model 和 mean.model。
variance.model 是一個命名列表，也許最感興趣的兩個元素是 model 和 garchOrder。model 是一個字符串，指定擬合哪一種類型的 GARCH 模型。包支持許多主要的 GARCH 模型（例如 EGARCH、IGARCH 等），對於「普通」GARCH 模型，要將其設置爲 sGARCH（或者只是忽略它，標準模型是默認的)。garchOrder 是模型中 ARCH 和 GARCH 部分的階數向量。
mean.model 容許擬合 ARMA-GARCH 模型，而且像 variance.model 同樣接受一個命名列表，最感興趣的參數是 armaOrder 和 include.mean。armaOrder 就像 garchOrder，它是一個指定 ARMA 模型階數的向量。include.mean 是一個布爾值，若是爲 true，則容許模型的 ARMA 部分具備非零均值。
在模擬過程時，咱們須要設置參數的值。這是經過 fixed.pars 參數完成的，該參數接受命名列表，列表的元素是數字。它們須要符合函數對於參數的約定。例如，若是咱們想設置 \(\text{GARCH}(1,1)\) 模型的參數，咱們列表元素的名稱應該是 alpha1 和 beta1。若是計劃是模擬一個模型，則應以這種方式設置模型中的每一個參數。
還有其餘有趣的參數，但我只關注這些，由於默認指定是 ARMA-GARCH 模型，ARMA 階數爲 \((1,1)\)，非零均值，而且 GARCH 模型的階數是 \((1, 1)\)。這不是我想要的普通 \(\text{GARCH}(1,1)\) 模型，因此我幾乎老是要修改它。
spec1 <- ugarchspec(
mean.model = list(
armaOrder = c(0,0), include.mean = FALSE),
fixed.pars = list(
"omega" = 0.2, "alpha1" = 0.2, "beta1" = 0.2))
spec2 <- ugarchspec(
mean.model = list(
armaOrder = c(0,0), include.mean = FALSE),
fixed.pars = list(
"omega" = 0.2, "alpha1" = 0.1, "beta1" = 0.7))
show(spec1)
## ## *---------------------------------* ## * GARCH Model Spec * ## *---------------------------------* ## ## Conditional Variance Dynamics ## ------------------------------------ ## GARCH Model : sGARCH(1,1) ## Variance Targeting : FALSE ## ## Conditional Mean Dynamics ## ------------------------------------ ## Mean Model : ARFIMA(0,0,0) ## Include Mean : FALSE ## GARCH-in-Mean : FALSE ## ## Conditional Distribution ## ------------------------------------ ## Distribution : norm ## Includes Skew : FALSE ## Includes Shape : FALSE ## Includes Lambda : FALSE
show(spec2)
## ## *---------------------------------* ## * GARCH Model Spec * ## *---------------------------------* ## ## Conditional Variance Dynamics ## ------------------------------------ ## GARCH Model : sGARCH(1,1) ## Variance Targeting : FALSE ## ## Conditional Mean Dynamics ## ------------------------------------ ## Mean Model : ARFIMA(0,0,0) ## Include Mean : FALSE ## GARCH-in-Mean : FALSE ## ## Conditional Distribution ## ------------------------------------ ## Distribution : norm ## Includes Skew : FALSE ## Includes Shape : FALSE ## Includes Lambda : FALSE
模擬一個 GARCH 過程
函數 ugarchpath() 模擬由 ugarchspec() 指定的 GARCH 模型。該函數首先須要由 ugarchspec() 建立的指定對象。參數 n.sim 和 n.start 分別指定過程的大小和預熱期的長度（分別默認爲 1000 和 0。我強烈建議將預熱期設置爲至少 500，但我設置爲 1000）。該函數建立的對象不只包含模擬序列，還包含殘差和 \(\sigma_t\)。
rseed 參數控制函數用於生成數據的隨機種子。請注意，此函數會有效地忽略 set.seed()，所以若是須要一致的結果，則須要設置此參數。
這些對象相應的 plot() 方法並不徹底透明。它能夠建立一些圖，當在命令行中對 uGARCHpath 對象調用 plot() 時，系統會提示用戶輸入與所需圖形對應的數字。這有時挺痛苦，因此不要忘記將所需的編號傳遞給 which 參數以免提示，設置 which = 2 將正好給出序列的圖。
old_par <- par()
par(mfrow = c(2, 2))
x_obj <- ugarchpath(
spec1, n.sim = 1000, n.start = 1000, rseed = 111217)
show(x_obj)
## ## *------------------------------------* ## * GARCH Model Path Simulation * ## *------------------------------------* ## Model: sGARCH ## Horizon: 1000 ## Simulations: 1 ## Seed Sigma2.Mean Sigma2.Min Sigma2.Max Series.Mean ## sim1 111217 0.332 0.251 0.915 0.000165 ## Mean(All) 0 0.332 0.251 0.915 0.000165 ## Unconditional NA 0.333 NA NA 0.000000 ## Series.Min Series.Max ## sim1 -1.76 1.62 ## Mean(All) -1.76 1.62 ## Unconditional NA NA
for (i in 1:4)
{
plot(x_obj, which = i)
}
par(old_par)
## Warning in par(old_par): graphical parameter "cin" cannot be set ## Warning in par(old_par): graphical parameter "cra" cannot be set ## Warning in par(old_par): graphical parameter "csi" cannot be set ## Warning in par(old_par): graphical parameter "cxy" cannot be set ## Warning in par(old_par): graphical parameter "din" cannot be set ## Warning in par(old_par): graphical parameter "page" cannot be set
# The actual series x1 <- x_obj@path$seriesSim plot.ts(x1)
擬合一個 \(\text{GARCH}(1,1)\) 模型
ugarchfit() 函數擬合 GARCH 模型。該函數須要指定和數據集。solver 參數接受一個字符串，說明要使用哪一個數值優化器來尋找參數估計值。函數的大多數參數管理數值優化器的接口。特別是，solver.control 能夠接受一個傳遞給優化器的參數列表。咱們稍後會更詳細地討論這個問題。
用於生成模擬數據的指定將不適用於 ugarchfit()，由於它包含其參數的固定值。在個人狀況下，我將須要建立第二個指定對象。
spec <- ugarchspec(
mean.model = list(armaOrder = c(0, 0), include.mean = FALSE))
fit <- ugarchfit(spec, data = x1)
show(fit)
## ## *---------------------------------* ## * GARCH Model Fit * ## *---------------------------------* ## ## Conditional Variance Dynamics ## ----------------------------------- ## GARCH Model : sGARCH(1,1) ## Mean Model : ARFIMA(0,0,0) ## Distribution : norm ## ## Optimal Parameters ## ------------------------------------ ## Estimate Std. Error t value Pr(>|t|) ## omega 0.000713 0.001258 0.56696 0.57074 ## alpha1 0.002905 0.003714 0.78206 0.43418 ## beta1 0.994744 0.000357 2786.08631 0.00000 ## ## Robust Standard Errors: ## Estimate Std. Error t value Pr(>|t|) ## omega 0.000713 0.001217 0.58597 0.55789 ## alpha1 0.002905 0.003661 0.79330 0.42760 ## beta1 0.994744 0.000137 7250.45186 0.00000 ## ## LogLikelihood : -860.486 ## ## Information Criteria ## ------------------------------------ ## ## Akaike 1.7270 ## Bayes 1.7417 ## Shibata 1.7270 ## Hannan-Quinn 1.7326 ## ## Weighted Ljung-Box Test on Standardized Residuals ## ------------------------------------ ## statistic p-value ## Lag[1] 3.998 0.04555 ## Lag[2*(p+q)+(p+q)-1][2] 4.507 0.05511 ## Lag[4*(p+q)+(p+q)-1][5] 9.108 0.01555 ## d.o.f=0 ## H0 : No serial correlation ## ## Weighted Ljung-Box Test on Standardized Squared Residuals ## ------------------------------------ ## statistic p-value ## Lag[1] 29.12 6.786e-08 ## Lag[2*(p+q)+(p+q)-1][5] 31.03 1.621e-08 ## Lag[4*(p+q)+(p+q)-1][9] 32.26 1.044e-07 ## d.o.f=2 ## ## Weighted ARCH LM Tests ## ------------------------------------ ## Statistic Shape Scale P-Value ## ARCH Lag[3] 1.422 0.500 2.000 0.2331 ## ARCH Lag[5] 2.407 1.440 1.667 0.3882 ## ARCH Lag[7] 2.627 2.315 1.543 0.5865 ## ## Nyblom stability test ## ------------------------------------ ## Joint Statistic: 0.9518 ## Individual Statistics: ## omega 0.3296 ## alpha1 0.2880 ## beta1 0.3195 ## ## Asymptotic Critical Values (10% 5% 1%) ## Joint Statistic: 0.846 1.01 1.35 ## Individual Statistic: 0.35 0.47 0.75 ## ## Sign Bias Test ## ------------------------------------ ## t-value prob sig ## Sign Bias 0.3946 6.933e-01 ## Negative Sign Bias 3.2332 1.264e-03 *** ## Positive Sign Bias 4.2142 2.734e-05 *** ## Joint Effect 28.2986 3.144e-06 *** ## ## ## Adjusted Pearson Goodness-of-Fit Test: ## ------------------------------------ ## group statistic p-value(g-1) ## 1 20 20.28 0.3779 ## 2 30 26.54 0.5965 ## 3 40 36.56 0.5817 ## 4 50 47.10 0.5505 ## ## ## Elapsed time : 2.60606
par(mfrow = c(3, 4))
for (i in 1:12)
{
plot(fit, which = i)
}
## ## please wait...calculating quantiles...
par(old_par)
## Warning in par(old_par): graphical parameter "cin" cannot be set ## Warning in par(old_par): graphical parameter "cra" cannot be set ## Warning in par(old_par): graphical parameter "csi" cannot be set ## Warning in par(old_par): graphical parameter "cxy" cannot be set ## Warning in par(old_par): graphical parameter "din" cannot be set ## Warning in par(old_par): graphical parameter "page" cannot be set
注意估計的參數和標準差？即便對於 1000 的樣本大小，估計也與「正確」數字相去甚遠，而且基於估計標準差的合理置信區間不包含正確的值。看起來我在上一篇文章中記錄的問題並無消失。
出於好奇，在 Prof. Santos 建議範圍的其餘指定會發生什麼？
x_obj <- ugarchpath(
spec2, n.start = 1000, rseed = 111317)
x2 <- x_obj@path$seriesSim
fit <- ugarchfit(spec, x2)
show(fit)
## ## *---------------------------------* ## * GARCH Model Fit * ## *---------------------------------* ## ## Conditional Variance Dynamics ## ----------------------------------- ## GARCH Model : sGARCH(1,1) ## Mean Model : ARFIMA(0,0,0) ## Distribution : norm ## ## Optimal Parameters ## ------------------------------------ ## Estimate Std. Error t value Pr(>|t|) ## omega 0.001076 0.002501 0.43025 0.66701 ## alpha1 0.001992 0.002948 0.67573 0.49921 ## beta1 0.997008 0.000472 2112.23364 0.00000 ## ## Robust Standard Errors: ## Estimate Std. Error t value Pr(>|t|) ## omega 0.001076 0.002957 0.36389 0.71594 ## alpha1 0.001992 0.003510 0.56767 0.57026 ## beta1 0.997008 0.000359 2777.24390 0.00000 ## ## LogLikelihood : -1375.951 ## ## Information Criteria ## ------------------------------------ ## ## Akaike 2.7579 ## Bayes 2.7726 ## Shibata 2.7579 ## Hannan-Quinn 2.7635 ## ## Weighted Ljung-Box Test on Standardized Residuals ## ------------------------------------ ## statistic p-value ## Lag[1] 0.9901 0.3197 ## Lag[2*(p+q)+(p+q)-1][2] 1.0274 0.4894 ## Lag[4*(p+q)+(p+q)-1][5] 3.4159 0.3363 ## d.o.f=0 ## H0 : No serial correlation ## ## Weighted Ljung-Box Test on Standardized Squared Residuals ## ------------------------------------ ## statistic p-value ## Lag[1] 3.768 0.05226 ## Lag[2*(p+q)+(p+q)-1][5] 4.986 0.15424 ## Lag[4*(p+q)+(p+q)-1][9] 7.473 0.16272 ## d.o.f=2 ## ## Weighted ARCH LM Tests ## ------------------------------------ ## Statistic Shape Scale P-Value ## ARCH Lag[3] 0.2232 0.500 2.000 0.6366 ## ARCH Lag[5] 0.4793 1.440 1.667 0.8897 ## ARCH Lag[7] 2.2303 2.315 1.543 0.6686 ## ## Nyblom stability test ## ------------------------------------ ## Joint Statistic: 0.3868 ## Individual Statistics: ## omega 0.2682 ## alpha1 0.2683 ## beta1 0.2669 ## ## Asymptotic Critical Values (10% 5% 1%) ## Joint Statistic: 0.846 1.01 1.35 ## Individual Statistic: 0.35 0.47 0.75 ## ## Sign Bias Test ## ------------------------------------ ## t-value prob sig ## Sign Bias 0.5793 0.5625 ## Negative Sign Bias 1.3358 0.1819 ## Positive Sign Bias 1.5552 0.1202 ## Joint Effect 5.3837 0.1458 ## ## ## Adjusted Pearson Goodness-of-Fit Test: ## ------------------------------------ ## group statistic p-value(g-1) ## 1 20 24.24 0.1871 ## 2 30 30.50 0.3894 ## 3 40 38.88 0.4753 ## 4 50 48.40 0.4974 ## ## ## Elapsed time : 2.841597
沒有更好。如今讓咱們看看當咱們使用不一樣的優化算法時會發生什麼。
rugarch 中的優化與參數估計
優化器的選擇
ugarchfit() 的默認參數很好地找到了我稱之爲模型 2 的適當參數（其中 \(\alpha = 0.1\) 和 \(\beta = 0.7\)），但不適用於模型 1（\(\alpha = \beta = 0.2\))。我想知道的是什麼時候一個求解器能擊敗另外一個求解器。
正如 Vivek Rao2 在 R-SIG-Finance 郵件列表中所說，「最佳」估計是最大化似然函數（或等效地，對數似然函數）的估計，在上一篇文章中我忽略了檢查對數似然函數值。在這裏，我將看到哪些優化程序致使最大對數似然。
下面是一個輔助函數，它簡化了擬合 GARCH 模型參數、提取對數似然、參數值和標準差的過程，同時容許將不一樣的值傳遞給 solver 和 solver.control。
evalSolverFit <- function(spec, data,
solver = "solnp",
solver.control = list())
{
# Calls ugarchfit(spec, data, solver, solver.control), and returns a vector
# containing the log likelihood, parameters, and parameter standard errors.
# Parameters are equivalent to those seen in ugarchfit(). If the solver fails
# to converge, NA will be returned
vec <- NA
tryCatch(
{
fit <- ugarchfit(
spec = spec,
data = data,
solver = solver,
solver.control = solver.control)
coef_se_names <- paste(
"se", names(fit@fit$coef), sep = ".")
se <- fit@fit$se.coef
names(se) <- coef_se_names
robust_coef_se_names <- paste(
"robust.se", names(fit@fit$coef), sep = ".")
robust.se <- fit@fit$robust.se.coef
names(robust.se) <- robust_coef_se_names
vec <- c(fit@fit$coef, se, robust.se)
vec["LLH"] <- fit@fit$LLH
},
error = function(w) { NA })
return(vec)
}
下面我列出將要考慮的全部優化方案。我只使用 solver.control，但可能有其餘參數能夠幫助數值優化算法，即 numderiv.control，它們做爲控制參數傳遞給負責標準差計算的數值算法。這利用了包含 numDeriv 的包，它執行數值微分。
solvers <- list(
# A list of lists where each sublist contains parameters to
# pass to a solver
list("solver" = "nlminb", "solver.control" = list()),
list("solver" = "solnp", "solver.control" = list()),
list("solver" = "lbfgs", "solver.control" = list()),
list("solver" = "gosolnp",
"solver.control" = list("n.restarts" = 100, "n.sim" = 100)),
list("solver" = "hybrid", "solver.control" = list()),
list("solver" = "nloptr", "solver.control" = list("solver" = 1)), # COBYLA
list("solver" = "nloptr", "solver.control" = list("solver" = 2)), # BOBYQA
list("solver" = "nloptr", "solver.control" = list("solver" = 3)), # PRAXIS
list("solver" = "nloptr",
"solver.control" = list("solver" = 4)), # NELDERMEAD
list("solver" = "nloptr", "solver.control" = list("solver" = 5)), # SBPLX
list("solver" = "nloptr",
"solver.control" = list("solver" = 6)), # AUGLAG+COBYLA
list("solver" = "nloptr",
"solver.control" = list("solver" = 7)), # AUGLAG+BOBYQA
list("solver" = "nloptr",
"solver.control" = list("solver" = 8)), # AUGLAG+PRAXIS
list("solver" = "nloptr",
"solver.control" = list("solver" = 9)), # AUGLAG+NELDERMEAD
list("solver" = "nloptr",
"solver.control" = list("solver" = 10)) # AUGLAG+SBPLX
)
tags <- c(
# Names for the above list
"nlminb",
"solnp",
"lbfgs",
"gosolnp",
"hybrid",
"nloptr+COBYLA",
"nloptr+BOBYQA",
"nloptr+PRAXIS",
"nloptr+NELDERMEAD",
"nloptr+SBPLX",
"nloptr+AUGLAG+COBYLA",
"nloptr+AUGLAG+BOBYQA",
"nloptr+AUGLAG+PRAXIS",
"nloptr+AUGLAG+NELDERMEAD",
"nloptr+AUGLAG+SBPLX"
)
names(solvers) <- tags
如今讓咱們進行優化計算選擇的交叉射擊（gauntlet），看看哪一個算法產生的估計在模型 1 生成的數據上達到最大的對數似然。遺憾的是，lbfgs 方法（Broyden-Fletcher-Goldfarb-Shanno 方法的低存儲版本）在這個序列上沒有收斂，因此我省略了它。
optMethodCompare <- function(data,
spec,
solvers)
{
# Runs all solvers in a list for a dataset
#
# Args:
# data: An object to pass to ugarchfit's data parameter containing the data
# to fit
# spec: A specification created by ugarchspec to pass to ugarchfit
# solvers: A list of lists containing strings of solvers and a list for
# solver.control
#
# Return:
# A matrix containing the result of the solvers (including parameters, se's,
# and LLH)
model_solutions <- lapply(
solvers,
function(s)
{
args <- s
args[["spec"]] <- spec
args[["data"]] <- data
res <- do.call(evalSolverFit, args = args)
return(res)
})
model_solutions <- do.call(
rbind, model_solutions)
return(model_solutions)
}
round(
optMethodCompare(
x1, spec, solvers[c(1:2, 4:15)]), digits = 4)
## omega alpha1 beta1 se.omega se.alpha1 se.beta1 robust.se.omega robust.se.alpha1 robust.se.beta1 LLH ## ------------------------- ------- ------- ------- --------- ---------- --------- ---------------- ----------------- ---------------- ---------- ## nlminb 0.2689 0.1774 0.0000 0.0787 0.0472 0.2447 0.0890 0.0352 0.2830 -849.6927 ## solnp 0.0007 0.0029 0.9947 0.0013 0.0037 0.0004 0.0012 0.0037 0.0001 -860.4860 ## gosolnp 0.2689 0.1774 0.0000 0.0787 0.0472 0.2446 0.0890 0.0352 0.2828 -849.6927 ## hybrid 0.0007 0.0029 0.9947 0.0013 0.0037 0.0004 0.0012 0.0037 0.0001 -860.4860 ## nloptr+COBYLA 0.0006 0.0899 0.9101 0.0039 0.0306 0.0370 0.0052 0.0527 0.0677 -871.5006 ## nloptr+BOBYQA 0.0003 0.0907 0.9093 0.0040 0.0298 0.0375 0.0057 0.0532 0.0718 -872.3436 ## nloptr+PRAXIS 0.2689 0.1774 0.0000 0.0786 0.0472 0.2444 0.0888 0.0352 0.2823 -849.6927 ## nloptr+NELDERMEAD 0.0010 0.0033 0.9935 0.0013 0.0039 0.0004 0.0013 0.0038 0.0001 -860.4845 ## nloptr+SBPLX 0.0010 0.1000 0.9000 0.0042 0.0324 0.0386 0.0055 0.0536 0.0680 -872.2736 ## nloptr+AUGLAG+COBYLA 0.0006 0.0899 0.9101 0.0039 0.0306 0.0370 0.0052 0.0527 0.0677 -871.5006 ## nloptr+AUGLAG+BOBYQA 0.0003 0.0907 0.9093 0.0040 0.0298 0.0375 0.0057 0.0532 0.0718 -872.3412 ## nloptr+AUGLAG+PRAXIS 0.1246 0.1232 0.4948 0.0620 0.0475 0.2225 0.0701 0.0439 0.2508 -851.0547 ## nloptr+AUGLAG+NELDERMEAD 0.2689 0.1774 0.0000 0.0786 0.0472 0.2445 0.0889 0.0352 0.2826 -849.6927 ## nloptr+AUGLAG+SBPLX 0.0010 0.1000 0.9000 0.0042 0.0324 0.0386 0.0055 0.0536 0.0680 -872.2736
根據最大似然準則，「最優」結果是由 gosolnp 實現的。結果有一個不幸的屬性——\(\beta \approx 0\)，這固然不是正確的，但至少 \(\beta\) 的標準差會建立一個包含 \(\beta\) 真值的置信區間。其中，個人首選估計是由 AUGLAG + PRAXIS 生成的，由於 \(\beta\) 彷佛是合理的，事實上估計都接近事實（至少在置信區間包含真值的意義上），但不幸的是，即便它們是最合理的，估計並無最大化對數似然。
若是咱們看一下模型 2，咱們會看到什麼？一樣，lbfgs 沒有收斂，因此我省略忽略了它。不幸的是，nlminb 也沒有收斂，所以也必須省略。
round(
optMethodCompare(
x2, spec, solvers[c(2, 4:15)]), digits = 4)
## omega alpha1 beta1 se.omega se.alpha1 se.beta1 robust.se.omega robust.se.alpha1 robust.se.beta1 LLH ## ------------------------- ------- ------- ------- --------- ---------- --------- ---------------- ----------------- ---------------- ---------- ## solnp 0.0011 0.0020 0.9970 0.0025 0.0029 0.0005 0.0030 0.0035 0.0004 -1375.951 ## gosolnp 0.0011 0.0020 0.9970 0.0025 0.0029 0.0005 0.0030 0.0035 0.0004 -1375.951 ## hybrid 0.0011 0.0020 0.9970 0.0025 0.0029 0.0005 0.0030 0.0035 0.0004 -1375.951 ## nloptr+COBYLA 0.0016 0.0888 0.9112 0.0175 0.0619 0.0790 0.0540 0.2167 0.2834 -1394.529 ## nloptr+BOBYQA 0.0010 0.0892 0.9108 0.0194 0.0659 0.0874 0.0710 0.2631 0.3572 -1395.310 ## nloptr+PRAXIS 0.5018 0.0739 0.3803 0.3178 0.0401 0.3637 0.2777 0.0341 0.3225 -1373.632 ## nloptr+NELDERMEAD 0.0028 0.0026 0.9944 0.0028 0.0031 0.0004 0.0031 0.0035 0.0001 -1375.976 ## nloptr+SBPLX 0.0029 0.1000 0.9000 0.0146 0.0475 0.0577 0.0275 0.1108 0.1408 -1395.807 ## nloptr+AUGLAG+COBYLA 0.0016 0.0888 0.9112 0.0175 0.0619 0.0790 0.0540 0.2167 0.2834 -1394.529 ## nloptr+AUGLAG+BOBYQA 0.0010 0.0892 0.9108 0.0194 0.0659 0.0874 0.0710 0.2631 0.3572 -1395.310 ## nloptr+AUGLAG+PRAXIS 0.5018 0.0739 0.3803 0.3178 0.0401 0.3637 0.2777 0.0341 0.3225 -1373.632 ## nloptr+AUGLAG+NELDERMEAD 0.0001 0.0000 1.0000 0.0003 0.0003 0.0000 0.0004 0.0004 0.0000 -1375.885 ## nloptr+AUGLAG+SBPLX 0.0029 0.1000 0.9000 0.0146 0.0475 0.0577 0.0275 0.1108 0.1408 -1395.807
這裏是 PRAXIS 和 AUGLAG + PRAXIS 給出了「最優」結果，只有這兩種方法作到了。其餘優化器給出了明顯糟糕的結果。也就是說，「最優」解在參數爲非零、置信區間包含正確值上是首選的。
若是咱們將樣本限制爲 100，會發生什麼？（lbfgs 仍然不起做用。）
round(
optMethodCompare(
x1[1:100], spec, solvers[c(1:2, 4:15)]), digits = 4)
## omega alpha1 beta1 se.omega se.alpha1 se.beta1 robust.se.omega robust.se.alpha1 robust.se.beta1 LLH ## ------------------------- ------- ------- ------- --------- ---------- --------- ---------------- ----------------- ---------------- --------- ## nlminb 0.0451 0.2742 0.5921 0.0280 0.1229 0.1296 0.0191 0.0905 0.0667 -80.6587 ## solnp 0.0451 0.2742 0.5921 0.0280 0.1229 0.1296 0.0191 0.0905 0.0667 -80.6587 ## gosolnp 0.0451 0.2742 0.5921 0.0280 0.1229 0.1296 0.0191 0.0905 0.0667 -80.6587 ## hybrid 0.0451 0.2742 0.5921 0.0280 0.1229 0.1296 0.0191 0.0905 0.0667 -80.6587 ## nloptr+COBYLA 0.0007 0.1202 0.8798 0.0085 0.0999 0.0983 0.0081 0.1875 0.1778 -85.3121 ## nloptr+BOBYQA 0.0005 0.1190 0.8810 0.0085 0.0994 0.0992 0.0084 0.1892 0.1831 -85.3717 ## nloptr+PRAXIS 0.0451 0.2742 0.5921 0.0280 0.1229 0.1296 0.0191 0.0905 0.0667 -80.6587 ## nloptr+NELDERMEAD 0.0451 0.2742 0.5920 0.0281 0.1230 0.1297 0.0191 0.0906 0.0667 -80.6587 ## nloptr+SBPLX 0.0433 0.2740 0.5998 0.0269 0.1237 0.1268 0.0182 0.0916 0.0648 -80.6616 ## nloptr+AUGLAG+COBYLA 0.0007 0.1202 0.8798 0.0085 0.0999 0.0983 0.0081 0.1875 0.1778 -85.3121 ## nloptr+AUGLAG+BOBYQA 0.0005 0.1190 0.8810 0.0085 0.0994 0.0992 0.0084 0.1892 0.1831 -85.3717 ## nloptr+AUGLAG+PRAXIS 0.0451 0.2742 0.5921 0.0280 0.1229 0.1296 0.0191 0.0905 0.0667 -80.6587 ## nloptr+AUGLAG+NELDERMEAD 0.0451 0.2742 0.5921 0.0280 0.1229 0.1296 0.0191 0.0905 0.0667 -80.6587 ## nloptr+AUGLAG+SBPLX 0.0450 0.2742 0.5924 0.0280 0.1230 0.1295 0.0191 0.0906 0.0666 -80.6587
round(
optMethodCompare(
x2[1:100], spec, solvers[c(1:2, 4:15)]), digits = 4)
## omega alpha1 beta1 se.omega se.alpha1 se.beta1 robust.se.omega robust.se.alpha1 robust.se.beta1 LLH ## ------------------------- ------- ------- ------- --------- ---------- --------- ---------------- ----------------- ---------------- ---------- ## nlminb 0.7592 0.0850 0.0000 2.1366 0.4813 3.0945 7.5439 1.7763 11.0570 -132.4614 ## solnp 0.0008 0.0000 0.9990 0.0291 0.0417 0.0066 0.0232 0.0328 0.0034 -132.9182 ## gosolnp 0.0537 0.0000 0.9369 0.0521 0.0087 0.0713 0.0430 0.0012 0.0529 -132.9124 ## hybrid 0.0008 0.0000 0.9990 0.0291 0.0417 0.0066 0.0232 0.0328 0.0034 -132.9182 ## nloptr+COBYLA 0.0014 0.0899 0.9101 0.0259 0.0330 0.1192 0.0709 0.0943 0.1344 -135.7495 ## nloptr+BOBYQA 0.0008 0.0905 0.9095 0.0220 0.0051 0.1145 0.0687 0.0907 0.1261 -135.8228 ## nloptr+PRAXIS 0.0602 0.0000 0.9293 0.0522 0.0088 0.0773 0.0462 0.0015 0.0565 -132.9125 ## nloptr+NELDERMEAD 0.0024 0.0000 0.9971 0.0473 0.0629 0.0116 0.0499 0.0680 0.0066 -132.9186 ## nloptr+SBPLX 0.0027 0.1000 0.9000 0.0238 0.0493 0.1308 0.0769 0.1049 0.1535 -135.9175 ## nloptr+AUGLAG+COBYLA 0.0014 0.0899 0.9101 0.0259 0.0330 0.1192 0.0709 0.0943 0.1344 -135.7495 ## nloptr+AUGLAG+BOBYQA 0.0008 0.0905 0.9095 0.0221 0.0053 0.1145 0.0687 0.0907 0.1262 -135.8226 ## nloptr+AUGLAG+PRAXIS 0.0602 0.0000 0.9294 0.0523 0.0090 0.0771 0.0462 0.0014 0.0565 -132.9125 ## nloptr+AUGLAG+NELDERMEAD 0.0000 0.0000 0.9999 0.0027 0.0006 0.0005 0.0013 0.0004 0.0003 -132.9180 ## nloptr+AUGLAG+SBPLX 0.0027 0.1000 0.9000 0.0238 0.0493 0.1308 0.0769 0.1049 0.1535 -135.9175
結果並不使人興奮。多個求解器得到了模型 1 生成序列的「最佳」結果，同時 \(\omega\) 的 95％ 置信區間（CI）不包含 \(\omega\) 的真實值，儘管其餘的 CI 將包含其真實值。對於由模型 2 生成的序列，最佳結果是由 nlminb 求解器實現的，但參數值不合理，標準差很大。至少 CI 將包含正確值。
從這裏開始，咱們不該再僅僅關注兩個序列，而是在兩個模型生成的許多模擬序列中研究這些方法的表現。這篇文章中的模擬對於個人筆記本電腦而言計算量太大，所以我將使用我院系的超級計算機來執行它們，利用其多核進行並行計算。
library(foreach)
library(doParallel)
logfile <- ""
# logfile <- "outfile.log"
# if (!file.exists(logfile)) {
# file.create(logfile)
# }
cl <- makeCluster(
detectCores() - 1, outfile = logfile)
registerDoParallel(cl)
optMethodSims <- function(
gen_spec,
n.sim = 1000,
m.sim = 1000,
fit_spec = ugarchspec(
mean.model = list(
armaOrder = c(0,0), include.mean = FALSE)),
solvers = list(
"solnp" = list(
"solver" = "solnp", "solver.control" = list())),
rseed = NA, verbose = FALSE)
{
# Performs simulations in parallel of GARCH processes via rugarch and returns
# a list with the results of different optimization routines
#
# Args:
# gen_spec: The specification for generating a GARCH sequence, produced by
# ugarchspec
# n.sim: An integer denoting the length of the simulated series
# m.sim: An integer for the number of simulated sequences to generate
# fit_spec: A ugarchspec specification for the model to fit
# solvers: A list of lists containing strings of solvers and a list for
# solver.control
# rseed: Optional seeding value(s) for the random number generator. For
# m.sim>1, it is possible to provide either a single seed to
# initialize all values, or one seed per separate simulation (i.e.
# m.sim seeds). However, in the latter case this may result in some
# slight overhead depending on how large m.sim is
# verbose: Boolean for whether to write data tracking the progress of the
# loop into an output file
# outfile: A string for the file to store verbose output to (relevant only
# if verbose is TRUE)
#
# Return:
# A list containing the result of calling optMethodCompare on each generated
# sequence
fits <- foreach(
i = 1:m.sim,
.packages = c("rugarch"),
.export = c(
"optMethodCompare", "evalSolverFit")) %dopar%
{
if (is.na(rseed))
{
newseed <- NA
}
else if (is.vector(rseed))
{
newseed <- rseed[i]
}
else
{
newseed <- rseed + i - 1
}
if (verbose)
{
cat(as.character(Sys.time()), ": Now on simulation ", i, "\n")
}
sim <- ugarchpath(
gen_spec, n.sim = n.sim, n.start = 1000,
m.sim = 1, rseed = newseed)
x <- sim@path$seriesSim
optMethodCompare(
x, spec = fit_spec, solvers = solvers)
}
return(fits)
}
# Specification 1 first
spec1_n100 <- optMethodSims(
spec1, n.sim = 100, m.sim = 1000,
solvers = solvers, verbose = TRUE)
spec1_n500 <- optMethodSims(
spec1, n.sim = 500, m.sim = 1000,
solvers = solvers, verbose = TRUE)
spec1_n1000 <- optMethodSims(
spec1, n.sim = 1000, m.sim = 1000,
solvers = solvers, verbose = TRUE)
# Specification 2 next
spec2_n100 <- optMethodSims(
spec2, n.sim = 100, m.sim = 1000,
solvers = solvers, verbose = TRUE)
spec2_n500 <- optMethodSims(
spec2, n.sim = 500, m.sim = 1000,
solvers = solvers, verbose = TRUE)
spec2_n1000 <- optMethodSims(
spec2, n.sim = 1000, m.sim = 1000,
solvers = solvers, verbose = TRUE)
如下是一組輔助函數，用於我要進行的分析。
optMethodSims_getAllVals <- function(param,
solver,
reslist)
{
# Get all values for a parameter obtained by a certain solver after getting a
# list of results via optMethodSims
#
# Args:
# param: A string for the parameter to get (such as "beta1")
# solver: A string for the solver for which to get the parameter (such as
# "nlminb")
# reslist: A list created by optMethodSims
#
# Return:
# A vector of values of the parameter for each simulation
res <- sapply(
reslist,
function(l)
{
return(l[solver, param])
})
return(res)
}
optMethodSims_getBestVals <- function(reslist,
opt_vec = TRUE,
reslike = FALSE)
{
# A function that gets the optimizer that maximized the likelihood function
# for each entry in reslist
#
# Args:
# reslist: A list created by optMethodSims
# opt_vec: A boolean indicating whether to return a vector with the name of
# the optimizers that maximized the log likelihood
# reslike: A bookean indicating whether the resulting list should consist of
# matrices of only one row labeled "best" with a structure like
# reslist
#
# Return:
# If opt_vec is TRUE, a list of lists, where each sublist contains a vector
# of strings naming the opimizers that maximized the likelihood function and
# a matrix of the parameters found. Otherwise, just the matrix (resembles
# the list generated by optMethodSims)
res <- lapply(
reslist,
function(l)
{
max_llh <- max(l[, "LLH"], na.rm = TRUE)
best_idx <- (l[, "LLH"] == max_llh) & (!is.na(l[, "LLH"]))
best_mat <- l[best_idx, , drop = FALSE]
if (opt_vec)
{
return(
list(
"solvers" = rownames(best_mat), "params" = best_mat))
}
else
{
return(best_mat)
}
})
if (reslike)
{
res <- lapply(
res,
function(l)
{
mat <- l$params[1, , drop = FALSE]
rownames(mat) <- "best"
return(mat)
})
}
return(res)
}
optMethodSims_getCaptureRate <- function(param,
solver,
reslist,
multiplier = 2,
spec,
use_robust = TRUE)
{
# Gets the rate a confidence interval for a parameter captures the true value
#
# Args:
# param: A string for the parameter being worked with
# solver: A string for the solver used to estimate the parameter
# reslist: A list created by optMethodSims
# multiplier: A floating-point number for the multiplier to the standard
# error, appropriate for the desired confidence level
# spec: A ugarchspec specification with the fixed parameters containing the
# true parameter value
# use_robust: Use robust standard errors for computing CIs
#
# Return:
# A float for the proportion of times the confidence interval captured the
# true parameter value
se_string <- ifelse(
use_robust, "robust.se.", "se.")
est <- optMethodSims_getAllVals(
param, solver, reslist)
moe_est <- multiplier * optMethodSims_getAllVals(
paste0(se_string, param), solver, reslist)
param_val <- spec@model$fixed.pars[[param]]
contained <- (param_val <= est + moe_est) & (param_val >= est - moe_est)
return(mean(contained, na.rm = TRUE))
}
optMethodSims_getMaxRate <- function(solver,
maxlist)
{
# Gets how frequently a solver found a maximal log likelihood
#
# Args:
# solver: A string for the solver
# maxlist A list created by optMethodSims_getBestVals with entries
# containing vectors naming the solvers that maximized the log
# likelihood
#
# Return:
# The proportion of times the solver maximized the log likelihood
maxed <- sapply(
maxlist,
function(l)
{
solver %in% l$solvers
})
return(mean(maxed))
}
optMethodSims_getFailureRate <- function(solver,
reslist)
{
# Computes the proportion of times a solver failed to converge.
#
# Args:
# solver: A string for the solver
# reslist: A list created by optMethodSims
#
# Return:
# Numeric proportion of times a solver failed to converge
failed <- sapply(
reslist,
function(l)
{
is.na(l[solver, "LLH"])
})
return(mean(failed))
}
# Vectorization
optMethodSims_getCaptureRate <- Vectorize(
optMethodSims_getCaptureRate, vectorize.args = "solver")
optMethodSims_getMaxRate <- Vectorize(
optMethodSims_getMaxRate, vectorize.args = "solver")
optMethodSims_getFailureRate <- Vectorize(
optMethodSims_getFailureRate, vectorize.args = "solver")
我首先爲固定樣本量和模型建立表：
- 全部求解器中，某個求解器達到最高對數似然的頻率
- 某個求解器未能收斂的頻率
- 基於某個求解器的解，95％ 置信區間包含每一個參數真實值的頻率（稱爲「捕獲率」，並使用穩健標準差）
solver_table <- function(reslist,
tags,
spec)
{
# Creates a table describing important solver statistics
#
# Args:
# reslist: A list created by optMethodSims
# tags: A vector with strings naming all solvers to include in the table
# spec: A ugarchspec specification with the fixed parameters containing the
# true parameter value
#
# Return:
# A matrix containing metrics describing the performance of the solvers
params <- names(spec1@model$fixed.pars)
max_rate <- optMethodSims_getMaxRate(
tags, optMethodSims_getBestVals(reslist))
failure_rate <- optMethodSims_getFailureRate(
tags, reslist)
capture_rate <- lapply(
params,
function(p)
{
optMethodSims_getCaptureRate(
p, tags, reslist, spec = spec)
})
return_mat <- cbind(
"Maximization Rate" = max_rate,
"Failure Rate" = failure_rate)
capture_mat <- do.call(cbind, capture_rate)
colnames(capture_mat) <- paste(
params, "95% CI Capture Rate")
return_mat <- cbind(
return_mat, capture_mat)
return(return_mat)
}
- Model 1, \(n=100\)
as.data.frame(
round(
solver_table(spec1_n100, tags, spec1) * 100,
digits = 1))
## Maximization Rate Failure Rate omega 95% CI Capture Rate alpha1 95% CI Capture Rate beta1 95% CI Capture Rate ## ------------------------- ------------------ ------------- -------------------------- --------------------------- -------------------------- ## nlminb 16.2 20.0 21.8 29.2 24.0 ## solnp 0.1 0.0 13.7 24.0 15.4 ## lbfgs 15.1 35.2 56.6 67.9 58.0 ## gosolnp 20.3 0.0 20.3 32.6 21.9 ## hybrid 0.1 0.0 13.7 24.0 15.4 ## nloptr+COBYLA 0.0 0.0 6.3 82.6 19.8 ## nloptr+BOBYQA 0.0 0.0 5.4 82.1 18.5 ## nloptr+PRAXIS 15.8 0.0 42.1 54.5 44.1 ## nloptr+NELDERMEAD 0.4 0.0 5.7 19.3 8.1 ## nloptr+SBPLX 0.1 0.0 7.7 85.7 24.1 ## nloptr+AUGLAG+COBYLA 0.0 0.0 6.1 84.5 19.9 ## nloptr+AUGLAG+BOBYQA 0.1 0.0 6.5 83.2 19.4 ## nloptr+AUGLAG+PRAXIS 22.6 0.0 41.2 54.6 44.1 ## nloptr+AUGLAG+NELDERMEAD 11.1 0.0 7.5 18.8 9.7 ## nloptr+AUGLAG+SBPLX 0.6 0.0 7.9 86.5 23.0
- Model 1, \(n=500\)
as.data.frame(
round(
solver_table(spec1_n500, tags, spec1) * 100,
digits = 1))
## Maximization Rate Failure Rate omega 95% CI Capture Rate alpha1 95% CI Capture Rate beta1 95% CI Capture Rate ## ------------------------- ------------------ ------------- -------------------------- --------------------------- -------------------------- ## nlminb 21.2 0.4 63.3 67.2 63.8 ## solnp 0.1 0.2 32.2 35.6 32.7 ## lbfgs 4.5 41.3 85.0 87.6 85.7 ## gosolnp 35.1 0.0 69.0 73.2 69.5 ## hybrid 0.1 0.0 32.3 35.7 32.8 ## nloptr+COBYLA 0.0 0.0 3.2 83.3 17.8 ## nloptr+BOBYQA 0.0 0.0 3.5 81.5 18.1 ## nloptr+PRAXIS 18.0 0.0 83.9 87.0 84.2 ## nloptr+NELDERMEAD 0.0 0.0 16.4 20.7 16.7 ## nloptr+SBPLX 0.1 0.0 3.7 91.4 15.7 ## nloptr+AUGLAG+COBYLA 0.0 0.0 3.2 83.3 17.8 ## nloptr+AUGLAG+BOBYQA 0.0 0.0 3.5 81.5 18.1 ## nloptr+AUGLAG+PRAXIS 21.9 0.0 80.2 87.4 83.4 ## nloptr+AUGLAG+NELDERMEAD 0.6 0.0 20.0 24.0 20.5 ## nloptr+AUGLAG+SBPLX 0.0 0.0 3.7 91.4 15.7
- Model 1, \(n=1000\)
as.data.frame(
round(
solver_table(spec1_n1000, tags, spec1) * 100,
digits = 1))
## Maximization Rate Failure Rate omega 95% CI Capture Rate alpha1 95% CI Capture Rate beta1 95% CI Capture Rate ## ------------------------- ------------------ ------------- -------------------------- --------------------------- -------------------------- ## nlminb 21.5 0.1 88.2 86.1 87.8 ## solnp 0.4 0.2 54.9 53.6 54.6 ## lbfgs 1.1 44.8 91.5 88.0 91.8 ## gosolnp 46.8 0.0 87.2 85.1 87.0 ## hybrid 0.5 0.0 55.0 53.6 54.7 ## nloptr+COBYLA 0.0 0.0 4.1 74.5 15.0 ## nloptr+BOBYQA 0.0 0.0 3.6 74.3 15.9 ## nloptr+PRAXIS 17.7 0.0 92.6 90.2 92.2 ## nloptr+NELDERMEAD 0.0 0.0 30.5 29.6 30.9 ## nloptr+SBPLX 0.0 0.0 3.0 82.3 11.6 ## nloptr+AUGLAG+COBYLA 0.0 0.0 4.1 74.5 15.0 ## nloptr+AUGLAG+BOBYQA 0.0 0.0 3.6 74.3 15.9 ## nloptr+AUGLAG+PRAXIS 13.0 0.0 83.4 93.9 86.7 ## nloptr+AUGLAG+NELDERMEAD 0.0 0.0 34.6 33.8 35.0 ## nloptr+AUGLAG+SBPLX 0.0 0.0 3.0 82.3 11.6
- Model 2, \(n=100\)
as.data.frame(
round(
solver_table(spec2_n100, tags, spec2) * 100,
digits = 1))
## Maximization Rate Failure Rate omega 95% CI Capture Rate alpha1 95% CI Capture Rate beta1 95% CI Capture Rate ## ------------------------- ------------------ ------------- -------------------------- --------------------------- -------------------------- ## nlminb 8.2 24.2 22.3 34.7 23.9 ## solnp 0.3 0.0 21.1 32.6 21.3 ## lbfgs 11.6 29.5 74.9 73.2 70.4 ## gosolnp 19.0 0.0 31.9 41.2 30.8 ## hybrid 0.3 0.0 21.1 32.6 21.3 ## nloptr+COBYLA 0.0 0.0 20.5 94.7 61.7 ## nloptr+BOBYQA 0.2 0.0 19.3 95.8 62.2 ## nloptr+PRAXIS 16.0 0.0 70.2 57.2 52.8 ## nloptr+NELDERMEAD 0.2 0.0 7.8 27.8 14.1 ## nloptr+SBPLX 0.1 0.0 24.9 91.0 65.0 ## nloptr+AUGLAG+COBYLA 0.0 0.0 21.2 95.1 62.5 ## nloptr+AUGLAG+BOBYQA 0.9 0.0 20.1 96.2 62.5 ## nloptr+AUGLAG+PRAXIS 38.8 0.0 70.4 57.2 52.7 ## nloptr+AUGLAG+NELDERMEAD 14.4 0.0 10.7 26.0 16.1 ## nloptr+AUGLAG+SBPLX 0.1 0.0 25.8 91.9 65.5
- Model 2, \(n=500\)
as.data.frame(
round(
solver_table(spec2_n500, tags, spec2) * 100,
digits = 1))
## Maximization Rate Failure Rate omega 95% CI Capture Rate alpha1 95% CI Capture Rate beta1 95% CI Capture Rate ## ------------------------- ------------------ ------------- -------------------------- --------------------------- -------------------------- ## nlminb 1.7 1.6 35.0 37.2 34.2 ## solnp 0.1 0.2 46.2 48.6 45.3 ## lbfgs 2.2 38.4 85.2 88.1 82.3 ## gosolnp 5.2 0.0 74.9 77.8 72.7 ## hybrid 0.1 0.0 46.1 48.5 45.2 ## nloptr+COBYLA 0.0 0.0 8.2 100.0 40.5 ## nloptr+BOBYQA 0.0 0.0 9.5 100.0 41.0 ## nloptr+PRAXIS 17.0 0.0 83.8 85.1 81.0 ## nloptr+NELDERMEAD 0.0 0.0 26.9 38.2 27.0 ## nloptr+SBPLX 0.0 0.0 8.2 100.0 40.2 ## nloptr+AUGLAG+COBYLA 0.0 0.0 8.2 100.0 40.5 ## nloptr+AUGLAG+BOBYQA 0.0 0.0 9.5 100.0 41.0 ## nloptr+AUGLAG+PRAXIS 77.8 0.0 84.4 85.4 81.3 ## nloptr+AUGLAG+NELDERMEAD 1.1 0.0 32.5 40.3 32.3 ## nloptr+AUGLAG+SBPLX 0.0 0.0 8.2 100.0 40.2
- Model 2, \(n=1000\)
as.data.frame(
round(
solver_table(spec2_n1000, tags, spec2) * 100,
digits = 1))
## Maximization Rate Failure Rate omega 95% CI Capture Rate alpha1 95% CI Capture Rate beta1 95% CI Capture Rate ## ------------------------- ------------------ ------------- -------------------------- --------------------------- -------------------------- ## nlminb 2.7 0.7 64.1 68.0 63.8 ## solnp 0.0 0.0 70.1 73.8 69.8 ## lbfgs 0.0 43.4 90.6 91.5 89.9 ## gosolnp 3.2 0.0 87.5 90.3 86.9 ## hybrid 0.0 0.0 70.1 73.8 69.8 ## nloptr+COBYLA 0.0 0.0 2.3 100.0 20.6 ## nloptr+BOBYQA 0.0 0.0 2.5 100.0 22.6 ## nloptr+PRAXIS 14.1 0.0 89.1 91.3 88.5 ## nloptr+NELDERMEAD 0.0 0.0 46.3 55.6 45.4 ## nloptr+SBPLX 0.0 0.0 2.2 100.0 19.5 ## nloptr+AUGLAG+COBYLA 0.0 0.0 2.3 100.0 20.6 ## nloptr+AUGLAG+BOBYQA 0.0 0.0 2.5 100.0 22.6 ## nloptr+AUGLAG+PRAXIS 85.5 0.0 89.1 91.3 88.5 ## nloptr+AUGLAG+NELDERMEAD 0.3 0.0 51.9 58.2 51.3 ## nloptr+AUGLAG+SBPLX 0.0 0.0 2.2 100.0 19.5
這些表已經揭示了不少信息。通常來講，NLOpt 中提供的 AUGLAG-PRAXIS 方法（使用主軸求解器的增廣拉格朗日方法）彷佛對模型 2 最有效，特別是對於大樣本；而對於模型 1，gosolnp 方法（使用 Yinyu Ye 的 solnp 求解器，但使用隨機初始化和重啓）彷佛在大樣本上勝出。
然而，更大的故事是任何方法都不能成爲「最佳」，特別是在樣本量較小的狀況下。有些優化器始終未能達到最大對數似然，沒有優化器可以始終如一地得到最佳結果。此外，不一樣的優化器彷佛在不一樣的模型下表現更好。真實世界的數據——真正的模型參數從未被知道——暗示了要嘗試每一個優化器（或至少那些有可能最大化對數似然的），並選擇產生最大對數似然的結果。沒有算法值得信賴，都沒法成爲首選算法。
如今讓咱們看一下參數估計的分佈圖。首先是輔助函數。
library(ggplot2)
solver_density_plot <- function(param,
tags,
list_reslist,
sample_sizes,
spec)
{
# Given a parameter, creates a density plot for each solver's distribution
# at different sample sizes
#
# Args:
# param: A string for the parameter to plot
# tags: A character vector containing the solver names
# list_reslist: A list of lists created by optMethodSimsf, one for each
# sample size
# sample_sizes: A numeric vector identifying the sample size corresponding
# to each object in the above list
# spec: A ugarchspec object containing the specification that generated the
# datasets
#
# Returns:
# A ggplot object containing the plot generated
p <- spec@model$fixed.pars[[param]]
nlist <- lapply(
list_reslist,
function(l)
{
optlist <- lapply(
tags,
function(t)
{
return(
na.omit(
optMethodSims_getAllVals(param, t, l)))
})
names(optlist) <- tags
df <- stack(optlist)
names(df) <- c("param", "optimizer")
return(df)
})
ndf <- do.call(rbind, nlist)
ndf$n <- rep(
sample_sizes, times = sapply(nlist, nrow))
ggplot(
ndf, aes(x = param)) +
geom_density(
fill = "black", alpha = 0.5) +
geom_vline(
xintercept = p, color = "blue") +
facet_grid(
optimizer ~ n, scales = "free_y")
}
開始畫圖。
- \(\omega\) 估計，model 1
solver_density_plot(
"omega", tags,
list(spec1_n100, spec1_n500, spec1_n1000),
c(100, 500, 1000), spec1)
- \(\alpha\) 估計，model 1
solver_density_plot(
"alpha1", tags,
list(spec1_n100, spec1_n500, spec1_n1000),
c(100, 500, 1000), spec1)
- \(\beta\) 估計，model 1
solver_density_plot(
"beta1", tags,
list(spec1_n100, spec1_n500, spec1_n1000),
c(100, 500, 1000), spec1)
請記住，只有 1000 個模擬序列，而且優化器爲每一個序列生成解，所以原則上優化器結果不該該是獨立的，但優化器運行得很是糟糕的時候，這些密度圖纔看起來是相同的。但即便優化器的表現不是很糟糕（就像 gosolnp、PRAXIS 和 AUGLAG-PRAXIS 方法的狀況同樣)，有證據代表估計 \(\omega\) 和 \(\alpha\) 的估計錯誤地接近 0，而且 \(\beta\) 的估計錯誤地接近 1。對於較小的樣本，這些錯誤更明顯。也就是說，對於更好的優化器，估計應該看起來幾乎無偏，特別是對於 \(\omega\) 和 \(\alpha\)，可是即便對於大樣本，它們的變更範圍也很大，特別是對於 \(\beta\) 的估計。可是，對於 AUGLAG-PRAXIS 優化器來講狀況並不是如此，它彷佛產生有偏的估計。
讓咱們看看模型 2 的圖。
- \(\omega\) 估計，model 2
solver_density_plot(
"omega", tags,
list(spec2_n100, spec2_n500, spec2_n1000),
c(100, 500, 1000), spec2)
- alpha$ 估計，model 2
solver_density_plot(
"alpha1", tags,
list(spec2_n100, spec2_n500, spec2_n1000),
c(100, 500, 1000), spec2)
- \(\beta\) 估計，model 2
solver_density_plot(
"beta1", tags,
list(spec2_n100, spec2_n500, spec2_n1000),
c(100, 500, 1000), spec2)
對於模型 2 來講估計器並無那麼費力，可是圖顯示仍然不太樂觀。PRAXIS 和 AUGLAG-PRAXIS 方法彷佛表現良好，但對於小樣本量來講遠非吸引人。
到目前爲止，個人實驗代表，從業者不該該依賴任何一個優化器，而是應該嘗試不一樣的優化器，並選擇具備最大對數似然的結果。假設咱們將此優化算法稱爲「最佳」優化器。這個優化器執行效果如何？
咱們來看看。
- Model 1，\(n=100\)
as.data.frame(
round(
solver_table(
optMethodSims_getBestVals(
spec1_n100, reslike = TRUE),
"best", spec1) * 100,
digits = 1))
## Maximization Rate Failure Rate omega 95% CI Capture Rate alpha1 95% CI Capture Rate beta1 95% CI Capture Rate ## ----- ------------------ ------------- -------------------------- --------------------------- -------------------------- ## best 100 0 49.5 63.3 52.2
- Model 1，\(n=500\)
as.data.frame(
round(
solver_table(
optMethodSims_getBestVals(
spec1_n500, reslike = TRUE),
"best", spec1) * 100,
digits = 1))
## Maximization Rate Failure Rate omega 95% CI Capture Rate alpha1 95% CI Capture Rate beta1 95% CI Capture Rate ## ----- ------------------ ------------- -------------------------- --------------------------- -------------------------- ## best 100 0 86 88.8 86.2
- Model 1，\(n=1000\)
as.data.frame(
round(
solver_table(
optMethodSims_getBestVals(
spec1_n1000, reslike = TRUE),
"best", spec1) * 100,
digits = 1))
## Maximization Rate Failure Rate omega 95% CI Capture Rate alpha1 95% CI Capture Rate beta1 95% CI Capture Rate ## ----- ------------------ ------------- -------------------------- --------------------------- -------------------------- ## best 100 0 92.8 90.3 92.4
- Model 2，\(n=100\)
as.data.frame(
round(
solver_table(
optMethodSims_getBestVals(
spec2_n100, reslike = TRUE),
"best", spec2) * 100,
digits = 1))
## Maximization Rate Failure Rate omega 95% CI Capture Rate alpha1 95% CI Capture Rate beta1 95% CI Capture Rate ## ----- ------------------ ------------- -------------------------- --------------------------- -------------------------- ## best 100 0 55.2 63.2 52.2
- Model 2，\(n=500\)
as.data.frame(
round(
solver_table(
optMethodSims_getBestVals(
spec2_n500, reslike = TRUE),
"best", spec2) * 100,
digits = 1))
## Maximization Rate Failure Rate omega 95% CI Capture Rate alpha1 95% CI Capture Rate beta1 95% CI Capture Rate ## ----- ------------------ ------------- -------------------------- --------------------------- -------------------------- ## best 100 0 83 86.3 80.5
- Model 2，\(n=1000\)
as.data.frame(
round(
solver_table(
optMethodSims_getBestVals(
spec2_n1000, reslike = TRUE),
"best", spec2) * 100,
digits = 1))
## Maximization Rate Failure Rate omega 95% CI Capture Rate alpha1 95% CI Capture Rate beta1 95% CI Capture Rate ## ----- ------------------ ------------- -------------------------- --------------------------- -------------------------- ## best 100 0 88.7 91.4 88.1
請記住，咱們經過 CI 捕獲率來評估「最佳」優化器的表現，該捕獲率應該在 95％ 左右。「最佳」優化器顯然具備良好的表現，但並不優於全部優化器。這使人失望。我曾但願「最佳」優化器具備捕獲率 95％ 的理想特性。除了較大的樣本量外，表現遠不及預期。標準差被低估，或者對於小樣本，正態分佈很難描述估計量的實際分佈（這意味着標準差乘以 2 不會致使置信區間具備所需置信水平）。
有趣的是，對於這種「最佳」估計器，兩種模型之間的表現沒有明顯差別。這啓示我，在實際數據中常見模型看似更好的結果可能正在利用優化器的誤差。
讓咱們看一下估計參數的分佈。
- \(\omega\) 估計，model 1
solver_density_plot(
"omega", "best",
lapply(
list(spec1_n100, spec1_n500, spec1_n1000),
function(l)
{
optMethodSims_getBestVals(
l, reslike = TRUE)
}),
c(100, 500, 1000), spec1)
- \(\alpha\) 估計，model 1
solver_density_plot(
"alpha1", "best",
lapply(
list(spec1_n100, spec1_n500, spec1_n1000),
function(l)
{
optMethodSims_getBestVals(
l, reslike = TRUE)
}),
c(100, 500, 1000), spec1)
- \(\beta\) 估計，model 1
solver_density_plot(
"beta1", "best",
lapply(
list(spec1_n100, spec1_n500, spec1_n1000),
function(l)
{
optMethodSims_getBestVals(
l, reslike = TRUE)
}),
c(100, 500, 1000), spec1)
- \(\omega\) 估計，model 2
solver_density_plot(
"omega", "best",
lapply(
list(spec2_n100, spec2_n500, spec2_n1000),
function(l)
{
optMethodSims_getBestVals(
l, reslike = TRUE)
}),
c(100, 500, 1000), spec2)
- \(\alpha\) 估計，model 2
solver_density_plot(
"alpha1", "best",
lapply(
list(spec2_n100, spec2_n500, spec2_n1000),
function(l)
{
optMethodSims_getBestVals(
l, reslike = TRUE)
}),
c(100, 500, 1000), spec2)
- \(\beta\) 估計，model 2
solver_density_plot(
"beta1", "best",
lapply(
list(spec2_n100, spec2_n500, spec2_n1000),
function(l)
{
optMethodSims_getBestVals(
l, reslike = TRUE)
}),
c(100, 500, 1000), spec2)
這些圖代表，「最佳」估計器仍然顯示出一些病態，即便它的表現不如其餘估計器差。不管模型選擇如何，我都沒有看到參數估計有偏的證據，但我不相信「最佳」估計器真正最大化對數似然，特別是對於較小的樣本量。\(\beta\) 的估計值特別糟糕。即便 \(\beta\) 的標準差應該很大，我也不認爲它應該像圖中揭示的那樣向 0 或 1 傾斜。
結論
我最初在一年前寫過這篇文章，直到如今才發表。中斷的緣由是由於我想獲得一個關於估計 GARCH 模型參數替代方法的文獻綜述。不幸的是，我從未完成過這樣的綜述，並且我已經決定發佈這篇文章。
那就是說，我會分享我正在讀的東西。Gilles Zumbach 的一篇文章試圖解釋爲何估計 GARCH 參數很難。他指出，求解器試圖最大化的準似然方程具備不良特性，例如非凸，且具備可能讓算法陷入困境的「平坦」區域。他提出了另外一種尋找 GARCH 模型參數的方法，在一個替代參數空間中找到最佳擬合（假設它具備比所使用 GARCH 模型的原始參數空間更好的屬性），而且使用例如矩方法估計其中一個參數，而沒有任何優化算法。另外一篇由 Fiorentini、Calzolari 和 Panattoni 撰寫的文章代表，GARCH 模型的解析梯度能夠明確計算，所以這裏看到的優化算法所使用的無梯度方法實際上並非必需的。因爲數值微分一般是一個難題，這能夠幫助確保不會引入致使這些算法沒法收斂的額外數值偏差。我還想探索其餘估計方法，看看它們是否可以以某種方式徹底避免數值技術，或具備更好的數值屬性，例如經過矩估計。我想閱讀 Andersen、Chung 和 Sørensen 撰寫的文章，以瞭解有關這種估計方法的更多信息。
然而，生活正在繼續，我沒有完成這篇綜述。項目繼續進行，基本上很好地避免了估計 GARCH 模型參數的問題。也就是說，我想從新審視這一點，或許能夠探索模擬退火等技術如何用於估計 GARCH 模型參數。
因此如今，若是你是一名從業者，在估計 GARCH 模型時你應該怎麼作？我想說不要理所固然地認爲你的包使用的默認估計算法會起做用。你應該探索不一樣的算法和不一樣的參數選擇，並使用致使最大對數似然的結果。我展現瞭如何以自動化方式完成這項工做，但你應該準備手動選擇最佳的模型（由對數似然肯定）。若是你不這樣作，你估計的模型實際上可能不是理論可行的模型。
我將在本文的最後再次說一遍，特別強調：不要將數值技術和結果視爲理所固然的！
sessionInfo()
## R version 3.4.2 (2017-09-28) ## Platform: i686-pc-linux-gnu (32-bit) ## Running under: Ubuntu 16.04.2 LTS ## ## Matrix products: default ## BLAS: /usr/lib/libblas/libblas.so.3.6.0 ## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0 ## ## 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] parallel stats graphics grDevices utils datasets methods ## [8] base ## ## other attached packages: ## [1] ggplot2_2.2.1 rugarch_1.3-8 printr_0.1 ## ## loaded via a namespace (and not attached): ## [1] digest_0.6.16 htmltools_0.3.6 ## [3] SkewHyperbolic_0.3-2 expm_0.999-2 ## [5] scales_0.5.0 DistributionUtils_0.5-1 ## [7] Rsolnp_1.16 rprojroot_1.2 ## [9] grid_3.4.2 stringr_1.3.1 ## [11] knitr_1.17 numDeriv_2016.8-1 ## [13] GeneralizedHyperbolic_0.8-1 munsell_0.4.3 ## [15] pillar_1.3.0 tibble_1.4.2 ## [17] compiler_3.4.2 highr_0.6 ## [19] lattice_0.20-35 labeling_0.3 ## [21] Matrix_1.2-8 KernSmooth_2.23-15 ## [23] plyr_1.8.4 xts_0.10-0 ## [25] spd_2.0-1 zoo_1.8-0 ## [27] stringi_1.2.4 magrittr_1.5 ## [29] reshape2_1.4.2 rlang_0.2.2 ## [31] rmarkdown_1.7 evaluate_0.10.1 ## [33] gtable_0.2.0 colorspace_1.3-2 ## [35] yaml_2.1.14 tools_3.4.2 ## [37] mclust_5.4.1 mvtnorm_1.0-6 ## [39] truncnorm_1.0-7 ks_1.11.3 ## [41] nloptr_1.0.4 lazyeval_0.2.1 ## [43] crayon_1.3.4 backports_1.1.1 ## [45] Rcpp_1.0.0
相關文章
- 1. R語言基於ARMA-GARCH過程的VaR擬合和預測
- 2. MLP,MAP,貝葉斯估計在NLP中參數估計
- 3. minepy 包——基於最大信息的非參數估計
- 4. R語言GARCH-DCC模型和DCC（MVT）建模估計
- 5. 關於程序包不存在問題
- 6. 參數的點估計問題與矩估計法
- 7. 關於參數估計
- 8. 機器學習中的參數估計問題
- 9. 《數學建模：基於R》——1.2 參數的區間估計與假設檢驗
- 10. 參數估計
- 更多相關文章...
- • 在Spring中使用Redis - Redis教程
- • Redis在Java Web中的應用 - Redis教程
- • TiDB 在摩拜單車在線數據業務的應用和實踐
- • ☆基於Java Instrument的Agent實現
相關標籤/搜索
每日一句
-
每一个你不满意的现在，都有一个你没有努力的曾经。
歡迎關注本站公眾號,獲取更多信息
相關文章