[Day7] 用 Python 实作 VAR 多变量时间序列预测

Hey guys, 第七篇就来实作一遍，「以传统统计方法」预测多变量时间序列吧

虽然 VAR 的准确度和复杂度不像我们後面天数要介绍的神经网络一样，但这个实作流程的重点，我觉得在於学习「更了解时间序列的特性和预测的关系」，比方说，平稳性、从时间序列中找自回归系数等。
对基本的东西更熟悉可以避免说，只是会使用神经网络，但不了解时间序列的特性。

像我一开始写这个系列也在挣扎说，要不要直接从机器学习回归开始讲，或是直接介绍各种神经网络、Attention is all you need? 不是吗XD
後来还是觉得应该从基础开始，即便多半是统计和数学，也至少确认自己是理解基本理论的。

好了，不废话了，开始我们的 Python 实作吧

今天的大纲如下：

• 资料集介绍
  • 预测目标
• 时间序列平稳性检测、原因
• 训练/测试集切分
• 自回归模型系数搜寻
• 建立模型
• 预测并判读结果

资料集介绍

我们使用 R 语言的一份范例资料集 Raotbl6，来自於套件库 urca (Unit Root and Cointegration Tests for Time Series Data，一个计量经济学的常用R套件)。这个资料集也是 Yash P. Mehra 曾在其 1994 年论文中使用的资料。

载入资料集

dataset = 'https://raw.githubusercontent.com/selva86/datasets/master/Raotbl6.csv'
df = pd.read_csv(dataset, parse_dates=['date'], index_col='date')
print(df.shape)
df.head()

• rgnp: 实际国民生产总值(real GNP)
• pgnp: 潜在国民生产总值(potential GNP)
• ulc: 单位劳动成本(unit labor cost)
• gdfco: 核心通膨率(core PCE)（剔除季节性食品及能源价格後的个人消费支出物价指数）：一个国家在不同时期内，个人消费支出总水平变动程度的经济指数
• gdf: 国民生产总值(GNP) 缩减指数
• gdfim: 进口缩减指数
• gdfcf: 食品价格通膨率
• gdfce: 能源价格通膨率

预测目标

了解这 8 种经济指标之间的关系，以及纳入指标间的相互影响，预测未来数值

时间序列平稳性检测

何谓平稳性（Stationarity）？

• 资料分布不会随着时间改变而改变，平均数与变异数维持固定，例如白噪音
• 大部分时间序列需要去掉趋势性、做 differencing 等动作，平稳性才会显现

为什麽需要检测时间序列的平稳性？

• 自回归模型中，仅以历史变量和当前变量预测未来；因此需要时间序列变量的基本特性能长时间维持不变，否则以过去预测未来的思路就不成立了。

检测平稳性的方法？

EX: Unit Root Tests

• 单位根检验用於测试时间序列是否非平稳，并具有单位根
• 虚无假设为，假定时间序列存在 unit root（表示是 non-stationary），计算 p-value 验证是否推翻虚无假设

以下我们使用 Augmented Dickey-Fuller Test 进行平稳性检测：

def adf_test(series, title=''):
    print(f'Augmented Dickey-Fuller Test: {title}')
    result = adfuller(series.dropna(), autolag='AIC') # .dropna() handles differenced data
    labels = ['ADF test statistic', 'p-value', 'Number of lags used', 'Number of observations used']
    out = pd.Series(result[0:4], index=labels)
    for k, v in result[4].items():
        out[f'critical value ({k})'] = v

    print(out)
    if result[1] <= 0.05:  # 有显着性，推翻虚无假设
        print("Data has no unit root and is stationary")
    else:
        print("Data has a unit root and is non-stationary")

for col in df.columns:
    adf_test(df[col], title=col)
    print()

几乎所有变量都是不平稳序列，differencing!

# differencing(1)
df_differenced = df.diff().dropna()

# do ADF test on differencing(1) dataframe
for col in df_differenced.columns:
    adf_test(df_differenced[col], title=col)
    print()

还是有部分时间序列变量不平稳，再做一次 differencing
（自动一点的方法可以写成当所有时间序列变量的 p-value 都 ≤ 0.05 时停止 differencing，保存做过差分的次数；这边为了快速示范，就没有特别写～）

# Second Differencing
df_differenced = df_differenced.diff().dropna()

# do ADF test on differencing(2) dataframe
for col in df_differenced.columns:
    adf_test(df_differenced[col], title=col)
    print()

所有时间序列变量都符合平稳性了！

训练/测试集切分

接着我们把要用来 forecasting 的部分切出来保存

steps = 4
train, test = df_differenced[0:-steps], df_differenced[-steps:]
print(train.shape) 
print(test.shape) 

自回归模型系数搜寻

接着搜寻自回归系数 p：

for i in range(1, 11):
    model = VAR(train)
    results = model.fit(i)
    print('Number of Lag  = ', i)
    print('AIC: ', results.aic)
    print('BIC: ', results.bic, '\n')

根据数据，选择 AIC 及 BIC 最小值所在的 Lag = 1

P.S. 这边不一定要用回圈找，也可以使用 VAR model 的 model.select_order(maxlags) 搜寻最佳 order(p)：

model.select_order(maxlags=12).summary()

建立模型

建立 order(p=1) 的自回归模型：

model_fitted = model.fit(1)
model_fitted.summary()

这边呈现了模型的拟合情况和各项系数

预测并判读结果

预测

lag_num = model_fitted.k_ar
test_x = df_differenced.values[-lag_num:]

# forecasting
pred = model_fitted.forecast(y=test_x, steps=steps)
# 加 "_2d" 是为了提醒，预测出的结果是经过 2 次 differencing 的
df_pred = pd.DataFrame(pred, index=df.index[-steps:], columns=df.columns + '_2d')

# 将预测结果还原
for col in df.columns:
    df_pred[f'{col}_1d'] = (df[col].iloc[-steps-1]-df[col].iloc[-steps-2]) + df_pred[f'{col}_2d'].cumsum()
    df_pred[f'{col}_forecast'] = df[col].iloc[-steps-1] + df_pred[f'{col}_1d'].cumsum()
    
test_original = df[-steps:]
test_original.index = pd.to_datetime(test_original.index)

将预测结果和实际画出比较：

fig, axes = plt.subplots(nrows=4, ncols=2, figsize=(15,10))
for i, ax in enumerate(axes.flatten()):
  ax.plot(test_original[df.columns[i]], color='blue', linewidth=1)
  ax.plot(df_forecast[f'{df.columns[i]}_Forecast'], color='red', linewidth=1)
  ax.set_title(df.columns[i])
  ax.spines['top'].set_alpha(0)
  ax.tick_params(labelsize=10)
  ax.legend(['actual', 'forecast'])
fig.tight_layout();

根据图表，我们发现在 rgnp 和 pgnp 的预测效果较好

你也可以引入 R2, RMSE, MSE 等指标来呈现成效好坏～

参考资料

• 资料序列stationary检定方法
• 如何理解时间序列的平稳性？
• Unit root test

好的，总结本篇教学重点：

1. 了解如何使用 python 套件建立自回归模型
2. 了解如何利用资讯指标搜寻最佳自回归系数
3. 了解何谓时间序列平稳性，以及为什麽传统计量经济上需要先做平稳化後，才能预测

感谢你的收看！明天见！