您可以捐助,支持我们的公益事业。

1元 10元 50元





认证码:  验证码,看不清楚?请点击刷新验证码 必填



  求知 文章 文库 Lib 视频 Code iProcess 课程 认证 咨询 工具 火云堂 讲座吧   成长之路  
会员   
 
   
 
  
每天15篇文章
不仅获得谋生技能
更可以追随信仰
 
     
   
 订阅
  捐助
时间序列预测:使用Python创建季节性ARIMA模型
 
156 次浏览     评价:  
 2018-12-4
 
编辑推荐:
本文来源个人博客,本文将要讨论关于预测的方法。有一种预测是跟时间相关的,而这种处理与时间相关数据的方法叫做时间序列模型。

加载数据,观察数据

%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import datetime
from dateutil.relativedelta import relativedelta
import seaborn as sns
import statsmodels.api as sm
from statsmodels.tsa.stattools import acf
from statsmodels.tsa.stattools import pacf
from statsmodels.tsa.seasonal import seasonal_decompose
df = pd.read_csv('portland-oregon-average-monthly.csv', parse_dates=['month'], index_col='month')
print(df.head())
df['riders'].plot(figsize=(12,8), title= 'Monthly Ridership', fontsize=14)
plt.savefig('month_ridership.png', bbox_inches='tight')

分解数据:

decomposition = seasonal_decompose(df['riders'], freq=12)
fig = plt.figure()
fig = decomposition.plot()
fig.set_size_inches(12, 6)

时间序列稳定化

测试方法:

