1 数据说明
“PriceIndex.xlsx”数据中包含了2018年1月到2021年12月的消费者价格指数（CPI）和工业生产者购进价格指数（IPI），我们采用VAR模型来研究CPI与IPI的相互作用。
2 导入模块
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.api import VAR
3 导入数据
data = pd.read_excel('D:\\Desktop\\PriceIndex.xlsx', sheet_name='Sheet1', na_values='n/a',index_col='日期')
4 单位根检验
data.plot()
plt.rcParams['font.sans-serif'] = ['SimHei'] # 解决中文显示乱码问题！
plt.rc('axes', unicode_minus=False) # 解决坐标轴负号显示乱码问题！
plt.show()
print(sm.tsa.adfuller(data['IPI'], regression='n'))
print(sm.tsa.adfuller(data['IPI'], regression='c'))
print(sm.tsa.adfuller(data['IPI'], regression='ct'))
(0.016766739091273176, 0.6898294359189003, 3, 44, {'1%': -2.6184270867768595, '5%': -1.94847612932006, '10%': -1.6118877404207361}, 82.86442949896774) (-3.5147708911582516, 0.007615751430375056, 1, 46, {'1%': -3.5812576580093696, '5%': -2.9267849124681518, '10%': -2.6015409829867675}, 80.70544149397) (-3.7111082329858274, 0.02163608562780306, 1, 46, {'1%': -4.170389571381606, '5%': -3.51066995808334, '10%': -3.18534353579354}, 79.07825688815427)
print(sm.tsa.adfuller(data['CPI'], regression='n'))
print(sm.tsa.adfuller(data['CPI'], regression='c'))
print(sm.tsa.adfuller(data['CPI'], regression='ct'))
(0.131783283622865, 0.7263638903735252, 4, 43, {'1%': -2.6196969497025417, '5%': -1.9486737067176476, '10%': -1.6117920603217326}, 70.23135564815142) (-5.607803624889886, 1.2228578626811272e-06, 1, 46, {'1%': -3.5812576580093696, '5%': -2.9267849124681518, '10%': -2.6015409829867675}, 60.75688768375825) (-5.533538592002369, 1.9775900163412986e-05, 1, 46, {'1%': -4.170389571381606, '5%': -3.51066995808334, '10%': -3.18534353579354}, 61.60721367469725)
5 构建VAR模型
var_model = VAR(data)
var_results = var_model.fit(maxlags=5, ic='aic') # 根据aic自动选择模型的滞后阶数
print(var_results.summary())
说明：如果想自己指定滞后期(例如滞后2期)，则可以使用var_results = var_model.fit(2)语句！
6 预测CPI和PPI
var_results.plot_forecast(3)
plt.show()
说明：如果想获取预测的数值，可以使用var_results.forecast(data.values[-results.k_ar:], 3)语句！
irf = var_results.irf(10)
irf.plot(orth=False)
plt.show()
fevd = var_results.fevd(10)
fevd.plot()
plt.show()