基于python的非平稳时间序列模型

前言

平稳时间序列指的是宽平稳时间序列，就是指时间序列的均值、方差和协方差等一二阶矩存在但不随时间改变，表现为时间的常数。若三个条件有一个不成立，那么就称该序列为非平稳时间序列。包括确定性趋势时间序列和随机性趋势时间序列。
要想把非平稳的时间序列转化为平稳的时间序列，需要去趋势和差分方法消除确定性趋势和随机性趋势。
实际数据分析中，一阶差分提取线性趋势、二阶或三阶等地阶差分提取曲线趋势，对于含有季节趋势的数据，通常选取差分的步长等于季节的周期可以较好地提取季节信息。
处理确定性趋势时间序列最好采用减去趋势部分的方法来去趋势，特别是对于曲线趋势明显的方差齐性序列来讲，采用此方法更好；而处理随机性趋势的时间序列最好使用差分运算来去趋势。
一般来讲，具有随机性趋势的非平稳时间序列在经过适当差分之后会变成一个平稳时间序列。此时，我们称这个非平稳序列为差分平稳序列。对差分平稳序列可以使用求和自回归移动平均模型(ARIMA)进行拟合。对非平稳时间序列的思路是转化为平稳序列，如果检验是非平稳的序列，对其进行差分运算，直到检验是平稳的；如果检验是平稳的，那么转入ARMA的建模步骤。
对序列作一阶差分可以消除线性趋势；作二阶、三阶等低阶差分可以消除曲线趋势，但是我们也知道这样做会有过差分风险，即人为造成的非平稳和不好的信息。在进行ARIMA建模时，一方面要看到差分运算能够充分提取确定性信息，另一方面也要看到差分运算的解释性不强，同时有过差分风险。
如果担心浪费残差信息。当残差中含有自相关系数时，可继续对残差序列建立自回归模型，这样即自然的提出了残差自回归模型。.
在进行残差自回归模型建模时，首先拟合趋势项和季节变化项，然后进行残差检验。当残差序列自相关不显著时，则建模结束；当残差序列自相关性显著时，再对残差进行建模。残差的建模步骤与ARIMA模型建模步骤一致。

代码

#一阶差分运算
elec_prod = np.loadtxt("elec_prod.txt")
index = pd.date_range(start="1974", end="2007", freq="M")
elec_df = pd.Series(elec_prod,index=index).diff()

fig = plt.figure(figsize=(12,4), dpi=150)
ax = fig.add_subplot(111)
ax.plot(elec_df, linestyle="-", color="blue")
ax.set_ylabel(ylabel="一阶差分序列值",fontsize=17)
ax.set_xlabel(xlabel="时 间",fontsize=17)
plt.xticks(fontsize=15); plt.yticks(fontsize=15)
fig.tight_layout(); plt.savefig(fname='fig/5_3.png')

#二阶差分运算
ningxia_gdp = np.loadtxt("ningxia_gdp.txt")
index = pd.date_range(start="2001", end="2021", freq="Y")
ningxia_gdp_df = pd.Series(ningxia_gdp,index=index).diff().diff()

fig = plt.figure(figsize=(12,4),dpi=150)
ax = fig.add_subplot(111)
ax.plot(ningxia_gdp_df, marker="o", linestyle="-", color="blue")
ax.set_ylabel(ylabel="二阶差分序列值", fontsize=17)
ax.set_xlabel(xlabel="时间", fontsize=17)
plt.xticks(fontsize=15); plt.yticks(fontsize=15)
fig.tight_layout(); plt.savefig(fname='fig/5_4.png')

#首先做一阶差分取消线性趋势，然后做四步差分提取季节趋势
Index = pd.date_range(start="2013", end="2017-06-30", freq="Q")
gdp_df = pd.read_csv('JDGDP.csv'); gdp_df.index = Index

cgdp = gdp_df.diff().diff(periods=4)
fig = plt.figure(figsize=(12,4),dpi=150)
ax = fig.add_subplot(111)
ax.plot(cgdp, marker="o", linestyle="-", color="blue")
ax.set_ylabel(ylabel="差分序列值", fontsize=17)
ax.set_xlabel(xlabel="时间", fontsize=17)
plt.xticks(fontsize=15); plt.yticks(fontsize=15)
fig.tight_layout(); plt.savefig(fname='fig/5_5.png')

#序列作一阶差分运算，对所得差分序列作出时序图、自相关图和偏自相关图
tr_data = np.loadtxt("tr_industry.txt")
Index = pd.date_range(start="1995", end="2015", freq="Y")
tr_ts = pd.Series(tr_data,index=Index)
tr_diff = tr_ts.diff()

fig = plt.figure(figsize=(12,6), dpi=150)
ax1 = fig.add_subplot(221)
ax1.plot(tr_ts,marker='o', linestyle='-', color='b')
ax1.set_ylabel(ylabel="第三产业增加值", fontsize=17)
ax1.set_xlabel(xlabel="图 5.7 第三产业增加值序列的时序图")
plt.xticks(fontsize=15); plt.yticks(fontsize=15)
ax2 = fig.add_subplot(222)
ax2.plot(tr_diff, marker='o', linestyle='-', color='b')
ax2.set_ylabel(ylabel="一阶差分序列", fontsize=17)
ax2.set_xlabel(xlabel="图 5.8 差分序列的时序图")
plt.xticks(fontsize=15); plt.yticks(fontsize=15)
ax3 = fig.add_subplot(223)
ACF(tr_diff[1:], lag=8)
ax4 = fig.add_subplot(224)
PACF(tr_diff[1:], lag=8, xlabel='lag',fname="fig/5_7.png")

import statsmodels.api as sm
tr_res = sm.tsa.SARIMAX(tr_ts, order=(0, 1, 1))
tr_est = tr_res.fit()
print(tr_est.summary())

acorr_ljungbox(tr_est.resid[1:],lags = [2,4,6,8],boxpierce=True, return_df=True)

tr_fore = tr_est.get_forecast()
confint = pd.concat([tr_fore.summary_frame(alpha=0.20),tr_fore.summary_frame().iloc[:,2:]],axis=1,ignore_index=False)
print(confint)

fig = plt.figure(figsize=(12,4),dpi=150)
ax = fig.add_subplot(111)
ax.plot(tr_ts, marker="o", linestyle="-", color="blue")
fcast1 = tr_est.get_forecast(2).summary_frame()
fcast1['mean'].plot(ax=ax, marker="o", color="red")
fcast2 = tr_est.get_forecast(steps=2).summary_frame(alpha=0.2)
ax.fill_between(fcast1.index, fcast1['mean_ci_lower'], fcast1['mean_ci_upper'], color='green', alpha=0.3)
ax.fill_between(fcast2.index, fcast2['mean_ci_lower'], fcast2['mean_ci_upper'], color='black', alpha=0.5)
ax.legend(["Real Values","Forecast"],loc="upper left",fontsize=13)
ax.set_ylabel(ylabel="第三产业增加值", fontsize=17)
ax.set_xlabel(xlabel="时 间", fontsize=17)
plt.xticks(rotation=360,fontsize=15); plt.yticks(fontsize=15)
fig.tight_layout(); plt.savefig(fname="fig/5_11.png")