经验首页 前端设计 程序设计 Java相关 移动开发 数据库/运维 软件/图像 大数据/云计算 其他经验
当前位置：技术经验 » 程序设计 » R语言 » 查看文章
老树新芽——矩估计遇到神经网络
来源：cnblogs　　作者：xuruilong100　　时间：2018/12/24 10:34:49　　对本文有异议

目录

老树新芽——矩估计遇到神经网络

摘要：本文将矩估计和神经网络结合，为估计 \(\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)\) 模型的神经矩估计包括以下几个步骤：

  1. 随机生成一组合理的参数 \((\omega, \alpha, \beta)\)；
  2. 模拟一组长度为 \(N\) 的样本；
  3. 计算样本的一组矩 \((M_{o_1}, \dots, M_{o_m})\)，其中 \(o_m\) 是矩的阶数；
  4. 重复第 1～3 步，得到若干矩和模拟参数的配对；
  5. 以矩为输入数据、模拟参数为输出数据训练神经网络；
  6. 对合理的参数的空间做网格切分；
  7. 从切分中选择一组合理的参数 \((\omega_{fix}, \alpha_{fix}, \beta_{fix})\)；
  8. 模拟 \(n\) 组长度为 \(N\) 的样本；
  9. 计算上述样本的矩 \((M_{o_1}^i, \dots, M_{o_m}^i), i = 1, \dots, n\)；
  10. 用神经网络预测上述矩对应的参数 \((\omega_{forecast}^i, \alpha_{forecast}^i, \beta_{forecast}^i), i = 1, \dots, n\)；
  11. 计算 \((\omega_{forecast}^i, \alpha_{forecast}^i, \beta_{forecast}^i)\) 的均值 \((\bar\omega_{forecast}, \bar\alpha_{forecast}, \bar\beta_{forecast})\)
  12. 遍历网格切分，得到若干 \((\omega_{fix}, \alpha_{fix}, \beta_{fix})\) 和 \((\bar\omega_{forecast}, \bar\alpha_{forecast}, \bar\beta_{forecast})\) 的配对；
  13. 用多元插值描述上述 \((\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)\) 参数需要遵循下面的步骤：

  1. 确保数据的长度接近 \(N\)；
  2. 计算一组矩 \((M_{o_1}, \dots, M_{o_m})\)；
  3. 用训练好的神经网络预测得到 \((\omega_{forecast}, \alpha_{forecast}, \beta_{forecast})\)；
  4. 通过多元插值计算最终的 \((\omega, \alpha, \beta)\)，这一步主要为了纠偏。

代码实现

载入相关包。

  1. # Packages ----
  3. library(ggplot2)
  4. library(fGarch)
  5. library(keras)
  6. 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)\)
  1. # Setup ----
  3. set.seed(123)
  5. M1 <- function(x)
  6. {
  7.     sum(abs(x)) / length(x)
  8. }
  10. M2 <- function(x)
  11. {
  12.     sum(x^2) / length(x)
  13. }
  15. M1_5 <- function(x)
  16. {
  17.     sum(sqrt(abs(x)^3)) / length(x)
  18. }
  20. GarchCoefGenerator <- function(max_omega = 0.5,
  21.                                up_limit = 1)
  22. {
  23.     omega <- runif(1, min = 0, max = max_omega)
  24.     alpha <- runif(1)
  25.     beta <- runif(1)
  27.     while (alpha + beta >= up_limit)
  28.     {
  29.         alpha <- runif(1)
  30.         beta <- runif(1)
  31.     }
  33.     return(
  34.         list(omega = omega, alpha = alpha, beta = beta))
  35. }
  37. train_length <- 3000
  38. val_length <- 1000
  39. test_length <- 1000
  41. train_m1 <- numeric(train_length)
  42. train_m2 <- numeric(train_length)
  43. train_m1_5 <- numeric(train_length)
  44. train_omega <- numeric(train_length)
  45. train_alpha <- numeric(train_length)
  46. train_beta <- numeric(train_length)
  48. val_m1 <- numeric(val_length)
  49. val_m2 <- numeric(val_length)
  50. val_m1_5 <- numeric(val_length)
  51. val_omega <- numeric(val_length)
  52. val_alpha <- numeric(val_length)
  53. val_beta <- numeric(val_length)
  55. test_m1 <- numeric(test_length)
  56. test_m2 <- numeric(test_length)
  57. test_m1_5 <- numeric(test_length)

