大数据分析笔记 (7) - 时间序列分析(Time Series Analysis)_stationary time series

总览

时间序列分析试图模型化一段时间内观测到的数据的底层结构。一个时间序列 (表示为 Y = a + bX) 是一个在时间上具有相同间隔值的有序序列。如图，图中提供了在12年内每个月的航空乘客的数量。

该方法可以应用在金融，经济，生物，工程，零售，制造等领域。
目标在于：

• 对时间序列的结构进行识别和建模 (identify and model)。
• 预测时间序列中未来的值。

一个时间序列由在时间上具有相同间隔值的有序序列组成。一个时间序列可以由以下部分组成: 趋势(Trend)，季节效应(Seasonality)，周期(Cyclic)和随机性(Random)。

趋势： 指在时间序列中的长期移动 (long-term movement)。其表示观测值是否随着时间变化增加或者降低。

季节效应： 描述了观测数据在时间上的固定的周期性波动 (fixed, periodic fluctuation)。

周期： 可以视为周期性，但是不是固定的波动(periodic but not fixed)。比如，零售时间序列可能常常遵循经济繁荣与萧条 (boom-bust cycles of the economy) 的周期。

随机性： 除了上述三部分剩下的部分。噪声 + 底层结构 可以被建模，以预测时间序列的未来值。

Box-Jenkins方法

1. 治理数据 (condition data) 与选择模型 (select a model)
   • 识别和考虑时间序列中的任何趋势 (trend) 与季节效应 (seasonality)。
   • 评估其余时间序列并确定一个合适的模型。
2. 估计模型参数。
3. 评估模型，如果有必要，返回到第一步。

ARIMA (自回归求和移动平均模型)

ARIMA: AutoRegressive Integrated Moving Average.

要应用ARIMA，我们需要把趋势和季节效应从时间序列中移除。这种时间序列叫做平稳时间序列 (stationary time series)。

若一个时间序列yt = {t = 1, 2, 3, …}是一个平稳时间序列，则：

• yt 的 期望值 (均值 mean) 对于所有的 t 都是常数 (constant)。
• yt 的 方差 (variance) 是有限的。
• 对于所有的t，yt 和 y t+h 的协方差值 (covariance) 只依赖于 h = 0,1,2,…。如 cov(3) = cov(y1, y4) = cov(y2, y5)

协方差可以用于衡量两个变量是如何一起变化的。如果两个变量相互独立，则它们的协方差是0。如果两个变量向着同一个方向变化，那么它们的协方差是正的。反之，如果两个变量向着相反方向变化，那么它们的协方差就是负的。
值得注意的是，对于 h=0，所有 t 的 cov(0) = cov(yt, yt) = var(yt)。
平稳时间序列例子(下图)：外观平整 flat looking(无趋势)，恒定方差 constant variance (类似散点 scattering)。但是，绘图无法洞察时间序列中变量的协方差及其底层结构。

自相关函数 (ACF)

ACF: Autocorrelation function
对于平稳时间序列，ACF的定义如下：

由于cov(yt, yt)是方差，ACF与两个变量的相关函数corr(yt, y t+h)类似，其值介于 -1 ~ 1 之间。因此，ACF(h) 的绝对值越接近1，那么 yt 就越能作为 y+h 的有效预测 (predictor)。

ACF中的h量被称为Lag，表示时间点 t 和 t + h 的差异。
当lag是0的时候，ACF提供每个点于其本身的相关性。所以ACF(0)总是等于1。
根据上面的ACF图，
当 lag=1 时，yt 和 y t -1的相关性大约是0.9，非常接近于1。因此y t-1 对 yt 似乎是一个不错的预测。
ACF(2) ≈ 0.8，所以 y t-2 对 yt 似乎也是个不错的模型。
自此，我们可以考虑这样一个模型，该模型将 yt 表述为 前8个项的线性和(因为 ACF(0 ≤ n ≤ 8) ＞0.6)。这类模型被称为 8阶自回归模型 (autoregressive model of order 8)。

自回归模型

对于一个平稳时间序列yt，t = 1, 2, 3,…，表示为AR( p )的p阶自回归模型如下：

δ 是一个以非零为中心的时间序列的常量。
Ф 是 j = 1, 2, 3, …, p 时的一个常量。
y t-j 是在时间t-j时的时间序列的值。
Фp ≠ 0
对于所有 t 值， ε 服从正态分布。
因此，时间序列中的一个特定的点可以被表述为时间序列中前p个值 y t-j (其中 j = 1, 2, …, p)的线性组合(linear combination)，外加一个随机的误差项 ε。在这个定义中，时间序列中的 ε 通常称为白噪声过程(white noise process)，用来表示时间序列中的随机，独立波动的部分。

如何在AR( p )中决定p的值？
对于AR(1)的模型，以 δ=0 为中心，我们有

基于上述公式，可知

替换 y t-1 之后可得

这个过程可以被不断重复，y t-h, h = 1, 2, 3…。
因此，即使在简单的 AR(1) 模型中，与较大的lag也存在着显著的自相关(significant autocorrelation)。

我们需要测量 yt 和 y t+h (h = 1, 2, 3, …) 之间的自相关系数，并在测量中除去 y t-1 到 y t-h-1 的影响。

部分自相关函数 (Partial autocorrelation function - PACF)

当AR( p )过程的PACF变为0时，这可以指示AR模型的阶数(order of the AR model)。因此，让我们检查数据的PACF，并找到PACF变为0之后的lag值。
也就是说，在使用线性回归移除 yt 与 y t+h 之间的变量对 yt 和 y t+1 的影响后，PACF就是剩余的相关系数。


观测：当 p ≥ 2 时，PACF变为 0 (急剧下降)，因此表明AR模型的阶为2 (非8)。这说明在移除了y t+1 和 y t+h+1 值的影响之后，y t+1 和 y t+h+1 的部分相关性也就相对小了。

