38、时间序列分析与预测：从模型构建到蒙特卡罗模拟

时间序列分析与预测：从模型构建到蒙特卡罗模拟

1. 时间序列模型初步处理

在对时间序列数据进行分析时，我们会发现其部分自相关图在滞后 12 处有一个峰值，且不同于自相关图，它在滞后 24 处没有峰值，因为周期性自相关在滞后 12 处已被考虑。虽然这似乎表明 AR(12) 模型可能合适，但它会产生大量需要学习的系数，尤其是在数据量相对较少的情况下。由于周期性是季节性的，我们可以通过二次差分来去除它。

我们已经对数据进行了一次差分，此时模型被称为自回归积分滑动平均（ARIMA）模型，差分阶数用参数 d 表示，完整的模型阶数可指定为 ARIMA(p, d, q)。以下是第二次差分以去除数据中强烈季节性的代码：

(defn ex-9-27 []
  (->> (airline-passengers)
       (difference)
       (difference 12)
       (autocorrelation)
       (take 15)
       (bar-plot)))

接下来是绘制部分自相关图的代码：

(defn ex-9-28 []
  (->> (airline-passengers)
       (difference)
       (difference 12)
       (partial-autocorrelation)
       (take 15)
       (bar-plot)))

经过处理后，强季节性周期占据了图表中大部分的显著性。两个图表在滞后 1 处都存在负自相关，自相关函数（ACF）在滞后 9 处有一个勉强显著的自相关。一般来说，正自相关最好通过在模型中添加 AR 项来处理，而负自相关通常最好通过添加 MA 项来处理。基于上述图表，一个合理的模型是 MA(1) 模型，但我们可以借此机会尝试捕捉 AR(9) 自相关，以展示如何为模型拟合大量参数。

2. 最大似然估计

在之前的一些优化问题中，我们通常用成本函数最小化来表达，例如在构建逻辑回归分类器时使用 Incanter 最小化逻辑成本函数，在进行批量和随机梯度下降时使用梯度下降最小化最小二乘成本函数。但优化也可以表示为最大化收益，最大似然估计就是通过最大化似然函数来找到模型的最佳参数。

假设给定模型参数 β 时观察值 x 的概率表示为 (P(x|\beta))，那么似然可以表示为 (L(x|\beta))，它衡量了在给定数据的情况下参数的概率。最大似然估计的目标是找到使观察到的数据最有可能出现的参数值。

在计算时间序列的似然之前，我们通过一个简单的抛硬币例子来说明这个过程。假设抛硬币 100 次，观察到 56 次正面（h）和 44 次反面（t），我们可以问什么值的 (P(h)) 能使观察到的数据最有可能出现。以下是相关代码：

(defn ex-9-29 []
  (let [data 56
        f (fn [p]
            (s/pdf-binomial data :size 100 :prob p))]
    (-> (c/function-plot f 0.3 0.8
                         :x-label "P"
                         :y-label "Likelihood")
        (i/view))))

这里我们使用二项分布来模拟抛硬币序列，关键是数据是固定的，我们绘制在不同二项分布参数下观察到该数据的不同概率。最有可能的二项分布参数是 (p = 0.56)。

对于时间序列，计算其参数的似然数学过程较为复杂，我们将使用 Clojure 库 Succession 来计算。通常我们使用对数似然而不是似然，这是为了数学计算方便，因为对数似然可以将乘积形式转换为求和形式：
(LL(\beta|x)=\sum_{i = 1}^{k}\log f(x_i|\beta))
以下是绘制简单 AR(2) 时间序列不同参数的对数似然图的代码：

(defn ex-9-30 []
  (let [init  (s/sample-normal 2)
        coefs [0 0.5]
        data  (take 100 (ar init coefs 0.2))
        f     (fn [coef]
                (log-likelihood [0 coef] 2 0 data))]
    (-> (c/function-plot f -1 1
                         :x-label "Coefficient"
                         :y-label "Log-Likelihood")
        (i/view))))

曲线的峰值对应于给定数据下参数的最佳估计。

3. 最大似然估计方法：Nelder - Mead 法

由于 ARMA 模型的参数数量较多，为了确定最大似然，我们使用在高维空间中表现良好的 Nelder - Mead 方法（也称为单纯形法）。在 n 维空间中，单纯形是一个具有 n + 1 个顶点的多面体。

