数据科学：时间序列（一）—— ARIMA

自回归差分移动平均(Autoregressive Integrated Moving Average, ARIMA)模型。

原理

平稳过程

平稳过程是一种特殊的随机过程，可分为强平稳过程和弱平稳过程：

强平稳过程（Strict(ly) stationary process）：时间序列中给定长度的任意两段子序列联合概率分布相同。可以推知强平稳过程 $x_t$ 的期望和方差不随t变化；
弱平稳过程：时间序列满足均值平稳性（stationary in mean）和二阶平稳性（secondary order stationary）
- 均值平稳性：均值不变，与t无关
- 二阶平稳性：$X(t)$ 与 $X(t-k)$ 的协方差只与滞后阶数k有关，与t无关，即相关系数只取决于时间间隔而与时间起点无关

到目前为止，时间序列分析基本上是以平稳时间序列为基础的，对于非平稳时间序列的统计分析，其方法和理论都很有局限性。强平稳过程在实际中几乎不可能被满足，通常情况下时间序列分析讨论的是弱平稳过程。但是非平稳序列一般都可以通过对数和差分方法转化为平稳序列，消除（或减小）时间序列的趋势和季节性。直观上看，平稳序列的时序图应该始终在一个常数值附近随机波动，而且波动的范围有明显的相似性特点。

自相关性

定义： $rt$ 和 $r{t-k}$ 的相关系数称为 $r_t$ 的间隔为 $k$ 的自相关系数。

在弱平稳假设下，这个间隔为 k 的自相关系数与时间 t 无关，而仅仅与间隔 k 有关，由 $\rho_k $ 表示。由相关系数的定义可知：

ACF 自相关图：X(t) 时间序列与 X(t-1), X(t-2),.., X(t-k)各阶之后时间序列之间的相关系数与滞后阶数 k 的关系。

如果 k 阶滞后序列 X(t-k)与原始序列 X(t)的相关系数在 [-0.2, 0.2] 之间，则称 X(t-k)与 X(t)不具有相关性。如果一个时间序列的任何一阶滞后序列与原始序列都不具有相关性，则该时间序列不具有自相关性
白噪声序列：均值为 0，且不具有自相关性的平稳序列

AR

p 阶自回归模型 $AR(p)$可以表示如下：

$y_{t} = c + \phi_{1}y_{t-1} + \phi_{2}y_{t-2} + \dots + \phi_{p}y_{t-p} + \varepsilon_{t}$

其中：

$c$：常数项；
$y_{t-k}$：前k期历史值；
${\varepsilon _{t}}$：白噪声；

对于一个 $AR(1)$ 模型而言：

当 $ϕ_1=0$ 时, $y_t$ 相当于白噪声；
当 $ϕ_1=1$ 并且 $c=0$ 时, $y_t$ 相当于随机游走模型；
当 $ϕ_1=1$ 并且 $c≠0$ 时, $y_t$ 相当于带漂移的随机游走模型；
当 $ϕ_1<0$ 时, $y_t$ 倾向于在正负值之间上下浮动;

我们通常将自回归模型的应用限制在平稳数据上，并且对回归系数也施加一些约束条件：

对于 $AR(1)$ 模型：$−1<ϕ1<1$；
对于 $AR(2)$ 模型：$−1<ϕ2<1$，$ϕ_1+ϕ_2<1$，$ϕ_2−ϕ_1<1$ 。
当 $p≥3$ 时，约束条件会更为复杂一些。R在估计自回归模型时可以解决这个问题

优点与限制：

必须具有自相关，自相关系数（$\varphi _{i}$）是关键。如果自相关系数(R)小于0.5，则不宜采用，否则预测结果极不准确；
自回归只能适用于预测与自身前期相关的经济现象，即受自身历史因素影响较大的经济现象，如矿的开采量，各种自然资源产量等；对于受社会因素影响较大的经济现象，不宜采用自回归，而应改采可纳入其他变数的向量自回归模型；

MA

不同于使用预测变量的历史值来进行回归，移动平均模型（moving average model）使用历史预测误差来建立一个类似回归的模型。q 阶移动平均模模型 $MA(q)$ 可以表示如下：

$y_{t} = c + \varepsilon_t + \theta_{1}\varepsilon_{t-1} + \theta_{2}\varepsilon_{t-2} + \dots + \theta_{q}\varepsilon_{t-q},$

其中：

