建模随手记3（2）---时间序列分析

---

接下来主要记录编码解决问题的过程。

在阅读过一些资料后，我认为时间序列分析主要分为五步：

1. 对序列进行稳定性检验，做相关图，观察p值，通过差分等手段获得平稳的时间序列。
2. 选取适合的模型，定阶数。
3. 模型求解。
4. 白噪声检验，优化。
5. 预测。

---

序列稳定性检验

主要参考：

http://xtf615.com/2017/03/08/Python%E5%AE%9E%E7%8E%B0%E6%97%B6%E9%97%B4%E5%BA%8F%E5%88%97%E5%88%86%E6%9E%90/

1. 可以直接肉眼观察法，大致估计是否满足平稳的性质。平均值是否随着时间改变，以及震动的速度是否随着时间改变。

data = pd.DataFrame([[1952, 100], [1953, 101.6], [1954, 103.3], [1955, 111.5], [1956, 116.5],
                     [1957, 120.1], [1958, 120.3], [1959, 100.6], [1960, 83.6], [1961, 84.7],
                     [1962, 88.7], [1963, 98.9], [1964, 111.9], [1965, 122.9], [1966, 131.9],
                     [1967, 134.2], [1968, 131.6], [1969, 132.2], [1970, 139.8], [1971, 142],
                     [1972, 140.5], [1973, 153.1], [1974, 159.2], [1975, 162.3], [1976, 159.1],
                     [1977, 155.1], [1978, 161.2], [1979, 171.5], [1980, 168.4], [1981, 180.4],
                     [1982, 201.6], [1983, 218.7], [1984, 247], [1985, 253.7], [1986, 261.4],
                     [1987, 273.2], [1988, 279.4]], columns=['year', 'gdp'])
plt.figure()
plt.plot(data['year'], data['gdp'])
plt.show()

观察图像可以明显发现，序列均值随着时间变化，因此原序列不满足弱稳定条件，为了时序列稳定我们进行差分处理：

data = data.set_index('year')
data_diff = data.diff(1)
plt.plot(data_diff, label='diff', color='red')

在经过差分之后，我们可以大致观察到，处理之后的序列进本为一个平稳序列。
2.为了定量的确定其平稳性，我们可以画出时序图同时计算p值进行检验。

data = data.set_index('year')
data_diff = data.diff(1)
data_diff.dropna(inplace=True)
plot_acf(data_diff, lags=20)

进行ADF检验：

def adf_test(timerseries):
    print('Result of Augment Dickey-Fuller Test:')
    dftest = adfuller(timerseries, autolag='AIC')
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic', 'p-value', '#Lags Used', 'Number of Observation'])
    for key, value in dftest[4].items():
        dfoutput['Critical Value (%s)' % key] = value
    return dfoutput

print(adf_test(data_diff['gdp']))

结果如下：

Result of Augment Dickey-Fuller Test:
Test Statistic           -3.156056
p-value                   0.022673
#Lags Used                0.000000
Number of Observation    35.000000
Critical Value (1%)      -3.632743
Critical Value (5%)      -2.948510
Critical Value (10%)     -2.613017

参数解释：

Test Static	代表检验统计量，也就是t统计量的值
p-value		带表p值检验的概率，小于0.05则拒绝原假设，那么序列为平稳序列
Lags used	使用的滞后k，计算过程中的延迟阶数
Number of Observation 观测样本数量

通过观察p值发现其通过5%检验，因此拒绝原假设，所以一阶差分后的序列平稳。

选取模型

我们这里使用ARMA模型，通过确定阶数p、q来确定具体模型。

定阶数方法分为BIC和AIC定阶，通过调用python的接口来实现寻找最小的BIC或AIC值。

也可通过直接观察自相关图和偏相关图的截尾和托尾来确定p、q的值。
拖尾：始终有非零取值，不会在k大于某个常数后就很在0附近波动。
截尾：在大于某个常数k后快速趋于0为k阶截尾

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

再结合自相关图和偏向管图分析：

发现自相关图和偏相关图都为托尾，初步确定使用ARMA模型，阶数为（1，1）。

下面计算BIC值来计算p，q的具体值。

import statsmodels.tsa.stattools as st
 order = st.arma_order_select_ic(data_diff, max_ar=3, max_ma=3, ic=['aic', 'bic', 'hqic'])
 print(order.bic_min_order)
 pmax = int(len(data_diff['gdp']) / 10)
 qmax = int(len(data_diff['gdp']) / 10)
 bic_matrix = []
 for p in range(pmax + 1):
     tmp = []
     for q in range(qmax + 1):
         try:
             tmp.append(ARMA(data_diff['gdp'], (p, q)).fit().bic)
         except:
             tmp.append(None)
     bic_matrix.append(tmp)
 bic_matrix = pd.DataFrame(bic_matrix)
 p, q = bic_matrix.stack().idxmin()
 print(p, q)

输出结果为 p=0，q=1，因此构建ARMA（0，1）模型

模型求解

模型求解十分简单，直接调用函数即可。

model = ARMA(data_diff, (0, 1)).fit()

先画出来看看：

plt.figure()
plt.title('gdp&year')
plt.plot(data_diff, label='Origin_diff', color='blue')
model = ARMA(data_diff, (0, 1)).fit()
predict = model.predict(start=0, end=40)
predict.index = 1953+predict.index
plt.plot(predict, label='Predict', color='red')
plt.show()

残差检验

虽然看起来拟合效果还是不错的，但是我们还需要进行残差检验。通过检验模型的残差序列是否为白噪声来判断模型的拟合程度。对残差序列进行ADF检验，同时画出自相关图进行观察。

resid = model.resid
plot_acf(resid)
print(adf_test(resid))
plt.show()

结果：

Result of Augment Dickey-Fuller Test:
Test Statistic          -5.702881e+00
p-value                  7.604291e-07
#Lags Used               0.000000e+00
Number of Observation    3.500000e+01
Critical Value (1%)     -3.632743e+00
Critical Value (5%)     -2.948510e+00
Critical Value (10%)    -2.613017e+00

通过观察ADF检验之后，p通过1%检验，因此可以显著拒绝原假设，所以该序列为稳定的。
通过观察自相关图，可以看到自相关性与k没有明显的趋势和周期。

之后对ARMA（1，1）进行检验之后发现，其拟合效果也很不错，同时通过了白噪声检验。

预测

可用以下两种方法获取预测结果，参数分别为起始值到终点值，第二个预测的步长。

predict = model.predict(start=0, end=45)
forecast = model.forecast(10)

待解决问题

1. 差分序列的恢复没有解决。
2. 还有许多其他模型没有学习。