使用Python构造VARIMA模型

简介

VARMA(p,q)结合了VAR和VMA模型，其中p是向量自回归(VAR)模型的滞后期数，q是VMA模型的移动平均的阶数。

VARMA是ARMA的推广，它将ARMA模型扩展到多个时间序列变量的情况，通过VAR和VMA的线性组合来描述多个时间序列变量之间的联合变化，适合描述多个时间序列变量之间的关系。 时间序列变量。

通过将 q 参数设置为 0，VARMA 模型可以像 VAR 模型一样工作；通过将 p 参数设置为 0，它也可以像 VMA 模型一样工作。VARMA 也不能处理非平稳金融时间序列数据。 矢量自回归积分移动平均（VARIMA）是一种经历差分过程的VARMA模型。

首先，应用约翰森检验（Johansen test ），结果表明英国的GDP、失业率和CPIH之间存在长期均衡关系。
因此，它们可以作为协变量来预测GDP。 正如ARIMA模型分析中提到的，GDP时间序列是不稳定的，因此必须在时间序列中实施一阶差分。 然后，还应用归一化过程。使用 MinMaxScaler() 函数后，数据将缩放到特定范围。 然后，将归一化后的数据按比例划分为训练集和测试集。
通过使用VARMAX功能，它将能够自动与AIC标准进行比较并找到最佳模型。 最优模型将具有最小的 AIC 值。 此外，预测的GDP值需要进行非标准化处理，以便与原始数据进行比较。

代码构建

首先导入需要用到的Python包：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.statespace.varmax import VARMAX
from sklearn.preprocessing import MinMaxScaler
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error

然后读取.csv文件的时序数据，这里使用了英国的GDP数据，CPI(通货膨胀率)和Unemployment rate(失业率)作为covariate(协变量)。

# 1. 读取csv时序数据
gdp_data = pd.read_csv('datasets/UK_GDP.csv')[["GDP"]]
inflation_data = pd.read_csv('datasets/UK_inflation.csv')[["Inflation"]]
unemployment_data = pd.read_csv('datasets/UK_unemployment.csv')[["Unemployment"]]
data_origin = gdp_data.copy()

接着对所有数据进行一阶差分，使其稳定（因为之前的博客已经对同数据进行过检测，并确定数据不稳定，所以要进行差分）。

gdp_data = gdp_data.diff().dropna()
inflation_data = inflation_data.diff().dropna()
unemployment_data = unemployment_data.diff().dropna()

然后对处理过的数据进行归一化。

scaler1 = MinMaxScaler()
scaled_gdp_data = pd.DataFrame(scaler1.fit_transform(gdp_data), columns=gdp_data.columns, index=gdp_data.index)
scaler2 = MinMaxScaler()
scaled_inflation_data = pd.DataFrame(scaler2.fit_transform(inflation_data), columns=inflation_data.columns, index=inflation_data.index)
scaler3 = MinMaxScaler()
scaled_unemployment_data = pd.DataFrame(scaler3.fit_transform(unemployment_data), columns=unemployment_data.columns, index=unemployment_data.index)

然后自动定阶，通过AIC找出最合适的参数。

merged_data = pd.concat([scaled_gdp_data, scaled_inflation_data, scaled_unemployment_data], axis=1)
train_size = int(len(merged_data))-3
train_data, test_data = merged_data[:train_size], merged_data[train_size:]
 
best_aic = np.inf
best_order = None
best_model = None
pq_range = range(2) # 取值范围
for p in pq_range:
    for q in pq_range:
        try:
            model = VARMAX(train_data, order=(p, q))
            result = model.fit()
            aic = result.aic
            if aic < best_aic:
                best_aic = aic
                best_order = (p, q)
                best_model = result
        except:
            continue
 
print("Best order:", best_order)
print("Best AIC:", best_aic)

使用VARIMA模型进行预测，打印预测值和真实值的对比图，并计算模型RMSE和MAPE指标。

gdp_predictions = best_model.forecast(steps=len(test_data))[['GDP']]
gdp_predictions = pd.DataFrame(gdp_predictions, columns=['GDP'], index=test_data.index-1)
gdp_predictions = scaler1.inverse_transform(gdp_predictions)
actual = scaler1.inverse_transform(test_data[['GDP']])
 
actual = np.array(gdp_data[-3:].cumsum() + data_origin.values[127])
predictions = gdp_predictions.cumsum() + data_origin.values[127]
 
plt.figure()
plt.plot(actual, label='Actual')
plt.plot(predictions, label='Predicted')
plt.legend()
plt.show()
 
rmse = np.sqrt(mean_squared_error(actual, predictions))
mape = mean_absolute_percentage_error(actual, predictions)
print(f"RMSE: {rmse}")
print(f"MAPE: {mape}")