在上一篇季节成分分离中，从图形直观得出基于销量的时间序列为非平稳序列。本文中，将采用
statmodels
包中的函数对样例时间序列进行量化分析。
00 数据预处理
方法同上一篇文章。
读取数据
import pandas as pd
df = pd.read_csv('https://raw.githubusercontent.com/camash/test_data_repository/master/carsales.csv')
season_series = df['日期']
conventional_series = df['传统汽车销量']
nev_series = df['新能源销量（万辆）x5']
total_series= conventional_series+nev_series
01 平稳性判定
使用增强型迪基-富勒检验(Augmented Dickey-Fuller test, 简称ADF)，判断时间序列是否平稳。ADF检验的H0假设是存在单位根，即序列是不平稳的。检验可以使用statmodels.tsa.stattools
包中的adfuller
函数，函数返回值为：
adf：t检验统计量 value: 假设检验结果 usedlag: 使用的滞后阶数 nobs: 用于ADF回归和计算临界值用到的观测值数目 critical values: 1%, 5%, 10%的统计边界值 icbest: 最大的信息准则值
# 首先对总销量进行处理
from statsmodels.tsa.stattools import adfuller as ADF
ADF(total_series)
(-1.022945057961515,
0.7448039049072146,
3,
61,
{'1%': -3.542412746661615,
'5%': -2.910236235808284,
'10%': -2.5927445767266866},
612.9697681954365)
因为结果的pvalue为75%，远大于5%，不能拒绝原假设，所以序列是不平稳的。
02 平稳性处理
平稳性是时间序列分析的前提条件，所以要先将不平稳的序列转换成平稳的序列。常规的方法包括对数变换、平滑、差分和分解等。这里我们使用差分法，方便使用后面的ARIMA模型进行拟合。
# 进行一阶差分
diff_1 = total_series.diff().dropna()
# 进行ADF检测
ADF(diff_1)
(-8.163678389532887,
9.001922543560421e-13,
2,
61,
{'1%': -3.542412746661615,
'5%': -2.910236235808284,
'10%': -2.5927445767266866},
602.6416910108027)
经过一阶差分之后，再使用增强型迪基-富勒检验后得到的pvalue<5%，所以一阶差分后的序列为平稳序列。
03 使用ARIMA模型
ARIMA模型由3个部分组成，AR模型，I模型，MA模型。其中，
AR代表P阶自回归过程，MA代表q阶移动平均过程，I为d阶差分过程。
通过上一小节，只要在一阶差分下，时间序列转化为平稳序列，因此可知。对于p和q的最优取值，可以通过查看自相关图的截尾和偏自相关图的拖尾确定p和q的值。也可以通过程序遍历所有的p和q的组合来确定最优的值。此处采用程序判断。
from statsmodels.tsa.arima_model import ARIMA
import numpy as np
# 制定最大的阶数，不超过观测值的10%
pmax = int(len(total_series)/10)
qmax = pmax
# 使用贝叶斯信息准则表示每个准则的拟合程序
bic_matrix = []
for p in range(pmax+1):
tmp = []
for q in range(qmax+1):
try:
tmp.append(ARIMA(total_series, (p,1,q)).fit().bic)
except:
tmp.append(np.nan)
bic_matrix.append(tmp)
# 转换为dataframe，取最小值的行列，即为最优的p和q值
df_bic_matrix = pd.DataFrame(bic_matrix)
p,q = df_bic_matrix.stack().idxmin()
此处得到使BIC最小的p和q为(0, 1)，代入模型。
model = ARIMA(total_series,(p,1,q)).fit()
# 获取模型报告
model.summary2()
"""
Results: ARIMA
==================================================================
Model: ARIMA BIC: 737.7186
Dependent Variable: D.y Log-Likelihood: -362.62
Date: 2020-03-07 22:25 Scale: 1.0000
No. Observations: 64 Method: css-mle
Df Model: 2 Sample: 1
Df Residuals: 62 5
Converged: 1.0000 S.D. of innovations: 69.434
No. Iterations: 8.0000 HQIC: 733.793
AIC: 731.2419
-------------------------------------------------------------------
Coef. Std.Err. t P>|t| [0.025 0.975]
-------------------------------------------------------------------
const 9.9229 2.2750 4.3616 0.0000 5.4639 14.3819
ma.L1.D.y -0.7564 0.0971 -7.7886 0.0000 -0.9468 -0.5661
---------------------------------------------------------------------------
Real Imaginary Modulus Frequency
---------------------------------------------------------------------------
MA.1 1.3220 0.0000 1.3220 0.0000
==================================================================
"""
使用拟合的模型对未来3季度进行预测。
# 结果的3个数组分别为预测结果、标准误差和置信区间
model.forecast(3)
(array([734.24507481, 744.16794371, 754.09081261]),
array([69.43427006, 71.46411844, 73.43788253]),
array([[598.15640621, 870.33374342],
[604.10084538, 884.23504204],
[610.15520776, 898.02641747]]))
04 使用相同的方式对新能源车进行处理
# 去除前期销售为零的数据
nev_series_clean = nev_series[40:]
ADF(nev_series_clean) # pvalue = 0.6
# 进行一阶差分
nev_diff_1 = nev_series_clean.diff().dropna()
ADF(nev_diff_1) # pvalue = 0.6
# 进行二阶差分
nev_diff_2 = nev_diff_1.diff().dropna()
ADF(nev_diff_2) # pvalue = 0.00016，此时序列为平稳序列
# 已知d为2，求出最优的p,q值
pmax = int(len(nev_series_clean)/10)
qmax = pmax
# 使用贝叶斯信息准则表示每个准则的拟合程序
bic_matrix = []
for p in range(pmax+1):
tmp = []
for q in range(qmax+1):
try:
tmp.append(ARIMA(nev_series_clean, (p,2,q)).fit().bic)
except:
tmp.append(np.nan)
bic_matrix.append(tmp)
# 转换为dataframe，取最小值的行列，即为最优的p和q值
df_bic_matrix = pd.DataFrame(bic_matrix)
p,q = df_bic_matrix.stack().idxmin()
此处得到使BIC最小的p和q为(0, 1)，代入模型。
nev_model = ARIMA(nev_series_clean,(p,2,q)).fit()
# 获取模型报告
nev_model.summary2()
"""
Results: ARIMA
====================================================================
Model: ARIMA BIC: 189.3947
Dependent Variable: D2.新能源销量（万辆）x5 Log-Likelihood: -89.994
Date: 2020-03-07 22:53 Scale: 1.0000
No. Observations: 23 Method: css-mle
Df Model: 2 Sample: 2
Df Residuals: 21 5
Converged: 1.0000 S.D. of innovations: 11.300
No. Iterations: 4.0000 HQIC: 186.845
AIC: 185.9882
--------------------------------------------------------------------
Coef. Std.Err. t P>|t| [0.025 0.975]
--------------------------------------------------------------------
const 0.0444 0.3332 0.1334 0.8952 -0.6086 0.6975
ma.L1.D2.新能源销量（万辆）x5 -1.0000 0.1086 -9.2068 0.0000 -1.2129 -0.7871
-----------------------------------------------------------------------------
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
MA.1 1.0000 0.0000 1.0000 0.0000
====================================================================
"""
使用拟合的模型对未来3季度进行预测。
nev_model.forecast(3)[0]
array([32.79069467, 35.72582405, 38.70538814])
05 预测燃油车销量
将预测的总销量减去新能源车的销量，得到燃油车的销量。
model.forecast(3)[0] - nev_model.forecast(3)[0]
array([701.45438015, 708.44211967, 715.38542447])
06 总结
使用的ARIMA进行时间序列的常规方法如下：
使用ADF判断序列是否平稳，否则进行一阶差分直至平稳，差分次数即差分阶数； 通过程序判断自回归(p)和移动平均阶数(q)的数值； 调用ARIMA模型。
文章转载自FishData，如果涉嫌侵权，请发送邮件至：contact@modb.pro进行举报，并提供相关证据，一经查实，墨天轮将立刻删除相关内容。
