时间序列预测分析(ARIMA)
2016-05-23 本文已影响5171人
26001a36aa12
使用pandas上传下载时间序列
pandas中有专门处理时间序列对象的库
import pandas as pd
import numpy as np
import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6
上传数据集和查看一些最初的行以及列的数据类型
data = pd.read_csv('AirPassengers.csv')
print data.head()
print '\n Data Types:'
print data.dtypes
结果如下：
注意：时间序列对象的读取和数据类型的对象读取不同。因此，为了将读取数据作为时间序列，我们必须通过特殊的参数读取csv指令。
dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m')
data = pd.read_csv('AirPassengers.csv', parse_dates='Month', index_col='Month',date_parser=dateparse)
print data.head()
- 参数含义：
- parse_dates：指定含有时间数据信息的列
- index_col：索引列，并且可以通过data.index检查索引数据类型
- date_parser：将输入字符串转换为可变的时间数据
结果如下：
还可以将列转化为序列对象：
ts = data[‘#Passengers’] ts.head(10)
结果如下：
然后就可以使用时间序列索引技术：
#1. Specific the index as a string constant:
ts['1949-01-01']
||
#2. Import the datetime library and use 'datetime' function:
from datetime import datetime
ts[datetime(1949,1,1)]
还有：
#1. Specify the entire range:
ts['1949-01-01':'1949-05-01']
#2. Use ':' if one of the indices is at ends:
ts[:'1949-05-01']
两种都会输出：
检验时间序列稳定性
稳定性条件：
- 恒定的平均数
- 恒定的方差
- 不随时间变化的自协方差
稳定性测试方法：
- 绘制滚动统计 ：我们可以绘制移动平均数和移动方差，观察它是否随着时间变化。但是，这更多的是一种视觉技术。
- DF检验： 这是一种检查数据稳定性的统计测试。
打印滚动统计量和DF统计量：
在此定义了一个需要时间序列作为输入的函数，用于生成结果。为了保持单元和平均数相似,绘制标准差来代替方差。
from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
#Determing rolling statistics
rolmean = pd.rolling_mean(timeseries, window=12)
rolstd = pd.rolling_std(timeseries, window=12)
#Plot rolling statistics:
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(block=False)
#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
运行结果如下图：
对时间序列进行建模
两种消除趋势和季节性的方法：差分和分解
差分
举例：一阶差分
ts_log_diff = ts_log - ts_log.shift()
plt.plot(ts_log_diff)
结果如下：
验证如下：
ts_log_diff.dropna(inplace=True)
test_stationarity(ts_log_diff)
验证结果：
我们可以看到平均数和标准差随着时间有小的变化。同时，DF检验统计量小于10%的临界值，因此该时间序列在90%的置信区间上是稳定的。同样可以采取二阶或三阶差分在具体情况中获得更好的结果。
分解
分解是一种比较容易解释的手段
消除趋势和季节性：
from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(ts_log)
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
plt.subplot(411)
plt.plot(ts_log, label='Original')
plt.legend(loc='best')
plt.subplot(412)
plt.plot(trend, label='Trend')
plt.legend(loc='best')
plt.subplot(413)
plt.plot(seasonal,label='Seasonality')
plt.legend(loc='best')
plt.subplot(414)
plt.plot(residual, label='Residuals')
plt.legend(loc='best')
plt.tight_layout()
得到如下结果：
再接着检验残差稳定性：
ts_log_decompose = residual
ts_log_decompose.dropna(inplace=True)
test_stationarity(ts_log_decompose)
得到如下结果：
DF测试统计量明显低于1%的临界值，这样时间序列是非常接近稳定。当然也可以尝试高级的分解技术产生更好的结果。
预测时间序列（ARIMA）
1、计算ACF和PACF
代码如下：
#ACF and PACF plots:
from statsmodels.tsa.stattools import acf, pacf
lag_acf = acf(ts_log_diff, nlags=20)
lag_pacf = pacf(ts_log_diff, nlags=20, method='ols')
#Plot ACF:
plt.subplot(121)
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Autocorrelation Function')
#Plot PACF:
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Partial Autocorrelation Function')
plt.tight_layout()
结果如下：
2、确定ARIMA模型
在此例中，应该使用ARIMA(2,1,2)模型
代码如下：
from statsmodels.tsa.arima_model import ARIMA
model = ARIMA(ts_log, order=(2, 1, 2))
results_ARIMA = model.fit(disp=-1)
plt.plot(ts_log_diff)
plt.plot(results_ARIMA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_ARIMA.fittedvalues-ts_log_diff)**2))
预测结果：
3、与原始数据相比较检验
首先，作为一个独立的序列，存储预测结果,观察它。
predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True)
print predictions_ARIMA_diff.head()
结果如下：
注意,这些是从‘1949-02-01’开始,而不是第一个月。这是因为我们将第一个月份取为滞后值，一月前面没有可以减去的元素。
然后，将差分转换为对数尺度的方法是这些差值连续地添加到基本值。一个简单的方法就是首先确定索引的累计总和,然后将其添加到基本值。
predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum()
print predictions_ARIMA_diff_cumsum.head()
结果如下：
接下来,我们将它们添加到基本值。
predictions_ARIMA_log = pd.Series(ts_log.ix[0], index=ts_log.index)
predictions_ARIMA_log = predictions_ARIMA_log.add(predictions_ARIMA_diff_cumsum,fill_value=0)
predictions_ARIMA_log.head()
结果如下：
最后,将指数与原序列比较。
predictions_ARIMA = np.exp(predictions_ARIMA_log)
plt.plot(ts)
plt.plot(predictions_ARIMA)
plt.title('RMSE: %.4f'% np.sqrt(sum((predictions_ARIMA-ts)**2)/len(ts)))
结果如下：