$c$：常量；
$\varepsilon_t$ ：白噪声;

任何一个 $AR(p)$ 模型其实都是可以用一个 $MA(∞)$模型来表示：

$\begin{align*} y_t &= \phi_1y_{t-1} + \varepsilon_t\\ &= \phi_1(\phi_1y_{t-2} + \varepsilon_{t-1}) + \varepsilon_t\\ &= \phi_1^2y_{t-2} + \phi_1 \varepsilon_{t-1} + \varepsilon_t\\ &= \phi_1^3y_{t-3} + \phi_1^2\varepsilon_{t-2} + \phi_1 \varepsilon_{t-1} + \varepsilon_t\\ &\text{etc.} \end{align*}$

假定中 $−1<ϕ_1<1$ ， $ϕ^k_1$ 的值会随着 k 的增大而减小，最终可以得到下式：

$y_t = \varepsilon_t + \phi_1 \varepsilon_{t-1} + \phi_1^2 \varepsilon_{t-2} + \phi_1^3 \varepsilon_{t-3} + \cdots,$

如果我们对MA模型的系数加以限制，那么我们也可以反过来用AR模型来表示MA模型。这样的MA模型被称为可逆的(invertible)。也就是说，我们可以将任何一个可逆的 $MA(q)$ 模型表示为一个 $AR(∞)$ 模型。可逆模型不仅仅是让我们足以将任何MA模型转化为AR模型，它们还拥有很多非常棒的数学性质。

其他模型关于可逆性的限制和平稳性的限制类似：

对于MA(1)模型: $−1<θ_1<1$
对于MA(2)模型: $−1<θ_2<1$，$θ_2+θ_1>−1$，$θ_1−θ_2<1$
对于 $q≥3$的其他更复杂的情况，R会在在估计模型时解决这个问题

差分

延迟算子 B 是一个重要的标记，它被用于表示时间序列的延迟，当 B 被用于 $y_t$ 时，意味着将时间反向回溯一个单位时段：

$B^k y_{t} = y_{t - k}$

$y_t$ 的 d 阶差分可以表示为：

$(1 - B)^{d} y_{t}$

ARIMA

当我们将差分和自回归模型以及移动平均模型结合起来的时候，我们可以得到一个非季节性 ARIMA 模型。ARIMA 是 AutoRegressive Integrated Moving Average 的简称。（在这里Integrated指的是差分的逆过程) ARIMA模型的表示如下：

$\begin{equation} y'_{t} = c + \phi_{1}y'_{t-1} + \cdots + \phi_{p}y'_{t-p} + \theta_{1}\varepsilon_{t-1} + \cdots + \theta_{q}\varepsilon_{t-q} + \varepsilon_{t} \end{equation}$

上式中 $y′_t$ 是差分序列（它可能经过多次差分）。右侧的“预测变量”包括 $y_t$ 的延迟值和延迟的误差。我们将这个模型称为ARIMA( p,d,q ) 模型，参数含义如下：

p = 自回归模型阶数
d = 差分阶数
q = 移动平均模型阶数

我们之前讨论的很多模型其实都是ARIMA模型的特殊情况，如下表所示：

类型	参数
白噪声	ARIMA(0,0,0)
随机游走模型	ARIMA(0,1,0) with no constant
带漂移的随机游走模型	ARIMA(0,1,0) with a constant
自回归模型	ARIMA(p,0,0)
移动平均模型	ARIMA(0,0,q)

一旦我们开始组合不同模型来形成复杂的模型，延迟算子就会显得格外简便。ARIMA( p,d,q ) 可以被表示为:

