今天给大家分享python实现时间序列的案例。时间序列的模型主要分四种:自回归模型AR,移动回归模型MA,两者的结合移动自回归模型ARMA,以及差分过的差分移动自回归模型ARIMA。
1、AR模型:
Xt时刻的值等于自回归系数乘上对应时刻的数值,ut为时间序列的随机游走。
2、MA模型:
Xt时刻的数值为每个时刻的白噪声的系数的加权和。当自回归和移动回归结合就是ARMA。
3、ARMA模型:
自回归移动平均模型由两部分组成:自回归部分和移动平均部分,因此包含两个阶数,可以表示为ARMA(p,q),p是自回归阶数,q为移动平均阶数,回归方程表示为:
ARMA既体现了每个时刻数据对于预测数据的影响也提现了随机游走因素对于预测数据的影响,模型应用比较常见。
4、ARIMA模型:
差分移动自回归模型,由于前三个模型都有时间序列平稳的假设在,如果时间序列存在明显的上升或者下降趋势,模型预测的效果大大折扣,对于这些有明显下降或者上升趋势的数据集,可以使用差分的方式使它们平稳,之后使用ARMA拟合。假设模型经过d次差分通过了时间序列平稳的检验,ARMA的系数为p,q,ARIMA模型为ARIMA(p,d,q)。
案例:
from __future__ import print_function
import pandas as pd
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import statsmodels.api as sm
#时间序列
dta=[10930,10318,10595,10972,7706,6756,9092,10551,9722,10913,11151,8186,6422,
6337,11649,11652,10310,12043,7937,6476,9662,9570,9981,9331,9449,6773,6304,9355,
10477,10148,10395,11261,8713,7299,10424,10795,11069,11602,11427,9095,7707,10767,
12136,12812,12006,12528,10329,7818,11719,11683,12603,11495,13670,11337,10232,
13261,13230,15535,16837,19598,14823,11622,19391,18177,19994,14723,15694,13248,
9543,12872,13101,15053,12619,13749,10228,9725,14729,12518,14564,15085,14722,
11999,9390,13481,14795,15845,15271,14686,11054,10395]
dta=np.array(dta,dtype=np.float)
#生成时间序列并画图
dta=pd.Series(dta)
dta.plot(figsize=(12,8))
有上升趋势需要差分
#一阶差分
fig = plt.figure(figsize=(12,8))
ax1= fig.add_subplot(111)
diff1 = dta.diff(1)
diff1.plot(ax=ax1)
#二阶差分
fig = plt.figure(figsize=(12,8))
ax2= fig.add_subplot(111)
diff2 = dta.diff(2)
diff2.plot(ax=ax2)
可以观察一阶差分平稳,做单位根检验
#一阶单位根检验
#二阶单位根检验
#选择使用一阶差分的时间序列
diff1= dta.diff(1)
fig = plt.figure(figsize=(12,8))
ax1=fig.add_subplot(211)
ax2 = fig.add_subplot(212)
自相关图和偏相关图
根据上图,猜测有以下模型可以供选择:
1)ARMA(0,1)模型:即自相关图在滞后1阶之后缩小为0,且偏自相关缩小至0,则是一个阶数q=1的移动平均模型;
2.)ARMA(7,0)模型:即偏自相关图在滞后7阶之后缩小为0,且自相关缩小至0,则是一个阶层p=7的自回归模型
3.)ARMA(7,1)模型:即使得自相关和偏自相关都缩小至零。则是一个混合模型。
4) 还可以有其他供选择的模型 ARMA(8,0)
#不同的值的p,q拟合
arma_mod70 = sm.tsa.ARMA(dta,(7,0)).fit()
print(arma_mod70.aic,arma_mod70.bic,arma_mod70.hqic)
arma_mod30 = sm.tsa.ARMA(dta,(0,1)).fit()
print(arma_mod30.aic,arma_mod30.bic,arma_mod30.hqic)
arma_mod71 = sm.tsa.ARMA(dta,(7,1)).fit()
print(arma_mod71.aic,arma_mod71.bic,arma_mod71.hqic)
arma_mod80 = sm.tsa.ARMA(dta,(8,0)).fit()
print(arma_mod80.aic,arma_mod80.bic,arma_mod80.hqic)
看不同的信息量:
选择ARMA(8,0)
resid = arma_mod80.resid
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
plt.show()
残差的ACF和PACF图,可以看到序列残差基本为白噪声。
德宾-沃森检验:
print(sm.stats.durbin_watson(arma_mod80.resid.values))
QQplot:
print(stats.normaltest(resid))
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
fig = qqplot(resid, line='q', ax=ax, fit=True)
#plt.show()
预测:
predict_dta = arma_mod80.predict('2090', '2150', dynamic=True)
print(predict_dta)
fig, ax = plt.subplots(figsize=(12, 8))
ax = dta.ix['2000':].plot(ax=ax)
fig = arma_mod80.plot_predict('2090', '2150', dynamic=True, ax=ax, plot_insample=False)
plt.show()
时间序列的初步建模并预测了之后的数据。
领取专属 10元无门槛券
私享最新 技术干货