因为 ACF 和 PACF 是基于相关性的，因此正值负值都有可能。对于各种lag值，应该考虑函数大小的绝对值。

移动平均模型 (Moving Average Models - MA)

对于以 0 居中的时间序列 t，表示为 MA(q) 的 q阶移动平均模型公式如下：

在一个MA(q)模型中，时间序列的值是当前白噪声项 (white noise) 和先前 q 个白噪声项的线性组合。所以先前的随机冲击 (random shock) 直接影响到时间序列的当前值。
和AR( p ) 模型的两个不同：

• ε t-1 被直接传播到 yt。
• ε t 之影响当前和未来的 p 值。

有限的MA模型始终是平稳时间序列。如下图的 MA(3) 时间序列

该模型的ACF

在一个自回归模型中，ACF缓慢衰减 (slowly decays)。但是对于MA(3)模型，ACF在 lag = 3 以后突然截断。因此，ACF可以帮助识别 MA(q) 中的阶q。
原因在于：

由图可知，yt 表达式 和 y t-1 到 y t-3 表达式共享了特定的白噪声变量，因此这三个变量都与 yt 相关联。然而，yt 表达式 和 y t-4 没有共享 (have no overlapping) 白噪声变量。所以 yt 和 y t-4 之间的相关系数理论 (theoretical correlation) 上是0。

自回归移动平均模型 (ARMA - Autoregression Moving Average Model)

AR( p ) 和 MA(q) 经常被结合成自回归移动平均模型 ARMA(p, q)。

为了合理地应用ARMA模型，时间序列必须是一个平稳时间序列。然而，许多时间序列随时间呈现一定的趋势。对于非平稳时间序列，我们需要调整数据来移除这种趋势。

• 一种可能的转换是对时间序列执行回归分析(linear or higher-order regression model)，然后从每个观察到的y值中减去拟合回归线的值(remove trends)。
• 第二种选项是计算连续y值之间的差异 (difference between successive y-values)，这就是所谓的差分法 (differencing)。
  
• 如果仍然不平稳，则对时间序列再次进行差分。
  
  

ARIMA Model

自回归求和移动平均模型。这是将差分操作包含(集成)在了ARMA模型中。

ARIMA(p,d,q): 在应用了 d 次差分之后，ARMA(p, q) 模型才能应用到时间序列 yt 上。

通常还需要考虑时间序列中的季节效应模式 (seasonal pattern)。可选方案是季节自回归求和移动平均模型(Seasonal ARIMA model):
ARIMA(p, d, q) x (P, D, Q)s

其中，s表示季节时段，P是AR模型中项在s时段上的数量。D是在s时段上的差异；Q是MA模型中项在s时段上的数量。(s的典型取值：52是周数据，12是月数据，7是日数据)

一次时间序列分析示例

获取一个原时间序列(每月汽油产量 - 百万桶)。我们需要预测短期产量。

对汽油产量时间序列进行了一次差分后，结果如下图

对差分以后的时间序列的ACF和PACF图如下所示。(虚线指95%置信度时的上下边界，位于边界外的ACF或PACF的任何值都表示该值与 0 显著不同)

可以从图中看到，ACF图在lag为12，24，36，48时缓慢衰减是我们感兴趣的(代表了季节性)。这表明了一个每12个月的季节效应自回归模型。在PACF图中也能看到，lag为12时，PACF值很大，而lag等于24，36，48时又接近于0。因此我们会考虑一个周期等于12的季节 AR(1) 模型。
ARIMA(0, 1, 0) x (1, 0, 0)12

上图中ACF图表明在lag12，24，36，48时的自回归行为已经被季节效应 AR(1) 解决。唯一剩余的ACF值发生在当 lag1 的时候。

如上图，PACF是缓慢衰减的，而ACF却在lag1处陡然削减。因此对于差分序列上的ARMA模型的非季节部分，应该考虑MA(1)模型。

使用ARIMA(0, 1, 1) x (1, 0, 0)12


新的ARIMA模型分析的ACF和PACF图表明，不需要在ARIMA模型中考虑更多的项。

建立并验证ARIMA模型

1. 比较拟合的时间序列模型：
   模型中的系数通过 最大似然估计(Maximum Likelihood Estimation MLE) 来估计。R语言提供了几个度量 (AIC, AICc, BIC) 用来比较两个拟合模型的适当性。首选的模型应该具有最小的AIC, AICc, BIC值。
   
2. 正态性 (normality) 和常数方差 (constant variance)
   最后的模型检验步骤是检验残差的正态性假设。如图是一个均值近似为0，方差为常数的残差。
   
   下图的直方图和Q-Q图支持了“误差项是正态分布”的假设。
   
   
   如果正态性或者常数方差的假设不显示为真，那么有必要在拟合ARIMA模型前转换时间序列。一种常见的转换就是应用对数函数 (logarithm function)。

注意事项

• ARIMA模型的优势是，它可以仅仅基于历史事件序列数据 (historical data) 来分析感兴趣的变量。
• ARIMA建模通常会忽略任何额外的输入变量，因而其简化了预测过程 (forecasting process)。如果使用的是回归模型，能考虑其他输入变量的影响，但是同样的，当我们想要预测的时候也会要求有这些输入变量的预测值。
• 缺点：模型无法指明哪些潜在变量 (underlying variables) 对结果产生了影响。比如使用ARIMA建模预测未来的零售销售，拟合后的 ARIMA 模型无法指明需要做些什么才能增加销售量。
• 注意：使用时间序列分析时，需要警惕严重冲击 (severe shocks) 对系统的影响。同时，时间序列分析只能在短期预测 (short-term forecast) 中使用。