生成训练数据。

  1. # Train Data ----
  3. for (i in 1:train_length)
  4. {
  5.     garch_coef <- GarchCoefGenerator()
  6.     garch_sample <- garchSim(
  7.         garchSpec(
  8.             model = garch_coef),
  9.         n.start = 1000, n = 1000)
  11.     train_m1[i] <- M1(garch_sample)
  12.     train_m2[i] <- M2(garch_sample)
  13.     train_m1_5[i] <- M1_5(garch_sample)
  15.     train_omega[i] <- garch_coef$omega
  16.     train_alpha[i] <- garch_coef$alpha
  17.     train_beta[i] <- garch_coef$beta
  18. }
  20. trainX <- data.frame(
  21.     m1 = train_m1,
  22.     m2 = train_m2,
  23.     m1_5 = train_m1_5) %>%
  24.     as.matrix()
  26. trainY <- data.frame(
  27.     omega = train_omega,
  28.     alpha = train_alpha,
  29.     beta = train_beta) %>%
  30.     as.matrix()

生成验证数据。

  1. # Val Data ----
  3. for (i in 1:val_length)
  4. {
  5.     garch_coef <- GarchCoefGenerator()
  6.     garch_sample <- garchSim(
  7.         garchSpec(
  8.             model = garch_coef),
  9.         n.start = 1000, n = 1000)
  11.     val_m1[i] <- M1(garch_sample)
  12.     val_m2[i] <- M2(garch_sample)
  13.     val_m1_5[i] <- M1_5(garch_sample)
  15.     val_omega[i] <- garch_coef$omega
  16.     val_alpha[i] <- garch_coef$alpha
  17.     val_beta[i] <- garch_coef$beta
  18. }
  20. valX <- data.frame(
  21.     m1 = val_m1,
  22.     m2 = val_m2,
  23.     m1_5 = val_m1_5) %>%
  24.     as.matrix()
  26. valY <- data.frame(
  27.     omega = val_omega,
  28.     alpha = val_alpha,
  29.     beta = val_beta) %>%
  30.     as.matrix()

用 keras 构建神经网络。

  1. # Model Training ----
  3. model <- keras_model_sequential()
  5. model %>%
  6.     layer_dense(
  7.         units = 9,
  8.         input_shape = c(3),
  9.         activation = 'sigmoid') %>%
  10.     layer_dense(
  11.         units = 3,
  12.         activation = 'sigmoid') %>%
  13.     compile(
  14.         loss = 'mae',
  15.         optimizer = 'adam')
  17. summary(model)

实验采用了结构非常简单的神经网络结构，一个两层的多元感知机，超参数、激活函数和优化方法没有调优。

  1. __________________________________________________________________________
  2. Layer (type)                     Output Shape                 Param #     
  3. ==========================================================================
  4. dense_1 (Dense)                  (None, 9)                    36          
  5. __________________________________________________________________________
  6. dense_2 (Dense)                  (None, 3)                    30          
  7. ==========================================================================
  8. Total params: 66
  9. Trainable params: 66
  10. Non-trainable params: 0
  11. __________________________________________________________________________

训练。

  1. history <- model %>%
  2.     fit(
  3.         trainX,
  4.         trainY,
  5.         validation_data = list(
  6.             valX, valY),
  7.         batch_size = 50,
  8.         epochs = 1000)
  10. plot(history)

训练结果非常理想，训练集和验证集上的 loss 基本上同步衰减到 0.11 左右，并且从趋势上看大约已经接近了极限水平。

延续前面两篇博文的主题，将 \(\text{GARCH}(1,1)\) 的模型参数限定为 \((0.2,0.2,0.2)\)，生成测试数据，后面要测试估计值的稳定性，以及是否有偏。

  1. # Test ----
  3. for (i in 1:test_length)
  4. {
  5.     garch_coef <- list(
  6.         omega = 0.2, alpha = 0.2, beta = 0.2)
  8.     garch_sample <- garchSim(
  9.         garchSpec(
  10.             model = garch_coef),
  11.         n.start = 1000, n = 1000)
  13.     test_m1[i] <- M1(garch_sample)
  14.     test_m2[i] <- M2(garch_sample)
  15.     test_m1_5[i] <- M1_5(garch_sample)
  16. }
  18. predictX <- data.frame(
  19.     m1 = test_m1,
  20.     m2 = test_m2,
  21.     m1_5 = test_m1_5) %>%
  22.     as.matrix()

