1 基本定义
数据平稳性的图判断:
平稳时间序列的均值和方差都为常数,因此平稳时间序列的时序图应该围绕一条水平线上下波动,而且波动范围有界。
(a)非平稳:有明显的周期性,趋势性
平稳时间序列的序列值之间有短期相关性,则其表现特征是:自相关函数会很快地衰减到零附近
(b)非平稳:自相关函数衰减到零附近的速度比较慢
(c)非平稳:自相关图典型特征,三角对称关系(图1.13)
(d)非平稳 :自相关系数衰减到零的速度非常缓慢,而且呈现明显的周期规律(图1.14)
平稳AR(1)模型例图:
(1)自相关指数按负指数衰减到零附近;
(2)自相关函数呈现正负相见地衰减;
(3)自相关函数具有‘伪周期’的衰减特征;
(4)自相函数衰减速度非常快(短期相关性)。
随机游走序列:
(b)自回归系数
1=1,所以噪声
对xt的影响不随时间的推移而减弱;随机游走在一段时间具有上升或下降趋势,也就是运动方向具有一定的持续性,与(a)相比随机游走序列不围绕某个值上下波动。
(c)方差逐渐增大。
(d)序列值有明显的减小趋势。
(z)故 b,c,d都不是平稳序列。
2 序列预处理
2.1判断序列平稳性
(1)画时间序列图和自相关函数图
# 导入相关模块
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.stats.diagnostic import acorr_ljungbox
%matplotlib inline
plt.rcParams['font.family'] = 'simsun'
plt.rcParams['axes.unicode_minus'] = False
# 绘制时间序列图
Data = pd.read_csv('datas/Guwu.csv',index_col=0).squeeze('columns')
fig = plt.figure(figsize=(12,4),dpi=150)
ax = fig.add_subplot(111)
ax.plot(Data,marker='o',linestyle='-',color='red')
ax.set_xlabel(xlabel='year',fontsize=17)
ax.set_ylabel(ylabel='output',fontsize=17)
plt.xticks(fontsize=15);plt.yticks(fontsize=15)
plt.title('time_1')
plt.tight_layout()
# 绘制自相关图
def ACF(ts, lag):
ts_acf = acf(ts, nlags=lag, fft=False)
fig = plt.figure(figsize=(12,4),dpi=150)
ax = fig.add_subplot(111)
ax.vlines(x=list(range(lag+1)),ymin=np.zeros(lag+1),ymax=ts_acf,color='b',linewidth=2.0)
ax.axhline(y=0, linestyle=':', color='b')
ax.axhline(y=-1.96/np.sqrt(len(ts)),linestyle='--',color='r')
ax.axhline(y=1.96/np.sqrt(len(ts)),linestyle='--',color='r')
ax.set_xlabel(xlabel='lag',fontsize=17)
ax.set_xlabel(xlabel='lag',fontsize=17)
plt.title('ACF_2',fontsize=20)
plt.xticks(fontsize=15);plt.yticks(fontsize=15)
fig.tight_layout()
ACF(guwu_ts,lag=16)
(2) 白噪声检验
# 6,12阶白噪声检验
acorr_ljungbox(Data, return_df=True, lags=[6,12])
判断序列为平稳非白噪声序列,再进行选择模型
2.2选择合适模型
已经画完自相关图,再画偏自相关图(PACF)
# 绘制偏自相关图
def PACF(ts,lag=20,fname='',xlabel=''):
lag_acf = pacf(ts, nlags=lag) # 与自相关函数不同处
plt.vlines(x=[range(lag+1)], ymax=lag_acf, ymin=np.zeros(lag+1), linewidth=2.0, color='b')
plt.axhline(y=0,linestyle=':', color='b')
plt.axhline(y=-1.96/np.sqrt(len(ts)),linestyle='--', color='r')
plt.axhline(y = 1.96/np.sqrt(len(ts)), linestyle='--', color='r')
plt.xlabel(xlabel=xlabel, fontsize=17)
plt.ylabel(ylabel='pacf', fontsize=17)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.title('PACF_3',fontsize=20)
plt.tight_layout()
if not fname:
plt.savefig(fname=fname)
PACF(guwu_ts,lag=8)
根据自相关函数图,偏自相关函数图的拖尾性和截尾性选择合适模型
2.3估计模型参数
# 估计未知参数
guwu_est = ARIMA(guwu_ts,order=(1,0,0)).fit()
print(guwu_est.summary().tables[1])
根据估计值写出模型口径
2.4模型检验
# 残差检验
guwu_resid = guwu_est.resid
acorr_ljungbox(guwu_resid,lags=[6,12],return_df=True, boxpierce=True)
残差检验的p值都大于显著性水平,则拟合模型的残差序列为白噪声序列,然后再进行系数检验
2.5参数检验
from scipy.stats import t
# 过度拟合检验
def pt(x,df):
if x <= 0:
pv = t.cdf(x,int(df))
else:
pv = t.sf(x,int(df))
return pv
# ar.L1系数检验
t1 = 0.3681/0.109;x = pt(t1,15)
# 常数const检验
t0 = 0.8491 /0.055;z = pt(t0,15)
print(x,z)
2.6模型优化
模型优化,从上面自相关图,偏自相关图看出也可选择ARMA(1,2)当作模型,下面对该模型进行模型检验。
guwu_est2 = ARIMA(guwu_ts,order=(1,0,2)).fit()
guwu_resid2 = guwu_est2.resid
acorr_ljungbox(guwu_resid2,lags=[6,12],return_df=True,boxpierce=True)
模型残差检验通过,进一步提取AIC,BIC,AICc,QHIC的值与原来的进行比较,选择值小的
print(guwu_est.aic<guwu_est2.aic, guwu_est.bic<guwu_est2.bic, guwu_est.aicc<guwu_est2.aicc,guwu_est.hqic<guwu_est2.hqic)
2.7预测
from statsmodels.tsa.arima.model import ARIMA
#进行预测最终模型
guwu_pred = guwu_est.get_prediction(start=75,end=79)
confint1 = guwu_pred.conf_int(alpha=0.20) # 80%置信区间
confint2 = guwu_pred.conf_int(alpha=0.05) # 95%置信区间
confint1.columns = ['lower_80', 'upper_80'] # 重命名列
confint2.columns = ['lower_95', 'upper_95'] # 重命名列
confint = pd.concat([confint1, confint2], axis=1)
# 绘制预测图
fig = plt.figure(figsize=(12,4), dpi=150)
ax = fig.add_subplot(111)
ax.plot(guwu_ts[1:],color='y')
guwu_pred.predicted_mean.plot(ax=ax,color='b')
ax.fill_between(confint1.index,confint1.iloc[:,0],confint1.iloc[:,1],color='k',alpha=.1)
ax.fill_between(confint2.index,confint2.iloc[:,0],confint2.iloc[:,1],color='k',alpha=.1)
ax.set_xlabel(xlabel="time", fontsize=17)
plt.legend(loc=2,labels=['real','predicted'],fontsize=13)
plt.xticks(fontsize=15); plt.yticks(fontsize=15)
fig.tight_layout()