目录
老树新芽——矩估计遇到神经网络
摘要：本文将矩估计和神经网络结合，为估计 \(\text{GARCH}(1,1)\) 模型的参数提出了一个新的计算框架，并提供了代码实现和计算试验结果。
关键词：矩估计、神经网络、\(\text{GARCH}(1,1)\) 模型
问题
先前的两篇博客文章（文 1、文 2）指出了一般大小样本量的前提下 \(\text{GARCH}(1,1)\) 模型参数估计值的不稳定性，同时指出样本量极端大时估计值才趋于稳定。
但是那两篇博文仅限于指出问题，没有提供解决方案。
解决方案
通常来讲，统计模型的参数估计存在三种主要“套路”：
- 极大似然估计及其变种，通常问题最终被转化为一个数值最优化问题，前面博文中使用的 R 软件包即属于这种套路；
- 矩估计及其变种，通常问题最终被转化为一个解方程（组）的问题，方程（组）可能是非线性的，最终的估计结果可能是一个有偏估计；
- MCMC 及其变种，通常问题最终被转化为一个精心设计过的随机模拟问题。
本文的方法论遵循“矩估计”的套路，文章接下来的内容展示如何将矩估计中解（非线性）方程组这关键步转化成一个神经网络训练问题，并由此提出一个比较通用的“神经矩估计”计算框架。
\(\text{GARCH}(1,1)\) 模型的神经矩估计设想
一个 \(\text{GARCH}(1,1)\) 模型表示如下
\[ \epsilon_t = \varepsilon_t \sigma_t \\ \sigma_t^2 = \omega + \alpha \epsilon_{t-1}^2 + \beta \sigma_{t-1}^2 \\ \varepsilon_t \sim N(0,1); \alpha, \beta \ge 0 \]
\(\text{GARCH}(1,1)\) 模型的神经矩估计包括以下几个步骤：
- 随机生成一组合理的参数 \((\omega, \alpha, \beta)\)；
- 模拟一组长度为 \(N\) 的样本；
- 计算样本的一组矩 \((M_{o_1}, \dots, M_{o_m})\)，其中 \(o_m\) 是矩的阶数；
- 重复第 1～3 步，得到若干矩和模拟参数的配对；
- 以矩为输入数据、模拟参数为输出数据训练神经网络；
- 对合理的参数的空间做网格切分；
- 从切分中选择一组合理的参数 \((\omega_{fix}, \alpha_{fix}, \beta_{fix})\)；
- 模拟 \(n\) 组长度为 \(N\) 的样本；
- 计算上述样本的矩 \((M_{o_1}^i, \dots, M_{o_m}^i), i = 1, \dots, n\)；
- 用神经网络预测上述矩对应的参数 \((\omega_{forecast}^i, \alpha_{forecast}^i, \beta_{forecast}^i), i = 1, \dots, n\)；
- 计算 \((\omega_{forecast}^i, \alpha_{forecast}^i, \beta_{forecast}^i)\) 的均值 \((\bar\omega_{forecast}, \bar\alpha_{forecast}, \bar\beta_{forecast})\)
- 遍历网格切分，得到若干 \((\omega_{fix}, \alpha_{fix}, \beta_{fix})\) 和 \((\bar\omega_{forecast}, \bar\alpha_{forecast}, \bar\beta_{forecast})\) 的配对；
- 用多元插值描述上述 \((\bar\omega_{forecast}, \bar\alpha_{forecast}, \bar\beta_{forecast})\) 到 \((\omega_{fix}, \alpha_{fix}, \beta_{fix})\) 的经验关系。
经过上面的 13 步，一个神经矩估计框架就建好了。下面做一下解释：
- 其中第 1～4 步是为训练神经网络生成训练数据，目的是为第 5 步服务，随机生成模型参数是为了尽可能均匀地遍历参数空间。
- 第 5 步将“解方程组”转换为一个神经网络训练问题，这样的原因有两个：（1）矩 \(M_{o}\) 可能是参数 \((\omega, \alpha, \beta)\) 的非线性函数，甚至得不到显式表达式；（2）从数值计算的角度看，神经网络的训练可以看作是在近似一个多维空间中的函数1。通过训练神经网络就可以忽略 \(M_{o}\) 的表达形式，以一种机械的手段近似出 \((M_{o_1}, \dots, M_{o_m})\) 与 \((\omega, \alpha, \beta)\) 的关系，达到解方程的目的，需要注意的是这里得到的是比较粗糙的近似解，不是精确解；
- 由于第 5 步得到的估计关系可能是有偏的，且不是精确解，为了比较精细地描述未知的偏差，需要得到预测结果均值和真实值之间的插值关系，第 6～13 步正是为此。
如果要估计某一组数据的 \(\text{GARCH}(1,1)\) 参数需要遵循下面的步骤：
- 确保数据的长度接近 \(N\)；
- 计算一组矩 \((M_{o_1}, \dots, M_{o_m})\)；
- 用训练好的神经网络预测得到 \((\omega_{forecast}, \alpha_{forecast}, \beta_{forecast})\)；
- 通过多元插值计算最终的 \((\omega, \alpha, \beta)\)，这一步主要为了纠偏。
代码实现
载入相关包。
# Packages ----
library(ggplot2)
library(fGarch)
library(keras)
library(dplyr)配置必要的函数和对象，作为一次简化的实验，只计算三种矩，分别是：
- \(M_1 = \frac{1}{N}\sum_t|\epsilon_t| = f_1(\omega, \alpha, \beta)\)
- \(M_2 = \frac{1}{N}\sum_t|\epsilon_t|^2 = f_2(\omega, \alpha, \beta)\)
- \(M_{1.5} =\frac{1}{N}\sum_t|\epsilon_t|^{1.5} = f_{1.5}(\omega, \alpha, \beta)\)
# Setup ----
set.seed(123)
M1 <- function(x)
{
sum(abs(x)) / length(x)
}
M2 <- function(x)
{
sum(x^2) / length(x)
}
M1_5 <- function(x)
{
sum(sqrt(abs(x)^3)) / length(x)
}
GarchCoefGenerator <- function(max_omega = 0.5,
up_limit = 1)
{
omega <- runif(1, min = 0, max = max_omega)
a