from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
#Determing rolling statistics
rolmean = timeseries.rolling(window=12).mean()
rolstd = timeseries.rolling(window=12).std()
#Plot rolling statistics:
fig = plt.figure(figsize=(12, 8))
orig = plt.plot(timeseries, color='blue',label='Original')
mean = plt.plot(rolmean, color='red', label='Rolling Mean')
std = plt.plot(rolstd, color='black', label = 'Rolling Std')
plt.legend(loc='best')
plt.title('Rolling Mean & Standard Deviation')
plt.show()
#Perform Dickey-Fuller test:
print('Results of Dickey-Fuller Test:')
dftest = adfuller(timeseries, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
print(dfoutput)

测试序列稳定性:

test_stationarity(df['riders'])

看以看到整体的序列并没有到达稳定性要求,要将时间序列转为平稳序列,有如下几种方法:

Deflation by CPI

Logarithmic(取对数)

First Difference(一阶差分)

Seasonal Difference(季节差分)

Seasonal Adjustment

这里会尝试取对数、一阶查分、季节差分三种方法,先进行一阶差分,去除增长趋势后检测稳定性:

df['first_difference'] = df['riders'].

diff(1)
test_stationarity(df['first_difference'].

dropna(inplace=False))

可以看到图形上看上去变稳定了,但p-value的并没有小于0.05。再来看看12阶查分(即季节查分),看看是否稳定:

df['seasonal_difference'] = df

['riders'].diff(12)
test_stationarity(df['seasonal

_difference'].dropna(inplace=False))

从图形上,比一阶差分更不稳定(虽然季节指标已经出来了),我们再来将一阶查分和季节查分合并起来,再来看下:

df['seasonal_first_difference'] = df

['first_difference'].diff(12)
test_stationarity(df['seasonal_first

_difference'].dropna(inplace=False))

可以看到,在一阶差分后再加上季节差分,数据终于稳定下了。P-value小于了0.05,并且Test Statistic的值远远小于Critical Value (5%)的值。

在《Analysis of Financial Time Series》提到,对原始数据取对数有两个用处:第一个是将指数增长转为线性增长;第二可以平稳序列的方差。再来看看取对数的数据:

df['riders_log']= np.log(df['riders'])
test_stationarity(df['riders_log'])

取对数后,从图形上看有些稳定了,但p-value还是没有接近0.05,没能达稳定性要求。在此基础上我们再来做下一阶差分,去除增长趋势:

df['log_first_difference'] = df

['riders_log'].diff(1)
test_stationarity(df['log_first

_difference'].dropna(inplace=False))

还是没能满足需求,我们再来对对数进行季节性差分:

df['log_seasonal_difference'] = df

['riders_log'].diff(12)
test_stationarity(df['log_seasonal

_difference'].dropna(inplace=False))

还是没能满足需求,再来试试在对数基础上进行一阶差分+季节差分的效果:

df['log_seasonal_first_difference'] =

df['log_first_difference'].diff(12)
test_stationarity(df['log_seasonal_first

_difference'].dropna(inplace=False))

可以看到稳定了满足需求了,但是稳定性相比不取对数的情况下,并没有提升而是下降了。

绘制ACF和PACF图,找出最佳参数

fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(df['seasonal

_first_difference'].iloc[13:], lags=40, ax=ax1)

#从13开始是因为做季节性差分时window是12
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(df['seasonal_

first_difference'].iloc[13:], lags=40, ax=ax2)

fig = sm.graphics.tsa.plot_acf ( df [' seasonal _ first _ difference ' ] . iloc [ 13 : ], lags = 40 , ax = ax1 ) #从13开始是因为做季节性差分时window是12

一直搞不清楚这2张图怎么看。今天一起来学习下:(貌似下面的内容也是糊里糊涂的,后续还需更深入的学习)

如何确定差分阶数和常数项:

1.假如时间序列在很高的lag(10以上)上有正的ACF,则需要很高阶的差分

2.假如lag-1的ACF是0或负数或者非常小,则不需要很高的差分,假如lag-1的ACF是-0.5或更小,序列可能overdifferenced。

3.最优的差分阶数一般在最优阶数下标准差最小(但并不是总数如此)

4.模型不查分意味着原先序列是平稳的;1阶差分意味着原先序列有个固定的平均趋势;二阶差分意味着原先序列有个随时间变化的趋势

5.模型没有差分一般都有常数项;有两阶差分一般没常数项;假如1阶差分模型非零的平均趋势,则有常数项

如何确定AR和MA的阶数:

1.假如PACF显示截尾或者lag-1的ACF是正的(此时序列仍然有点underdifferenced),则需要考虑AR项;PACF的截尾项表明AR的阶数

2.假如ACF显示截尾或者lag-1的ACF是负的(此时序列有点overdifferenced),则需要加MA项,ACF的截尾项表明AR的阶数

3.AR和MA可以相互抵消对方的影响,所以假如用AR-MA模型去拟合数据,仍然需要考虑加些AR或MA项。尤其在原先模型需要超过10次的迭代去converge。

假如在AR部分有个单位根(AR系数和大约为1),此时应该减少一项AR增加一次差分

假如在MA部分有个单位根(MA系数和大约为1),此时应该减少一项AR减少一次差分

假如长期预测出现不稳定,则可能AR、MA系数有单位根

如何确定季节性部分:

1.假如序列有显著是季节性模式,则需要用季节性差分,但是不要使用两次季节性差分或总过超过两次差分(季节性和非季节性)

2.假如差分序列在滞后s处有正的ACF,s是季节的长度,则需要加入SAR项;假如差分序列在滞后s处有负的ACF,则需要加入SMA项,如果使用了季节性差异,则后一种情况很可能发生,如果数据具有稳定和合乎逻辑的季节性模式。如果没有使用季节性差异,前者很可能会发生,只有当季节性模式不稳定时才适用。应该尽量避免在同一模型中使用多于一个或两个季节性参数(SAR + SMA),因为这可能导致过度拟合数据和/或估算中的问题。

创建模型并进行预测

经过上面的分析,可以独处需要制作的模型参数为(0,1,0)x(1,1,1,12):(说实话我也不清楚这个参数是怎么确定的,还需要继续学习, SARIMAX怎么用也还没学习明白)

mod = sm.tsa.statespace.SARIMAX(df['riders'], trend='n', order=(0,1,0), seasonal_order=(1,1,1,12))
results = mod.fit()
print(results.summary())


Statespace Model Results
Dep. Variable: riders No. Observations: 114
Model: SARIMAX(0, 1, 0)x(1, 1, 1, 12) Log Likelihood -501.340
Date: Fri, 16 Nov 2018 AIC 1008.680
Time: 14:13:15 BIC 1016.526
Sample: 01-01-1973 HQIC 1011.856
- 06-01-1982
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.S.L12 0.3236 0.186 1.739 0.082 -0.041 0.688
ma.S.L12 -0.9990 43.079 -0.023 0.981 -85.433 83.435
sigma2 984.7897 4.23e+04 0.023 0.981 -8.19e+04 8.39e+04
Ljung-Box (Q): 36.56 Jarque-Bera (JB): 4.81
Prob(Q): 0.63 Prob(JB): 0.09
Heteroskedasticity (H): 1.48 Skew: 0.38
Prob(H) (two-sided): 0.26 Kurtosis: 3.75

Warnings:
[1] Covariance matrix calculated using the outer

product of gradients (complex-step).

预测值与真实值比较:

df['forecast'] = results.predict(start = 102, end= 114, dynamic= True)
df[['riders', 'forecast']].plot(figsize=(12, 6))
plt.savefig('ts_df_predict.png', bbox_inches='tight')

npredict =df['riders']['1982'].shape[0]
fig, ax = plt.subplots(figsize=(12,6))
npre = 12
ax.set(title='Ridership', xlabel='Date', ylabel='Riders')
ax.plot(df.index[-npredict-npre+1:], df.ix[-npredict-npre+1:, 'riders'], 'o', label='Observed')
ax.plot(df.index[-npredict-npre+1:], df.ix[-npredict-npre+1:, 'forecast'], 'g', label='Dynamic forecast')
legend = ax.legend(loc='lower right')
legend.get_frame().set_facecolor('w')
# plt.savefig('ts_predict_compare.png', bbox_inches='tight')

为了产生针对未来的预测值,我们加入新的时间:

start = datetime.datetime.strptime("1982-07-01", "%Y-%m-%d")
date_list = [start + relativedelta(months=x) for x in range(0,12)]
future = pd.DataFrame(index=date_list, columns= df.columns)
df = pd.concat([df, future])

在新加入的时间上来预测未来值:

df['forecast'] = results.predict(start = 114, end = 125, dynamic= True)
df[['riders', 'forecast']].ix[-24:].plot(figsize=(12, 8))
#plt.savefig('ts_predict_future.png', bbox_inches='tight')

   
156 次浏览  评价: 差  订阅 捐助
相关文章

手机软件测试用例设计实践
手机客户端UI测试分析
iPhone消息推送机制实现与探讨
Android手机开发(一)
相关文档

Android_UI官方设计教程
手机开发平台介绍
android拍照及上传功能
Android讲义智能手机开发
相关课程

Android高级移动应用程序
Android系统开发
Android应用开发
手机软件测试
每天2个文档/视频
扫描微信二维码订阅
订阅技术月刊
获得每月300个技术资源
 
 

关于我们 | 联系我们 | 京ICP备10020922号 京公海网安备110108001071号