时间序列可以被预测，主要基于以下事实：我们可以部分掌握影响该时间序列的因素的变化情况。换句话说，对时间序列进行预测，其实就是利用各种理论和工具，对观察到的时间序列进行“抽丝剥茧”，以试图掌握其变化的本质，从而对未来的表现进行预测。
谈论时间序列预测，就绕不开ARIMA模型；而应用ARIMA模型的关键，在于确定参数 p p p、 d d d和 q q q的值（此处不考虑季节性因素的影响）。我整理了一下网上对于如何解决这个问题的有关资料，总结在此。
本篇为第一篇，主要分析ACF与PACF。
1、相关性与相关系数
在统计学中，相关性可以用来衡量两个随机变量之间的线性关系的强度和方向。例如，从经验上可以知道，空气温度和居民用电量之间存在着相关性。
为了能够量化描述这种相关性，人们发明了相关系数。在众多相关系数中，皮尔逊相关系数（Pearson correlation coefficient）是最常用的一个，它是两个变量的协方差与其标准差的乘积之比，取值区间为 [ − 1 , 1 ] [-1,1] [−1,1]。
本文后续中提到的相关系数，均指皮尔逊相关系数。
2、ACF
从字面上的意思理解，自相关函数（Autocorrelation Function）就是用来计算时间序列自身的相关性的函数。
对于一个时间序列 x t x_t xt，根据相关系数的计算公式易得： c o r r ( x t , x t ) = 1 corr(x_t,x_t)=1 corr(xt,xt)=1。这也很好理解，两个完全相同的序列的相关性当然为1。
实际上，在应用自相关函数时，其输入分别为原始的时间序列 x t x_t xt及其 k k k阶滞后序列 x t − k x_{t-k} xt−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}，于是， c o r r ( x t , x t − 1 ) corr(x_t,x_{t-1}) corr(xt,xt−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 k k不断增加，那么就会得到一系列对应的相关系数，将这些相关系数绘制成图，就得到了一个时间序列的自相关系数图。
下图展示了一个价格时间序列，所需数据可以从这里下载。
用如下代码，可以绘制该时间序列的ACF图：
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
df = 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 = α 0 + ∑ i < t α i x i x_{t}=\alpha_0+\sum_{i<t}\alpha_i x_i xt=α0+∑i<tαixi。具体地， x t = α 0 + α 1 x t − 1 x_t=\alpha_0+\alpha_1x_{t-1} xt=α0+α1xt−1表示过去一个时间步的值会影响当前的值， α 0 \alpha_0 α0是截距项， α 1 \alpha_1 α1则表示 x t − 1 x_{t-1} xt−1对当前值影响的程度。
这类似于 y = k x + b y=kx+b y=kx+b的一元线性模型，但需要注意前提条件： x t x_t xt只受其前一步的影响，而更长时间步之前的诸如 x t − 2 x_{t-2} xt−2、 x t − 3 x_{t-3} xt−3等则不会对 x t x_{t} xt产生影响（或者说产生的影响很小以至于可以被忽略）。
如果更长时间步之前的值，例如 x t − 2 x_{t-2} xt−2，也会影响当前值，那么则需要修改模型为： x t = α 0 + α 1 x t − 1 + α 2 x t − 2 x_t=\alpha_0+\alpha_1x_{t-1}+\alpha_2x_{t-2} xt=α0+α1xt−1+α2xt−2。类似地，还可以扩展到时间步为 p p p的情况：
x t = α 0 + ∑ j = 1 p α j x t − j x_t=\alpha_0+\sum_{j=1}^{p}\alpha_jx_{t-j} xt=α0+j=1∑pαjxt−j
这就是 A R ( p ) AR(p) AR(p)模型。
在应用 A R ( p ) AR(p) AR(p)时，需要时刻注意其成立的前提条件：当前时间步的值仅受过去 p p p个时间步的值的影响。也就是说，过去 p p p个时间步的值所包含的信息就（在很大程度上）决定了当前时间步的值。
接下来，一个自然而然的问题就是：如何确定 p p p的值？
3.2 确定阶数的关键
通过上一节的分析，可以得出如下结论：确定 A R AR AR模型中 p p p的值，其实就是在判断至少需要几个时间步的值就可以预测当前值。
为了能够进行判断，接下来还需要一个判断的标准，用于说明为什么选用过去 p p p个时间步的值就可以进行预测而不是 p − 1 p-1 p−1或 p + 1 p+1 p+1。
以 A R ( 2 ) AR(2) AR(2)为例，我们认为 x t − 1 x_{t-1} xt−1和 x t − 2 x_{t-2} xt−2含有决定 x t x_t xt的值的信息，而 x t − 3 x_{t-3} xt−3及更之前的值所含的信息则可以忽略不计。如果用相关性的概念来描述这个例子，则上一句话就变成了“ x t − 1 x_{t-1} xt−1和 x t − 2 x_{t-2} xt−2为与 x t x_t xt的相关性最大的两个，且足够用于预测 x t x_t xt”。
巧了，上文讨论的ACF刚好就是用于计算原始序列和其滞后序列之间相关性的，那相关系数在置信区间之外的前几个序列是否就对应着此处的 p p p呢？
答案是否定的。
举个例子来说明为什么不能直接这样做：假设时间序列确由一个 A R ( 1 ) AR(1) AR(1)模型生成的，也就是说， x t − 1 x_{t-1} xt−1影响 x t x_t xt， x t − 2 x_{t-2} xt−2影响 x t − 1 x_{t-1} xt−1，如果按照ACF的计算方法，那么 x t − 2 x_{t-2} xt−2很可能会通过 x t − 1 x_{t-1} xt−1来“影响” x t x_{t} xt。这就会导致无法判断 x t − 1 x_{t-1} xt−1中是否已经包含了足够的信息来预测 x t x_{t} xt，还是说必须再加上 x t − 2 x_{t-2} xt−2。
也就是说，按照ACF计算得到的 x t x_{t} xt与 x t − 2 x_{t-2} xt−2的相关性并非是“独立的”——因为 x t − 1 x_{t-1} xt−1的存在，导致ACF无法计算 x t − 2 x_{t-2} xt−2到底给预测 x t x_{t} xt带来了多少有用的信息。
因此，为了计算 x t − p x_{t-p} xt−p给 x t x_{t} xt带来了多少“纯粹”的影响，需要发明一种新的计算相关系数的方法。
3.3 PACF
PACF，指偏自相关函数（Partial Autocorrelation Function），是一种计算 x t x_{t} xt和 x t − k x_{t-k} xt−k的相关系数时可以移除 x t − 1 x_{t-1} xt−1到 x t − ( k − 1 ) x_{t-(k-1)} xt−(k−1)影响的方法。
仍以 A R ( 2 ) AR(2) AR(2)为例，来看PACF是如何计算的。用 ϕ ( 2 , 2 ) \phi(2,2) ϕ(2,2)来表示 x t x_{t} xt和 x t − 2 x_{t-2} xt−2的偏自相关系数，则其计算公式如下：
ϕ ( 2 , 2 ) = c o r r ( x t − x ^ t , x t − 2 − x ^ t − 2 ) \phi(2,2)=corr(x_{t}-\hat{x}_t,x_{t-2}-\hat{x}_{t-2}) ϕ(2,2)=corr(xt−x^t,xt−2−x^t−2)
其中， x ^ t \hat{x}_t x^t（ x ^ t − 2 \hat{x}_{t-2} x^t−2）表示通过最小化均方误差建立 x t x_{t} xt（ x t − 2 x_{t-2} xt−2）对 x t − 1 x_{t-1} xt−1的回归方程后所做出的预测值。那么， ϕ ( 2 , 2 ) \phi(2,2) ϕ(2,2)实际上计算的就是两个回归方程的残差序列的相关系数。
我尝试理解了一下这种计算方式的含义，写在这里，仅供参考：
x t − x ^ t x_{t}-\hat{x}_t xt−x^t是用 x t − 1 x_{t-1} xt−1来预测 x t x_{t} xt时的残差，可以表示 x t x_{t} xt中无法被 x t − 1 x_{t-1} xt−1解释的那部分的信息量；类似地， x t − 2 − x ^ t − 2 x_{t-2}-\hat{x}_{t-2} xt−2−x^t−2表示 x t − 2 x_{t-2} xt−2中无法被 x t − 1 x_{t-1} xt−1解释的那部分的信息量。如果它们的相关系数很大，即 x t x_{t} xt中除去 x t − 1 x_{t-1} xt−1影响后的信息与 x t − 2 x_{t-2} xt−2中除去 x t − 1 x_{t-1} xt−1影响后的信息具有很高的相关性，从而也就意味着， x t − 2 x_{t-2} xt−2中除去 x t − 1 x_{t-1} xt−1影响后的信息与 x t x_t xt中除去 x t − 1 x_{t-1} xt−1影响后的信息具有高相关性（ x t − 2 x_{t-2} xt−2对 x t x_t xt的“纯粹”的影响很大），据此可以推断， x t − 2 x_{t-2} xt−2中含有决定 x t x_{t} xt值的信息。
读者可能会有疑问： x t − 2 x_{t-2} xt−2比 x t − 1 x_{t-1} xt−1发生的时间要早，那用后者来预测前者的实际意义是什么呢？这个问题表明，PACF是一种利用数学技术来计算两个序列“纯粹”相关性的方法，并非整个过程都对应着直观的现实理解。
理解了上述关于2阶PACF的计算，不难拓展到 p p p阶的情况，这里直接给出公式：
ϕ 1 , 1 = c o r r ( x t , x t − 1 ) , f o r p = 1 , ϕ k , k = c o r r ( x t − x ^ t , x t − p − x ^ t − p ) , f o r p ≥ 2 \phi_{1,1}=corr(x_t,x_{t-1}),\ for\ p=1, \\ \phi_{k,k}=corr(x_t-\hat x_{t},x_{t-p}-\hat x_{t-p}),\ for\ p\ge 2 ϕ1,1=corr(xt,xt−1), for p=1,ϕk,k=corr(xt−x^t,xt−p−x^t−p), for p≥2
其中， x ^ t \hat x_{t} x^t为 x t x_t xt对 { 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)}做回归后的预测值， x ^ t − p \hat x_{t-p} x^t−p为 x t − p x_{t-p} xt−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 LinearRegression
df['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:, :]) # 残差序列1
lr.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:, :]) # 残差序列2
print(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