时间序列可以被预测，主要基于以下事实：我们可以部分掌握影响该时间序列的因素的变化情况。换句话说，对时间序列进行预测，其实就是利用各种理论和工具，对观察到的时间序列进行“抽丝剥茧”，以试图掌握其变化的本质，从而对未来的表现进行预测。

谈论时间序列预测，就绕不开ARIMA模型；而应用ARIMA模型的关键，在于确定参数$p$、$d$和$q$的值（此处不考虑季节性因素的影响）。我整理了一下网上对于如何解决这个问题的有关资料，总结在此。

本篇为第一篇，主要分析ACF与PACF。

1、相关性与相关系数

在统计学中，相关性可以用来衡量两个随机变量之间的线性关系的强度和方向。例如，从经验上可以知道，空气温度和居民用电量之间存在着相关性。

为了能够量化描述这种相关性，人们发明了相关系数。在众多相关系数中，皮尔逊相关系数（Pearson correlation coefficient）是最常用的一个，它是两个变量的协方差与其标准差的乘积之比，取值区间为$[-1,1]$。

本文后续中提到的相关系数，均指皮尔逊相关系数。

2、ACF

从字面上的意思理解，自相关函数（Autocorrelation Function）就是用来计算时间序列自身的相关性的函数。

对于一个时间序列$x_t$，根据相关系数的计算公式易得：$corr(x_t,x_t)=1$。这也很好理解，两个完全相同的序列的相关性当然为1。

实际上，在应用自相关函数时，其输入分别为原始的时间序列$x_t$及其$k$阶滞后序列$x_{t-k}$。具体地，假设原始序列为 { x 1 , x 2 , . . . , x n } \{x_1,x_2,...,x_n\} {x1​,x2​,...,xn​}，那么其1阶滞后序列为 { x 1 , x 2 , . . . , x n − 1 } \{x_1,x_2,...,x_{n-1}\} {x1​,x2​,...,xn−1​}，于是，$corr(x_t,x_{t-1})$就变成了
c o r r ( x t , x t − 1 ) = c o r r ( { x 2 , . . . , x n } , { x 1 , x 2 , . . . , x n − 1 } ) corr(x_t,x_{t-1})=corr(\{x_2,...,x_{n}\},\{x_1,x_2,...,x_{n-1}\}) corr(xt​,xt−1​)=corr({x2​,...,xn​},{x1​,x2​,...,xn−1​})
注意，为了使两个序列的长度一致，对原始序列进行了截取。

可以想见，如果将滞后数$k$不断增加，那么就会得到一系列对应的相关系数，将这些相关系数绘制成图，就得到了一个时间序列的自相关系数图。

下图展示了一个价格时间序列，所需数据可以从这里下载。

用如下代码，可以绘制该时间序列的ACF图：

import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacfdf = pd.read_excel('coal-price-mixed.xlsx')
plot_acf(df['price'], lags=45, auto_ylims=True)
plt.show()

得到的图像如下：

在上图中，横坐标表示滞后的阶数，纵坐标表示对应的滞后序列与原始序列的相关系数。可以看出，随着滞后阶数的增加，滞后序列与原始序列的相关性也在不断地降低。图中的蓝色区域表示置信区间，它被用来标识相关系数是否具有统计显著性。简单来说，如果相关系数落在置信区间内，则表明对应的两个序列的相关系数并不能代表其真实相关性。

即使是两个完全不相干的白噪声序列，由于随机性的影响，其相关系数也不可能全都为0，因此，需要使用置信区间来过滤掉那些由于随机性造成的“伪相关”。关于如何确定置信区间的内容，则不在本文讨论的范围内。

3、自回归时间序列与PACF

3.1 AR模型