预测结果。

  1. predictY <- model %>%
  2.     predict(predictX) %>%
  3.     as.data.frame()
  5. predictY$idx <- 1:nrow(predictY)
  7. colnames(predictY) <- c('omega', 'alpha', 'beta', 'idx')
  9. head(predictY)
  1.       omega     alpha      beta idx
  2. 1 0.1946828 0.1991726 0.2187224   1
  3. 2 0.1677427 0.2192231 0.2021861   2
  4. 3 0.1989923 0.2052666 0.2165215   3
  5. 4 0.1739261 0.2159361 0.2049208   4
  6. 5 0.1766528 0.2010799 0.2137308   5
  7. 6 0.2087296 0.1955684 0.2248684   6

目前预测结果看起来还不错。

验证估计结果

用红线代表估计值的均值，蓝线代表理论值。

\(\omega\)

  1. fig_scatter <- ggplot(
  2.     data = predictY) +
  3.     geom_point(
  4.         mapping = aes(
  5.             x = idx, y = omega)) +
  6.     geom_hline(yintercept = 0.2, color = 'blue') +
  7.     geom_hline(yintercept = mean(predictY$omega), color = 'red') +
  8.     ylim(c(0.1, 0.3))
  10. fig_scatter
  12. fig_density <- ggplot(
  13.     data = predictY) +
  14.     geom_density(
  15.         mapping = aes(
  16.             x = omega)) +
  17.     geom_vline(xintercept = 0.2, color = 'blue') +
  18.     geom_vline(xintercept = mean(predictY$omega), color = 'red') +
  19.     xlim(c(0.1, 0.3))
  21. fig_density
  23. plot_grid(
  24.     fig_scatter, fig_density,
  25.     nrow = 2, labels = c('Scatter', 'Density'))

\(\omega\) 的估计值非常稳定，变动区间也很小。估计是有偏的，估计值会系统性地低估真实值，差距在 0.015 左右。

\(\alpha\)

  1. fig_scatter <- ggplot(
  2.     data = predictY) +
  3.     geom_point(
  4.         mapping = aes(
  5.             x = idx, y = alpha)) +
  6.     geom_hline(yintercept = 0.2, color = 'blue') +
  7.     geom_hline(yintercept = mean(predictY$alpha), color = 'red') +
  8.     ylim(c(0.1, 0.3))
  10. fig_scatter
  12. fig_density <- ggplot(
  13.     data = predictY) +
  14.     geom_density(
  15.         mapping = aes(
  16.             x = alpha)) +
  17.     geom_vline(xintercept = 0.2, color = 'blue') +
  18.     geom_vline(xintercept = mean(predictY$alpha), color = 'red') +
  19.     xlim(c(0.1, 0.3))
  21. fig_density
  23. plot_grid(
  24.     fig_scatter, fig_density,
  25.     nrow = 2, labels = c('Scatter', 'Density'))

\(\alpha\) 的估计值非常稳定，变动区间也很小。估计是有偏的，估计值会系统性地高估真实值，差距在 0.02 左右。

\(\beta\)

  1. fig_scatter <- ggplot(
  2.     data = predictY) +
  3.     geom_point(
  4.         mapping = aes(
  5.             x = idx, y = beta)) +
  6.     geom_hline(yintercept = 0.2, color = 'blue') +
  7.     geom_hline(yintercept = mean(predictY$beta), color = 'red') +
  8.     ylim(c(0.1, 0.3))
  10. fig_scatter
  12. fig_density <- ggplot(
  13.     data = predictY) +
  14.     geom_density(
  15.         mapping = aes(
  16.             x = beta)) +
  17.     geom_vline(xintercept = 0.2, color = 'blue') +
  18.     geom_vline(xintercept = mean(predictY$beta), color = 'red') +
  19.     xlim(c(0.1, 0.3))
  21. fig_density
  23. plot_grid(
  24.     fig_scatter, fig_density,
  25.     nrow = 2, labels = c('Scatter', 'Density'))

\(\beta\) 的估计值非常稳定，变动区间也很小。估计是有偏的，估计值会系统性地高估真实值，差距在 0.006 左右，已经很小了。

一个悬而未决的理论问题

从工程的角度看，前面的实践是没问题的，但有一个理论问题悬而未决，那就是样本矩的存在性。一个潜在的论证思路如下：

  1. 研究 \(\text{GARCH}(1,1)\) 过程 \(\epsilon_t\) 矩的存在条件；
  2. 论证过程 \(|\epsilon_t|^o(1 \le o \le 2)\) 的遍历性；
  3. 结合以上两点，用 \(\epsilon_t\) 矩的存在性加上 \(|\epsilon_t|^o\) 的遍历性证明样本矩的存在性。