单纯形优化的优点是它不需要在每个点计算梯度来下降（或上升）到更优位置。Nelder - Mead 方法通过单纯形上每个测试点的目标函数行为进行外推，用通过其余点质心反射得到的点替换最差的点。如果新点比当前最佳点更好，则沿此线指数拉伸单纯形；如果新点改善不大，可能是跨越了山谷，此时收缩单纯形朝向可能更好的点。

单纯形方法在 Incanter 中未实现，但可以在 Apache Commons Math 库中使用。以下是使用该库进行 Nelder - Mead 优化的代码：

(defn objective-function [f]
  (ObjectiveFunction. (reify MultivariateFunction
                        (value [_ v]
                          (f (vec v))))))

(defn arma-model [p q ys]
  (let [m (+ p q)
        f (fn [params]
            (sc/log-likelihood params p q ys))
        optimal (.optimize (SimplexOptimizer. 1e-10 1e-10)
                           (into-array
                            OptimizationData
                            [(MaxEval. 100000)
                             (objective-function f)
                             GoalType/MAXIMIZE
                             (->> (repeat 0.0)
                                  (take m)
                                  (double-array)
                                  (InitialGuess.))
                             (->> (repeat 0.1)
                                  (take m)
                                  (double-array)
                                  (NelderMeadSimplex.))]))
        point (-> optimal .getPoint vec)
        value (-> optimal .getValue)]
    {:ar (take p point)
     :ma (drop p point)
     :ll value}))

(defn ex-9-31 []
  (->> (airline-passengers)
       (i/log)
       (difference)
       (difference 12)
       (arma-model 9 1)))

运行上述代码后，我们可以得到模型的最大似然估计参数。

4. 赤池信息准则（AIC）选择模型

在评估多个模型时，最大似然估计值最大的模型似乎是最好的，但它没有考虑模型的复杂性，一般来说，更简单的模型更可取。赤池信息准则（AIC）是一种比较模型的方法，它奖励由似然函数评估的拟合优度，但包含一个与参数数量有关的惩罚项，以防止过拟合。

AIC 的计算公式为：(AIC = 2k - 2\log L)，其中 (k) 是模型的参数数量，(L) 是似然函数。以下是在 Clojure 中计算 AIC 的代码：

(defn aic [coefs p q ys]
  (- (* 2 (+ p q 1))
     (* 2 (log-likelihood coefs p q ys))))

如果要生成多个模型并选择最佳模型，我们应该选择 AIC 最低的模型。

以下是一个简单的流程图展示整个过程：

graph TD;
    A[时间序列数据] --> B[差分处理];
    B --> C[绘制自相关和部分自相关图];
    C --> D[确定模型类型];
    D --> E[最大似然估计];
    E --> F[计算 AIC 选择模型];

通过以上步骤，我们完成了从时间序列数据的初步处理到模型选择的过程，为后续的时间序列预测奠定了基础。在后续内容中，我们将介绍如何使用这些模型进行时间序列预测，并通过蒙特卡罗模拟来估计预测的置信区间。

时间序列分析与预测：从模型构建到蒙特卡罗模拟

5. 时间序列预测

在确定了模型的参数估计后，我们终于可以使用模型进行时间序列预测了。我们已经有了一个 arma 函数，它能够根据一些种子数据以及模型参数 p 和 q 生成自回归滑动平均序列。种子数据将是我们从航空数据中测量得到的 y 值，而 p 和 q 的值则是我们使用 Nelder - Mead 方法计算得到的参数。

以下是将这些参数代入 ARMA 模型并生成 y 的预测序列的代码：

(defn ex-9-32 []
  (let [data (i/log (airline-passengers))
        diff-1  (difference 1 data)
        diff-12 (difference 12 diff-1)
        forecast (->> (arma (take 9 (reverse diff-12))
                       []
                       (:ar params)
                       (:ma params) 0)
                      (take 100)
                      (undifference 12 diff-1)
                      (undifference 1 data))]
    (->> (concat data forecast)
         (i/exp)
         (timeseries-plot))))

上述代码生成的图表中，到时间切片 144 的线是原始序列，之后的线是我们的预测序列。预测结果看起来符合预期，指数增长趋势持续，并且有规律的季节性高峰和低谷模式也继续存在。

然而，预测结果过于规则，与 1 到 144 点的原始序列不同，预测中没有噪声。为了使预测更现实，我们需要添加一些噪声。为了确定合理的噪声量，我们可以查看过去预测的误差。为避免误差累积，我们应该一次预测一个时间步，并观察预测值与实际值之间的差异。