前面已经提到，要对时间序列进行预测，就要从该时间序列的历史变化中找到影响未来值的因素。“过去影响未来”，用公式来说明就是$x_t={\alpha }_0+{\sum }_{i<t}{\alpha }_i{x}_i$。具体地，$x_t={\alpha }_0+{\alpha }_1{x}_{t-1}$表示过去一个时间步的值会影响当前的值，${\alpha }_0$是截距项，${\alpha }_1$则表示$x_{t-1}$对当前值影响的程度。

这类似于$y=kx+b$的一元线性模型，但需要注意前提条件：$x_t$只受其前一步的影响，而更长时间步之前的诸如$x_{t-2}$、$x_{t-3}$等则不会对$x_t$产生影响（或者说产生的影响很小以至于可以被忽略）。

如果更长时间步之前的值，例如$x_{t-2}$，也会影响当前值，那么则需要修改模型为：$x_t={\alpha }_0+{\alpha }_1{x}_{t-1}+{\alpha }_2{x}_{t-2}$。类似地，还可以扩展到时间步为$p$的情况：
$$x_t={\alpha }_0+\sum _{j=1}^p{\alpha }_j{x}_{t-j}$$
这就是$AR(p)$模型。

在应用$AR(p)$时，需要时刻注意其成立的前提条件：当前时间步的值仅受过去$p$个时间步的值的影响。也就是说，过去$p$个时间步的值所包含的信息就（在很大程度上）决定了当前时间步的值。

接下来，一个自然而然的问题就是：如何确定$p$的值？

3.2 确定阶数的关键

通过上一节的分析，可以得出如下结论：确定$AR$模型中$p$的值，其实就是在判断至少需要几个时间步的值就可以预测当前值。

为了能够进行判断，接下来还需要一个判断的标准，用于说明为什么选用过去$p$个时间步的值就可以进行预测而不是$p-1$或$p+1$。

以$AR(2)$为例，我们认为$x_{t-1}$和$x_{t-2}$含有决定$x_t$的值的信息，而$x_{t-3}$及更之前的值所含的信息则可以忽略不计。如果用相关性的概念来描述这个例子，则上一句话就变成了“$x_{t-1}$和$x_{t-2}$为与$x_t$的相关性最大的两个，且足够用于预测$x_t$”。

巧了，上文讨论的ACF刚好就是用于计算原始序列和其滞后序列之间相关性的，那相关系数在置信区间之外的前几个序列是否就对应着此处的$p$呢？

答案是否定的。

举个例子来说明为什么不能直接这样做：假设时间序列确由一个$AR(1)$模型生成的，也就是说，$x_{t-1}$影响$x_t$，$x_{t-2}$影响$x_{t-1}$，如果按照ACF的计算方法，那么$x_{t-2}$很可能会通过$x_{t-1}$来“影响”$x_t$。这就会导致无法判断$x_{t-1}$中是否已经包含了足够的信息来预测$x_t$，还是说必须再加上$x_{t-2}$。

也就是说，按照ACF计算得到的$x_t$与$x_{t-2}$的相关性并非是“独立的”——因为$x_{t-1}$的存在，导致ACF无法计算$x_{t-2}$到底给预测$x_t$带来了多少有用的信息。

因此，为了计算$x_{t-p}$给$x_t$带来了多少“纯粹”的影响，需要发明一种新的计算相关系数的方法。

3.3 PACF

PACF，指偏自相关函数（Partial Autocorrelation Function），是一种计算$x_t$和$x_{t-k}$的相关系数时可以移除$x_{t-1}$到$x_{t-(k-1)}$影响的方法。

仍以$AR(2)$为例，来看PACF是如何计算的。用$\varphi (2,2)$来表示$x_t$和$x_{t-2}$的偏自相关系数，则其计算公式如下：
$$\varphi (2,2)=corr(x_t-{\hat{x}}_t,x_{t-2}-{\hat{x}}_{t-2})$$
其中，${\hat{x}}_t$（${\hat{x}}_{t-2}$）表示通过最小化均方误差建立$x_t$（$x_{t-2}$）对$x_{t-1}$的回归方程后所做出的预测值。那么，$\varphi (2,2)$实际上计算的就是两个回归方程的残差序列的相关系数。

