时序分析 32

时序预测 从ARIMA到SARIMAX

    上一篇文章中我们介绍了SARIMA模型，相比于ARIMA模型，SARIMA考虑了季节性因素。读过上篇文章的细心读者会发现，python的statsmodels对sarima的支持对象名为SARIMAX。事实上，确实存在一个时序分析模型称为SARIMAX(Seasonal Auto-Regressive Integrated Moving Average with eXogenous factors)，对比于SARIMA，所多的这个’X’代表exogenous，也就是说SARIMAX可考虑其他外部变量来协助对该时序数据进行建模。这种情况也是经常会出现的。

    本篇博文会结合一个例子来重新梳理一下时序分析预测中，从ARIMA -> SARIMA -> SARIMAX的整个过程。顺便说一下，另外一个时序模型名为ARIMAX；不言而喻，该模型就是在ARIMA模型的基础上增加了可考虑外部变量。

    本篇博文中所采用的数据是沃尔玛的销售数据集，可以从kaggle(https://www.kaggle.com/competitions/demand-forecasting-kernels-only)上获得。 此数据集中有5年的销售数据，分位train.csv和test.csv。本节我们的任务就是对此数据集进行时序建模并预测。

理论回顾

    上篇关于SARIMA的文章中，笔者并没有详细介绍SARIMA背后的数学模型。本文我们对比ARIMA模型简单介绍一下。这里面的数学比较复杂，读者不用太过在意，对此无兴趣的读者可略过这段内容。

加性模型和乘性模型

我们先简单介绍一下加性模型和乘性模型。

• 加性模型：$x_t=Trend+Seasonal+Random$
• 乘性模型：$x_t=Trend\ast Seasonal\ast Random$
  这里的 Random 项在时序分解中经常被认为是不规则项，不一定就是随机因素。
  加性模型一般适用于季节性波动随着时间保持稳定；而乘性模型相对适用于季节波动随着时间会变大。

$ARIMA(1,1,1)$:

   ARIMA是一个加性模型(Additive Model)，ARIMA(1,1,1)可以被下式描述，
$$\Delta y_t=c+{\varphi }_1\Delta y_{t-1}+{\theta }_1{ϵ}_{t-1}+{ϵ}_t$$
上式中 $\Delta$ 是一阶差分操作符， ${ϵ}_t\sim N(0,{\sigma }^2)$
,上式可以改写成$$(1-{\varphi }_{1}L)\Delta y_t=c+(1+{\theta }_{1}L){ϵ}_t$$
$L$是滞后算子(lag operator)，其运算规则为$L{y}_t=y_{t-1}，L^2{y}_t=y_{t-2}$

$SARIMA(p,d,q)\times (P,D,Q{)}_s$

   大多数SARIMA模型是一个乘性模型(Multiplicative Model)，可以泛化为下式:
$${\varphi }_p(L){\tilde{\varphi }}_P(L^s){\Delta }^d{\Delta }_s^D{y}_t=A(t)+{\theta }_q(L){\tilde{\theta }}_Q(L^s){ϵ}_t$$
或者
$${\varphi }_p(L){\tilde{\varphi }}_P(L^s)y_t^{\ast }=A(t)+{\theta }_q(L){\tilde{\theta }}_Q(L^s){ϵ}_t$$
$$y_t^{\ast }={\Delta }^d{\Delta }_s^D{y}_t$$

这里，

• ${\varphi }_p(L)$ 是非季节性自回归滞后多项式 ：捕获非季节性自回归因素。
• ${\tilde{\varphi }}_P(L^s)$是季节性自回归滞后多项式：捕获季节性回归因素。
• ${\Delta }^d{\Delta }_s^D{y}_t$ 是时序数据 $d$ 阶差分；季节性 $D$ 阶差分：提供时序平稳化功能。
• $A(t)$ 为趋势多项式：包括截距项
• ${\theta }_q(L)$ 是非季节性移动平均滞后多项式：
• ${\tilde{\theta }}_Q(L^s)$ 是季节性移动平均滞后多项式

假设我们有一个 $SARIMA(2,1,0)\times (1,1,0{)}_{1}2$ 包含截距项， 我们有

• ${\varphi }_p(L)=(1-{\varphi }_{1}L-{\varphi }_2{L}^2)$
• ${\tilde{\varphi }}_P(L^s)=(1-{\varphi }_1{L}^{1}2)$
• $d=1,D=1,s=12$ 意味着 $y_t^{\ast }$ 是经过一阶差分和季节性（12个月）差分所得到的
• $A(t)=c$ 意味着趋势多项式就是一个截距项。
• ${\theta }_q(L)={\tilde{\theta }}_Q(L^s)=1$ 无移动平均效果

那么，SARIMA泛化公式就会变为
$$(1-{\varphi }_{1}L-{\varphi }_2{L}^2)(1-{\tilde{\varphi }}_1{L}^{12})\Delta {\Delta }_{12}{y}_t=c+{ϵ}_t$$
$=>$
$$(1-{\varphi }_{1}L-{\varphi }_2{L}^2-{\tilde{\varphi }}_1{L}^{12}+{\varphi }_1{\tilde{\varphi }}_1{L}^{13}+{\varphi }_2{\tilde{\varphi }}_1{L}^{14})y_t^{\ast }=c+{ϵ}_t$$
$=>$
$$y_t^{\ast }=c+{\varphi }_1{y}_{t-1}^{\ast }+{\varphi }_2{y}_{t-2}^{\ast }+{\tilde{\varphi }}_1{y}_{t-12}^{\ast }-{\varphi }_1{\tilde{\varphi }}_1{y}_{t-13}^{\ast }-{\varphi }_2{\tilde{\varphi }}_1{y}_{t-14}^{\ast }+{ϵ}_t$$

   我们可以看到，最后的结果也类似于加性模型，但是系数实际上是非季节性和季节性参数的乘性组合。

SARIMAX

   虽然本篇文章旨在引出SARIMAX，但其基本原理与ARIMAX类同，只是多了季节性因素而已。为了简单介绍引入外部变量的原理和方法，我们在理论部分着重介绍ARMAX，先忽略季节因素和差分（假设时序平稳)。
   模型名称中的’X’，主要是指外部变量。也就是说，你需要实时知道 $t$ 时刻的外部变量的值 $x_t$。
   ARMAX模型最简单的方案就是把外部变量 $x_t$ 当成一个依赖变量，加入到ARMA的回归方程中，如下:
$$y_t=\beta x_t+{\varphi }_1{y}_{t-1}+\cdots +{\varphi }_p{y}_{t-p}-{\theta }_1{z}_{t-1}-\cdots -{\theta }_q{z}_{t-q}-z_t$$
这里的减号其实并没有关系，只是为了后续推导便利。
   $x_t$ 就是外部变量，$\beta$ 为其回归系数。这种方式看上去很直接，但是最大的缺点是回归系数 $\beta$ 很难解释。 $\beta$ 的值不能解释为 $x_t$ 增加1，$y_t$ 所增加的值。回归方程中的滞后项表示 $\beta$ 只能被解释为条件于响应变量前面的值，这非常不直观。

   我们使用滞后算子来改写一下，
$$\varphi (L)y_t=\beta x_t+\theta (L)z_t\text{or}{y}_t=\frac{\beta }{\varphi (L)}{x}_t+\frac{\theta (L)}{\varphi (L)}{z}_t,$$
这里，$\varphi (L)=1-{\varphi }_{1}L-\cdots -{\varphi }_p{L}^p$，且$\theta (L)=1-{\theta }_{1}L-\cdots -{\theta }_q{L}^q$

注意到自回归系数(AR)混入了外部变量系数和误差项系数。
我们可以考虑回归ARMA模型的误差，
$$y_t=\beta x_t+n_t$$
$$n_t={\varphi }_1{n}_{t-1}+\cdots +{\varphi }_p{n}_{t-p}-{\theta }_1{z}_{t-1}-\cdots -{\theta }_q{z}_{t-q}+z_t$$
这种形势下，回归系数就可以以其通用解释方法来解释。以滞后算子表示如下，
$$y_t=\beta x_t+\frac{\theta (L)}{\varphi (L)}{z}_t.$$

Box和Jenkis带来一个使用转换函数的泛化模型，
$$y_t=\frac{\beta (L)}{v(L)}{x}_t+\frac{\theta (L)}{\varphi (L)}{z}_t.$$
这种形式允许外部变量或称为协变量的作用恨到（通过$\beta (L)$，或者使协变量的衰减效果很大(通过$V(L)$)

有时我们将这种回归方式称为动态回归模型(dynamic regression models).
对于非平稳时序，我们只需要将公式中的 $\varphi (B)$ 替换成 ${\nabla }^d\varphi (B)$，这里 $\nabla =(1-B)$，模型就转变成了ARIMAX; 最后我们只需要增加季节性因子就变成了SARIMAX，模型方程如下：
$$y_t={\beta }_t{x}_t+u_t$$
$${\varphi }_p(L){\tilde{\varphi }}_P(L^s){\Delta }^d{\Delta }_s^D{u}_t=A(t)+{\theta }_q(L){\tilde{\theta }}_Q(L^s){ϵ}_t$$

未完，待续…