$\begin{equation} \begin{array}{c c c c} (1-\phi_1B - \cdots - \phi_p B^p) & (1-B)^d y_{t} &= &c + (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t\\ {\uparrow} & {\uparrow} & &{\uparrow}\\ \text{AR(p)} & \text{d differences} & & \text{MA(q)}\\ \end{array} \end{equation}$

实践

ARIMA 模型使用的一般流程：

根据时间序列的散点图、自相关函数和偏自相关函数图识别其平稳性。
对非平稳的时间序列数据进行平稳化处理。直到处理后的自相关函数和偏自相关函数的数值非显著非零。
根据所识别出来的特征建立相应的时间序列模型。平稳化处理后，若偏自相关函数是截尾的，而自相关函数是拖尾的，则建立AR模型；若偏自相关函数是拖尾的，而自相关函数是截尾的，则建立MA模型；若偏自相关函数和自相关函数均是拖尾的，则序列适合ARMA模型。
参数估计，检验是否具有统计意义。
假设检验，判断（诊断）残差序列是否为白噪声序列。
利用已通过检验的模型进行预测。

时间序列预处理

平稳性检验

时序图检验：根据平稳时间序列的均值和方差都为常数的性质，平稳序列的时序图显示该序列值始终在一个常数附近随机波动，而且波动的范围有界；如果有明显的趋势性或者周期性，那它通常不是平稳序列
自相关图检验：平稳序列具有短期相关性，这个性质表明对平稳序列而言通常只有近期的序列值对现时值得影响比较明显，间隔越远的过去值对现时值得影响越小。随着延迟期数k的增加，平稳序列的自相关系数会比较快的衰减趋向于零，并在零附近随机波动，而非平稳序列的自相关系数衰减的速度比较慢。
单位根检验：最常用，原理待补充
白噪声检验：一般是构造检验统计量来检验序列的纯随机性，常用的检验统计量有Q统计量、LB统计量，由样本各延迟期数的自相关系数可以计算得到检验统计量，然后计算出对应的p值，如果p值显著大于显著性水平a，则表示该序列不能拒绝纯随机的原假设，可以停止对该序列的分析

平稳化处理

取对数
取差分

参数选择

自相关图和偏自相关图

如果数据来自于 ARIMA(p, d, 0)或者ARIMA(0, d, q)模型，则自相关图和偏自相关图在判定 p 或者 q 的取值时非常有用。如果 p 和 q 都为正，则这些图在寻找最合适的 p 和 q 值时不再有用。

截尾：自相关函数（ACF）或偏自相关函数（PACF）在某阶后在0附近随机波动
拖尾：始终有非零取值，不会在k大于某个常数后就恒等于零

可以通过自相关和偏自相关图的截尾还是拖尾来选择模型和阶数：

AR模型：自相关系数拖尾，偏自相关系数截尾；
MA模型：自相关系数截尾，偏自相关函数拖尾；
ARMA模型：自相关函数和偏自相关函数均拖尾。

如果你的图片成这样，估计十有八九是一个ARMA模型了。自相关7阶拖尾（n从7开始缩至置信区间），偏自相关2阶拖尾。

极大似然估计

AIC信息准则（Akaike Information Criterion）在选择用于回归模型的变量时非常有用，同样在确定ARIMA模型阶数时也可以发挥很大作用。它可以写作:

$\text{AIC} = -2 \log(L) + 2(p+q+k+1),$

这里的 L 是数据的似然函数,当 $c≠0$ 时 $k=1$ 且当 $c=0$ 时 $k=0$。需要注意到是式子中的最后一项是模型中参数的个数（包括残差方差 $σ2$ ）

对于ARIMA模型而言，修正过的AIC值可以被表示为：

$\text{AICc} = \text{AIC} + \frac{2(p+q+k+1)(p+q+k+2)}{T-p-q-k-2},$

BIC信息准则（Bayesian Information Criterion）的表示方式如下:

$\text{BIC} = \text{AIC} + [\log(T)-2](p+q+k+1).$

通过最小化AIC，AICc或者BIC我们都可以得到最优模型。其中我们倾向于使用AICc。有一点需要格外注意：这些信息准则在选择模型的合适的差分阶数（ d ）时效果并不好，只能被用于选择 p 和 q 的取值。这是因为差分改变了似然函数计算所使用数据，这会使得不同差分阶数的模型的AIC值无法比较。因此我们需要用其他方法先选择出合适的 d ，然后再通过 AICc 来选择 p 和 q 的取值。

参数选择

如果您有一个SARIMAX（p，d，q）x（p，d，q，s）模型，快速的封套计算表明您至少需要：

d+d*s+max（3*q+1，3*q*s+1，p，p*s）+1个观测值

（0，0，0）x（1，0，0，52）=>0+0+52+1=53

（0，0，0）x（1，1，0，52）=>0+52+52+1=105

（0，0，0）x（2，1，0，52）=>0+52+104+1=157

以及：

（0，0，0）x（0，0，1，52）=>0+0+3*52+1+1=158

（0，0，0）x（0，1，1，52）=>0+52+3*52+1+1=210

（0，0，0）x（1，1，1，52）=>0+52+3*51+1+1=210

哈哈哈