\(\text{GARCH}(1,1)\) 模型矩的统计性质

\(\text{GARCH}(1,1)\) 过程 \(\epsilon_t\) 存在 \(2m\) 阶矩的充分必要条件是2

\[ \mu(\alpha, \beta, m) = \sum_{j=0}^m C_m^j a_j \alpha^j \beta^{m-j} < 1 \]

其中，\(m\) 是正整数；\(a_0=1\)，\(a_j = \prod_{i=1}^j (2i-1), j = 1,\dots,m\)。

如果要求 \(\epsilon_t\) 有 4 阶矩，即 \(m = 2\)，那么

\[ (\alpha + \beta)^2 + 2\alpha^2 < 1 \]

模型参数的取值范围比较狭小，仅局限在原点为中心的一个椭圆内。

如果放松要求，让 \(\epsilon_t\) 仅有 2 阶矩，即 \(m = 1\)，那么

\[ \alpha + \beta < 1 \]

这已经是普通 \(\text{GARCH}(1,1)\) 模型能达到的最大范围了，如果 \(\alpha + \beta = 1\)，模型将变为没有二阶矩的 \(\text{I-GARCH}(1,1)\) 模型。

\(\text{GARCH}(1,1)\) 模型的遍历性

目前还没有找到合适的文献，可能需要在 \(\alpha + \beta < 1\) 的基础上在近一步缩小取值范围。

如果读者中有人对此问题有较深入的见解，请不吝赐教，我的邮箱是：xuruilong100@163.com

结论

受限于个人家用电脑的计算能力，我的数值实验只能做到这一步，毕竟框架的后半部分是一个计算密集型任务。

我对目前的试验结果很满意，神经矩估计在一般大小样本上表现出的稳定性远超基于极大似然方法论的经典方法，对于工程实践来说，现在到达的计算准确度基本够用了。

对于样本矩的存在性还需要进一步的研究。

如果有足够的计算能力验证框架的后半部分是有效且实用的，神经矩估计方法有潜力推广成为一个估计复杂统计（随机）模型参数的普适框架。

参考文献

  1. 《Neural Networks and Deep Learning》（Michael Nielsen），第四章。?

  2. 《协整理论与波动模型（第二版）》（张士英、樊智），定理 6-4。?

 友情链接：直通硅谷 　点职佳 　北美留学生论坛

R语言热门文章

  • R语言初学者的一些常见报错指南
  • Installation failed: Timeout was reached: Operation timed out after 10000 milliseconds with 0 out of 0 bytes received
  • R语言数据重塑知识点总结
  • R语言dataframe的常用操作总结
  • 在 R 中估计 GARCH 参数存在的问题（基于 rugarch 包）
  • 《R语言实战（第2版）》高清中文版PDF+高清英文版PDF+源代码
  • PCA
  • R语言科学计数法介绍:digits和scipen设置方式
  • R语言绘制带误差线的条形图
  • 通过R语言统计考研英语（二）单词出现频率

    • R语言推荐文章

  • 随机森林R语言预测工具
  • 地理探测器R语言实现：geodetector
  • R语言初学者的一些常见报错指南
  • R语言学习之线图的绘制详解
  • R语言学习笔记缺失数据的Bootstrap与Jackknife方法
  • R语言关于数据帧的知识点详解
  • 机器学习 - 相关概念与实现流程
  • PCA
  • R语言构建蛋白质网络并实现GN算法
  • 面试题7 指定一个m*n的整数矩阵 如果(r, c)位置元素的元素是0 则将行 r 和列 c 整个变成0 - 雄霸天下-无人能挡

    • 本站QQ群：前端 618073944 | Java 606181507 | Python 626812652 | C/C++ 612253063 | 微信 634508462 | 苹果 692586424 | C#/.net 182808419 | PHP 305140648 | 运维 608723728

    W3xue 的所有内容仅供测试，对任何法律问题及风险不承担任何责任。通过使用本站内容随之而来的风险与本站无关。
    关于我们  |  意见建议  |  捐助我们  |  报错有奖  |  广告合作、友情链接（目前9元/月）请联系QQ：27243702 沸活量
    皖ICP备17017327号-2　皖公网安备34020702000426号