在 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实现
相关标签/搜索
每日一句
-
每一个你不满意的现在，都有一个你没有努力的曾经。
欢迎关注本站公众号,获取更多信息
相关文章