我尝试理解了一下这种计算方式的含义，写在这里，仅供参考：

$x_t-{\hat{x}}_t$是用$x_{t-1}$来预测$x_t$时的残差，可以表示$x_t$中无法被$x_{t-1}$解释的那部分的信息量；类似地，$x_{t-2}-{\hat{x}}_{t-2}$表示$x_{t-2}$中无法被$x_{t-1}$解释的那部分的信息量。如果它们的相关系数很大，即$x_t$中除去$x_{t-1}$影响后的信息与$x_{t-2}$中除去$x_{t-1}$影响后的信息具有很高的相关性，从而也就意味着，$x_{t-2}$中除去$x_{t-1}$影响后的信息与$x_t$中除去$x_{t-1}$影响后的信息具有高相关性（$x_{t-2}$对$x_t$的“纯粹”的影响很大），据此可以推断，$x_{t-2}$中含有决定$x_t$值的信息。

读者可能会有疑问：$x_{t-2}$比$x_{t-1}$发生的时间要早，那用后者来预测前者的实际意义是什么呢？这个问题表明，PACF是一种利用数学技术来计算两个序列“纯粹”相关性的方法，并非整个过程都对应着直观的现实理解。

理解了上述关于2阶PACF的计算，不难拓展到$p$阶的情况，这里直接给出公式：
$$\begin{gathered} {\varphi }_{1,1}=corr(x_t,x_{t-1}),forp=1, \\ {\varphi }_{p,p}=corr(x_t-{\hat{x}}_t,x_{t-p}-{\hat{x}}_{t-p}),forp\ge 2 \end{gathered}$$
其中，${\hat{x}}_t$为$x_t$对 { x t − 1 , x t − 2 , . . . , x t − ( p − 1 ) } \{x_{t-1},x_{t-2},...,x_{t-(p-1)}\} {xt−1​,xt−2​,...,xt−(p−1)​}做回归后的预测值，${\hat{x}}_{t-p}$为$x_{t-p}$对 { x t − 1 , x t − 2 , . . . , x t − ( p − 1 ) } \{x_{t-1},x_{t-2},...,x_{t-(p-1)}\} {xt−1​,xt−2​,...,xt−(p−1)​}做回归后的预测值。回归方程的参数均由最小化均方误差确定。

3.4 绘制PACF图

采用与2中相同的数据集，可以通过如下代码绘制该时间序列的PACF图：

plot_pacf(df['price'], lags=45, auto_ylims=True, method='ols')
plt.show()

得到的图像如下：

3.5 其他

可以通过以下代码来逐步计算PACF值：

from sklearn.linear_model import LinearRegressiondf['price-shift-1'] = df['price'].shift(1)
df['price-shift-2'] = df['price'].shift(2)lr = LinearRegression()
lr.fit(df[['price-shift-1']].iloc[1:, :], df[['price']].iloc[1:, :])
res1 = df[['price']].iloc[1:, :] - lr.predict(df[['price-shift-1']].iloc[1:, :])  # 残差序列1lr.fit(df[['price-shift-1']].iloc[2:, :], df[['price-shift-2']].iloc[2:, :])
res2 = df[['price-shift-2']].iloc[2:, :] - lr.predict(df[['price-shift-1']].iloc[2:, :])  # 残差序列2print(np.corrcoef(res1['price'].values[1:], res2['price-shift-2'].values)[1, 0])
# -0.9028955361943247

将上述结果与pacf的结果进行对比：

from statsmodels.tsa.api import pacf
print(pacf(df['price'], method='ols-inefficient'))[2]  # p为2时的偏自相关系数
# -0.9029804339992743

参考

• https://timeseriesreasoning.com/contents/partial-auto-correlation
• Partial autocorrelation function