在 R 中估计 GARCH 参数存在的问题（基于 rugarch 包）

本文翻译自《Problems in Estimating GARCH Parameters in R (Part 2; rugarch)》

原文链接：https://ntguardian.wordpress.com/2019/01/28/problems-estimating-garch-parameters-r-part-2-rugarch/

导论

这是一篇本应早就写完的博客文章。一年前我写了一篇文章，关于在 R 中估计 GARCH(1, 1) 模型参数时遇到的问题。我记录了参数估计的行为（重点是 \(\beta\)），以及使用 fGarch 计算这些估计值时发现的病态行为。我在 R 社区呼吁帮助，包括通过 R Finance 邮件列表发送我的博客文章。

反馈没有让我感到失望。你可以看到一些邮件列表反馈，并且一些来自 Reddit 的评论也很有帮助，但我认为我得到的最佳反馈来自于我自己的电子邮件。

Dr. Brian G. Peterson 作为 R 金融社区的一员，给我发送了一些发人深思的电子邮件。首先，他告诉我 fGarch 不再是处理 GARCH 模型的首选方案。RMetrics 套件包（包括 fGarch）由 ETH Zürich 的 Diethelm Würtz 教授维护。他在 2016 年的车祸中丧生。

Dr. Peterson 建议我研究另外两个用于 GARCH 建模的现代软件包，rugarch（适用于单变量 GARCH 模型）和 rmgarch（适用于多变量 GARCH 模型）。之前我没有听说过这些包（我之所以知道 fGarch 的原因是因为它在由 Shumway 和 Stoffer 编写的时间序列教科书——Time Series Analysis and Its Applications with R Examples 中），所以我非常感谢这个建议。由于我现在对单变量时间序列感兴趣，所以我研究了 rugarch。该软件包似乎具有比 fGarch 更多的功能和函数，这可以解释为什么它似乎更难以使用。然而，包的 vignette 很有帮助，值得打印出来。

Dr. Peterson 对我提出的应用也有一些有趣的评论。他认为，日内数据应优于日间数据，并且模拟数据（包括模拟 GARCH 过程）具有在实际数据中看不到的特质。获取日间数据的便利性（特别是亚洲金融危机期间的 USD/JPY，这是我正在研究的检验统计量的预期应用）激发了我对日间数据的兴趣。不过，他的评论可能会让我重新考虑这个应用。¹（我也许应该试图通过 EUR/USD 来检测 2010 年欧元区金融危机。为此，我可以从 HistData.com 获得免费的日内数据。）但是，如果对于小样本而言不能信任标准差的估计，我们的检验统计量仍然会遇到麻烦，因为它涉及小样本的参数估计。

他还警告说，模拟数据表现出在实际数据中看不到的行为。这可能是真的，但模拟数据很重要，因为它可以被认为是统计学家的最佳情景。另外，生成模拟数据的过程的属性是先验已知的，包括生成参数的值，以及哪些假设（例如序列中是否存在结构变化）是真的。这允许对估计器和检验进行健全的检查。这对现实世界来说是不可能的，因为我们没有所需的先验知识。

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）。(可能有也可能没有，我没有检查过。如果有人知道，请分享。)

我的导师联系了另一位 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 Rao² 在 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. 当我最初写这篇文章时，我的导师和他的前学生开发了一个检验统计量，应该检测时间序列中的早期或晚期变点，包括 GARCH 模型参数的变化。我对我们所撰写论文的贡献包括证明，当应用于真实世界数据时，检验统计量比其他检验统计量更快地检测到结构变化。为了说服审阅者，我们的检验统计量应检测到另一个统计量在获得更多数据之前无法检测到的变化。这意味着变化应该存在，但不会太强，以至于两个统计数据都可以立即通过微小的 \(p\)-value 检测到变化。?

2. 我链接到的 LinkedIn 上的个人资料可能是也可能不是正确的人，我猜这是基于列出的职业和历史。如果我找错了人，我很抱歉。?