利用python进行时间序列分析——季节性ARIMA

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

%matplotlib inline
sns.set(color_codes=True)

首先对知识点进行补充:

注:知识点补充参考了《Analysis of Financial Time Series》 RUEY S. TSAY和zhuanlan.zhihu.com/p/21

1.季节性指数

所谓季节指数就是用简单平均法计算的周期内各时期季节性影响的相对数。

x_{ij}=\overline{x}\cdot S_{j}+I_{ij}

首先计算周期内各期平均数 \overline{x}_{k}=\frac{\sum_{i=1}^{n}x_{ik}}{n},k=1,2,..,m

然后计算总平均数 \overline{x}_{k}=\frac{\sum_{i=1}^{n}\sum_{k=1}^{m}x_{ik}}{nm}

再计算季节指数 S_{k}=\frac{\overline {x}_{k}}{\overline{x}},k=1,2,...,m

季节指数反映了该季度与总平均值之间的一种比较稳定的关系:

如果比值大于1,说明该季度的值常常会高于总平均值

如果比值小于1,说明该季度的值常常低于总平均值

如果序列的季节指数都近似为1,就说明该序列没有明显的季节性

2.季节性差分

在有些应用中,季节性是次要的,我们需要把它从数据中消除,这个过程叫季节调整,其中季节性差分化是一种常见的方法。

正规差分化其形式为: \Delta x_{t}=x_{t}-x_{t-1}

我们将它推广,如果一个序列具有周期性,且周期为s,则季节性差分化为:

\Delta_{s} x_{t}=x_{t}-x_{t-s}

3.多重季节性模型

经过了季节性差分和正规差分后,序列成为了平稳时间序列,则我们可以用MA模型对多重差分后的模型建模。则我们有模型:

(1-B^{s})(1-B)x_{t}=(1-\theta B)(1-\Theta B^{s})a_{t}

(1-B^{s})(1-B)x_{t}=a_{t}-\theta a_{t-1}-\Theta a_{t-s}+\theta \Theta a_{t-s-1},\left | \theta \right |<1,\left | \Theta \right |<1

其中,B为向后推移算子,代表前一时刻的值。 B的s次方则为前s时刻的值。 {at}为白噪声序列,s为序列的周期。上述模型为航空模型。

另外,序列经过了季节性差分和正规差分消除序列相关性(序列平稳),则可以认为是间隔s和间隔1的周期共同影响,因此该模型称为多重季节性MA模型

4.参数估计

AR部分是模型简单的正规差分和季节性差分的组合,MA部分有两个参数,于是主要关注MA部分。

w_{t}= (1-\theta B)(1-\Theta B^{s})a_{t}=a_{t}-\theta a_{t-1}-\Theta a_{t-s}+\theta \Theta a_{t-s-1}

我们很容易得到:

E(w_{t})=0 Var(w_{t})=(1+\theta ^{2})(1+\Theta ^{2})\sigma _{a}^{2} Cov(w_{t},w_{t-1})=-\theta(1+\Theta ^{2})\sigma _{a}^{2} Cov(w_{t},w_{t-s+1})=\theta\Theta\sigma _{a}^{2} Cov(w_{t},w_{t-s})=-\Theta(1+\theta^{2})\sigma _{a}^{2} Cov(w_{t},w_{t-s-1})=\theta\Theta\sigma _{a}^{2} Cov(w_{t},w_{t-l})=0, l\neq 0,1,s-1,s,s+1

w_{t} 序列的ACF如下:

\rho_{1} =\frac{-\theta }{1+\theta ^{2}},\rho_{s} =\frac{-\Theta }{1+\Theta ^{2}},\rho_{s-1}=\rho_{s+1}=\rho_{1}\rho_{s}=\frac{\theta\Theta }{(1+\theta ^{2})(1+\Theta ^{2})}

\rho_{l}=0,l>0,l\neq1,s-1,s,s+1

将ACF值以上公式可以求出参数的值。


下面做一个具体的例子:Seasonal ARIMA with Python是对此文的翻译,此外这篇增加了些了理论Statistical forecasting: notes on regression and time series analysis,尤其是在阶数选取上,里面做了个很好的总结。


数据下载

数据为波特兰公共交通系统每月的骑车人数。时间从1973年1月到1982年6月

分析过程



观察数据

数据的整体趋势是什么?有没季节性趋势?这对选择什么模型很重要。假如没有季节性趋势,可以使用ARIMA。假如使用每日数据,有太大波动性性很难发现趋势时,可以使用resample进行每月重采样,或者使用rolling_mean