以下是添加噪声（ sigma 为 0.02）的代码：

(defn ex-9-33 []
  (let [data (i/log (airline-passengers))
        diff-1  (difference 1 data)
        diff-12 (difference 12 diff-1)
        forecast (->> (arma (take 9 (reverse diff-12))
                       []
                       (:ar params)
                       (:ma params) 0.02)
                      (take 10)
                      (undifference 12 diff-1)
                      (undifference 1 data))]
    (->> (concat data forecast)
         (i/exp)
         (timeseries-plot))))

运行上述代码后，我们可以得到一个带有一定波动性的预测图表。通过多次运行模拟，我们可以了解不同可能结果的多样性。接下来，我们希望能够确定预测的置信区间，即所有未来序列的上下期望，包括噪声。

6. 蒙特卡罗模拟进行预测

虽然存在解析方法来计算时间序列的预期未来值和置信区间，但我们将通过模拟来得到这些值。通过研究大量预测的变化，我们可以得到模型预测的置信区间。

蒙特卡罗方法是一种常用的统计工具，用于解决解析上难以处理的问题。它是在曼哈顿计划期间由 John Von Neumann 和 Stanislaw Ulam 提出的，用于研究中子通过辐射屏蔽的传播特性，并以摩纳哥的蒙特卡罗赌场命名。

我们已经为时间序列预测的蒙特卡罗模拟奠定了基础，只需要运行数百次模拟并收集结果。以下是运行 1000 次模拟并计算 95% 置信区间的代码：

(defn ex-9-34 []
  (let [data (difference (i/log (airline-passengers)))
        init (take 12 (reverse data))
        forecasts (for [n (range 1000)]
                    (take 20
                          (arma init [0.0028 0.0028]
                                (:ar params1)
                                (:ma params1)
                                0.0449)))
        forecast-mean (map s/mean (i/trans forecasts))
        forecast-sd (-> (map s/sd (i/trans forecasts))
                        (i/div 2)
                        (i/mult 1.96))
        upper (->> (map + forecast-mean forecast-sd)
                   (concat data)
                   (undifference 0)
                   (i/exp))
        lower (->> (map - forecast-mean forecast-sd)
                   (concat data)
                   (undifference 0)
                   (i/exp))
        n (count upper)]
    (-> (c/xy-plot   (range n) upper
                     :x-label "Time"
                     :y-label "Value"
                     :series-label "Upper Bound"
                     :legend true)
        (c/add-lines (range n) lower
                     :series-label "Lower Bound")
        (i/view))))

上述代码生成的图表中，上下边界为我们的时间序列未来预测提供了置信区间。

总结

时间序列分析涉及对离散时间序列的分析，即按固定时间间隔进行的顺序观测。通过将时间序列分解为趋势、季节性和周期性等组件，可以更轻松地对其进行建模。ARMA 模型进一步将序列分解为自回归和滑动平均组件，每个组件都在某种程度上由序列的过去值决定。

Clojure 自然的递归函数定义和惰性序列能力非常适合生成此类序列。通过将序列的每个值确定为前一个值的函数，我们实现了一个递归 ARMA 生成器，能够模拟测量序列并向前预测。

我们还学习了最大似然估计，它将优化问题的解决方案重新定义为在给定数据的情况下生成最大似然的方案。通过 Apache Commons Math 库的 Nelder - Mead 方法，我们可以估计最大似然参数。

最后，我们展示了如何通过向前播放序列进行预测，并使用蒙特卡罗模拟来估计序列的未来误差。整个过程可以总结为以下表格：
|步骤|描述|
|----|----|
|数据预处理|对时间序列数据进行差分处理，去除季节性|
|模型选择|通过自相关和部分自相关图确定模型类型，使用最大似然估计和 AIC 选择最佳模型|
|预测|使用 ARMA 模型进行预测，添加噪声使预测更现实|
|置信区间估计|使用蒙特卡罗模拟计算预测的置信区间|

以下是整个时间序列分析与预测过程的流程图：

graph TD;
    A[时间序列数据] --> B[差分处理];
    B --> C[绘制自相关和部分自相关图];
    C --> D[确定模型类型];
    D --> E[最大似然估计];
    E --> F[计算 AIC 选择模型];
    F --> G[时间序列预测];
    G --> H[添加噪声];
    H --> I[蒙特卡罗模拟计算置信区间];

通过这些步骤，我们可以对时间序列数据进行全面的分析和预测，为实际应用提供有力的支持。