df = pd.read_csv('portland-oregon-average-monthly-.csv', index_col=0)
df.index.name=None  #将index的name取消
df.reset_index(inplace=True)
df.drop(df.index[114], inplace=True)
start = datetime.datetime.strptime("1973-01-01", "%Y-%m-%d") #把一个时间字符串解析为时间元组
date_list = [start + relativedelta(months=x) for x in range(0,114)]  #从1973-01-01开始逐月增加组成list
df['index'] =date_list
df.set_index(['index'], inplace=True)
df.index.name=None
df.columns= ['riders']
df['riders'] = df.riders.apply(lambda x: int(x)*100)
df.riders.plot(figsize=(12,8), title= 'Monthly Ridership', fontsize=14)
#plt.savefig('month_ridership.png', bbox_inches='tight')


图中发现riders有上升趋势,且有季节性

同时可以使用seasonal_decompose函数进行分析,可以看出季节性非常明显

decomposition = seasonal_decompose(df.riders, freq=12)  
fig = plt.figure()  
fig = decomposition.plot()  
fig.set_size_inches(15, 8)


#可以分别获得趋势、季节性和随机性
trend = decomposition.trend
seasonal = decomposition.seasonal 
residual = decomposition.resid

残差值实际上取消了数据的趋势和季节性,使得这些值与时间无关。可以尝试使用外生变量来对残差进行建模,但是如果尝试将预测的残差值转换回有意义的数值,可能会非常棘手。

对数据进行平稳化

平稳序列是指均值,协方差都不是时间的函数,只与lag有关

为什么需要平稳呢?在线性回归中,我们假设所有观测值是独立的。但在时间序列里,观测是时间依赖的。但假如序列是平稳的,对独立随机变量很好的结果对平稳序列也同样成立,因此我们也可对平稳时间序列进行回归。可以用Dickey-Fuller test检验序列是否平稳

金融领域很多东西之所以难以估计,就是因为其经常突变,根本就不是平稳的


如何深入理解时间序列分析中的平稳性?www.zhihu.com图标


from statsmodels.tsa.stattools import adfuller   #Dickey-Fuller test
def test_stationarity(timeseries):

    #Determing rolling statistics
    rolmean = pd.rolling_mean(timeseries, window=12)
    rolstd = pd.rolling_std(timeseries, window=12)

    #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')  #autolag : {‘AIC’, ‘BIC’, ‘t-stat’, None}
    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)


Results of Dickey-Fuller Test:
Test Statistic                  -1.536597
p-value                          0.515336
#Lags Used                      12.000000
Number of Observations Used    101.000000
Critical Value (1%)             -3.496818
Critical Value (5%)             -2.890611
Critical Value (10%)            -2.582277
dtype: float64

可以看出时间序列不平稳,因此我们要将时间序列转为平稳序列,有如下几种方法:


  1. Deflation by CPI
  2. Logarithmic
  3. First Difference
  4. Seasonal Difference
  5. Seasonal Adjustment


更多信息可以从这里得到


我们首先用1阶差分消除序列的整体趋势

df['first_difference'] = df.riders - df.riders.shift(1)    #也可以使用diff()
test_stationarity(df.first_difference.dropna(inplace=False))


Results of Dickey-Fuller Test:
Test Statistic                  -1.938696
p-value                          0.314082
#Lags Used                      11.000000
Number of Observations Used    101.000000
Critical Value (1%)             -3.496818
Critical Value (5%)             -2.890611
Critical Value (10%)            -2.582277
dtype: float64

看以看出改善了序列平稳,但是还是不平稳,接下来使用季节性差分消除序列的季节性

df['seasonal_difference'] = df.riders - df.riders.shift(12)  
test_stationarity(df.seasonal_difference.dropna(inplace=False))


Results of Dickey-Fuller Test:
Test Statistic                 -2.469741
p-value                         0.123011
#Lags Used                      3.000000
Number of Observations Used    98.000000
Critical Value (1%)            -3.498910
Critical Value (5%)            -2.891516
Critical Value (10%)           -2.582760
dtype: float64

同样看出改善了序列平稳,但是还是不平稳,接下来在1阶差分的基础上使用季节性差分

df['seasonal_first_difference'] = df.first_difference - df.first_difference.shift(12)  
test_stationarity(df.seasonal_first_difference.dropna(inplace=False))


Results of Dickey-Fuller Test:
Test Statistic                -9.258520e+00
p-value                        1.427874e-15
#Lags Used                     0.000000e+00
Number of Observations Used    1.000000e+02
Critical Value (1%)           -3.497501e+00
Critical Value (5%)           -2.890906e+00
Critical Value (10%)          -2.582435e+00
dtype: float64

终于平稳了!!

下面尝试下取对数看看

df.riders_log= df.riders.apply(lambda x: np.log(x))  
test_stationarity(df.riders_log)


Results of Dickey-Fuller Test:
Test Statistic                  -1.677830
p-value                          0.442570
#Lags Used                      12.000000
Number of Observations Used    101.000000
Critical Value (1%)             -3.496818
Critical Value (5%)             -2.890611
Critical Value (10%)            -2.582277
dtype: float64

不平稳,在对数基础上1阶差分?

df['log_first_difference'] = df.riders_log - df.riders_log.shift(1)  
test_stationarity(df.log_first_difference.dropna(inplace=False))


Results of Dickey-Fuller Test:
Test Statistic                  -2.047539
p-value                          0.266126
#Lags Used                      11.000000
Number of Observations Used    101.000000
Critical Value (1%)             -3.496818
Critical Value (5%)             -2.890611
Critical Value (10%)            -2.582277
dtype: float64

还不平稳,在对数基础上做季节性差分!

df['log_seasonal_difference'] = df.riders_log - df.riders_log.shift(12)  
test_stationarity(df.log_seasonal_difference.dropna(inplace=False))


Results of Dickey-Fuller Test:
Test Statistic                  -1.919681
p-value                          0.322860
#Lags Used                       0.000000
Number of Observations Used    101.000000
Critical Value (1%)             -3.496818
Critical Value (5%)             -2.890611
Critical Value (10%)            -2.582277
dtype: float64

也不平稳,继续在对数、1阶差分基础上季节性差分!

df['log_seasonal_first_difference'] = df.log_first_difference - df.log_first_difference.shift(12)  
test_stationarity(df.log_seasonal_first_difference.dropna(inplace=False))


Results of Dickey-Fuller Test:
Test Statistic                -8.882112e+00
p-value                        1.309452e-14
#Lags Used                     0.000000e+00
Number of Observations Used    1.000000e+02
Critical Value (1%)           -3.497501e+00
Critical Value (5%)           -2.890906e+00
Critical Value (10%)          -2.582435e+00
dtype: float64

平稳了~看来做对数对平稳性列没啥用~ ~

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


画ACF\PACF寻找最优的p和q

下图是季节性1阶差分的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)


由于差异序列的ACF在滞后期12(一年后)为负,我应该在模型中使用SMA。尝试不同的项,我发现增加一个SAR项可以提高1982年预测的准确性。通过包含这个项,我可能会过度拟合我的模型。我创建了一个函数,它使用所有可能的参数组合拟合模型,使用这些模型预测多个时间段的结果,然后选择具有最小平方误差总和的模型。


补充:

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


  1. 假如时间序列在很高的lag(10以上)上有正的ACF,则需要很高阶的差分
  2. 假如lag-1的ACF是0或负数或者非常小,则不需要很高的差分,假如lag-1的ACF是-0.5或更小,序列可能overdifferenced。BEWARE OF OVERDIFFERENCING
  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。BEWARE OF USING MULTIPLE AR TERMS AND MULTIPLE MA TERMS IN THE SAME MODEL.
  4. 假如在AR部分有个单位根(AR系数和大约为1),此时应该减少一项AR增加一次差分
  5. 假如在MA部分有个单位根(MA系数和大约为1),此时应该减少一项AR减少一次差分
  6. 假如长期预测出现不稳定,则可能AR、MA系数有单位根

如何确定季节性部分:


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

这部分参考people.duke.edu/~rnau/a

建立模型

综上模型参数选择((0,1,0)x(1,1,1,12), 其实知道参数后,建模是非常容易的

mod = sm.tsa.statespace.SARIMAX(df.riders, trend='n', order=(0,1,0), seasonal_order=(0,1,1,12))
results = mod.fit()
print(results.summary())
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())

做预测

下面简单比较下预测值和真实值

df['forecast'] = results.predict(start = 102, end= 114, dynamic= True)  
df[['riders', 'forecast']].plot(figsize=(12, 8)) 
#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')

参考文献:

Seasonal ARIMA with Python

Statistical forecasting: notes on regression and time series analysis

金融时间序列入门【三】--- 季节模型

《Analysis of Financial Time Series》 RUEY S. TSAY

编辑于 2018-04-06

文章被以下专栏收录