ARIMA--基于python实现的代码

数据文件可在github：
http://github.com/aarshayj/Analytics_Vidhya/tree/master/Articles/Time_Series_Analysis 中下载

#1.导入包
import pandas as pd
import numpy as np
import matplotlib.pylab as plt
from matplotlib.pylab import rcParams    #设定画布大小
from statsmodels.tsa.stattools import adfuller      #判断时序数据稳定性的第二种方法
from statsmodels.tsa.seasonal import seasonal_decompose    #分解(decomposing) 可以用来把时序数据中的趋势和周期性数据都分离出来
from statsmodels.tsa.stattools import acf, pacf   #自相关   偏相关
from statsmodels.tsa.arima_model import ARIMA   #ARIMA模型    (p,d,q)
#p--代表预测模型中采用的时序数据本身的滞后数(lags) ,也叫做AR/Auto-Regressive项
#d--代表时序数据需要进行几阶差分化，才是稳定的，也叫Integrated项。
#q--代表预测模型中采用的预测误差的滞后数(lags)，也叫做MA/Moving Average项
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

#2.设定画布的大小
rcParams['figure.figsize'] = 15, 6

#3.读取数据
#data=pd.read_csv('C:\\Users\\Desktop\\AirPassengers.csv')
#print(data.head)
#print ('\n Data types:')
#print(data.dtypes)
#month是object  我们需要将month的类型变为datetime  同时作为index
dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m')
#---其中parse_dates 表明选择数据中的哪个column作为date-time信息，指定日期在哪一列
#---index_col 告诉pandas以哪个column作为 index
#--- date_parser 使用一个function(本文用lambda表达式代替)，使一个string转换为一个datetime变量
data = pd.read_csv('C:\\Users\\Desktop\\AirPassengers.csv', parse_dates=['Month'], index_col='Month',date_parser=dateparse)
#print(data.head)
#print(data.index)

#4.检查时序数据的稳定性
#均值，方差，自方差与时间无关，数据是稳定的
#python判断时序数据稳定性的两种方法
#Rolling statistic-- 即每个时间段内的平均的数据均值和标准差情况。
#Dickey-Fuller Test(迪基-福勒检验) -- 这个比较复杂，大致意思就是在一定置信水平下，对于时序数据假设 Null hypothesis: 非稳定。
#if 通过检验值(statistic)< 临界值(critical value)，则拒绝null hypothesis，即数据是稳定的；反之则是非稳定的。
def test_stationarity(timeseries):
    
    #这里以一年为一个窗口，每一个时间t的值由它前面12个月（包括自己）的均值代替，标准差同理。
    rolmean = pd.rolling_mean(timeseries,window=12)   #均值
    rolstd = pd.rolling_std(timeseries, window=12)    #标准差
    
    #plot rolling statistics:
    fig = plt.figure()#画子图
    fig.add_subplot()  #fig.add_subplot(111)    参数1：子图总行数    参数2：子图总列数   参数3：子图位置
    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 standard deviation')    #标准差
    
    plt.legend(loc = 'best')   #图例
    plt.title('Rolling Mean & Standard Deviation')   #标题
    plt.show(block=False)
    
    
    #Dickey-Fuller test:
    
    print ('Results of Dickey-Fuller Test:')
    dftest = adfuller(timeseries,autolag = 'AIC')
    #dftest的输出前一项依次为检测值，p值，滞后数，使用的观测数，各个置信度下的临界值
    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 = data['#Passengers']
test_stationarity(ts)

#可以看到，数据的rolling均值/标准差具有越来越大的趋势，是不稳定的。
#且DF-test可以明确的指出，在任何置信度下，数据都不是稳定的。

#5.让时序数列变成稳定

#让数据变得不稳定的原因主要有俩：

#趋势（trend）-数据随着时间变化。比如说升高或者降低。
#季节性(seasonality)-数据在特定的时间段内变动。比如说节假日，或者活动导致数据的异常。
#由于原数据值域范围比较大，为了缩小值域，同时保留其他信息，常用的方法是对数化，取log。

ts_log = np.log(ts)
#检测和去除趋势
#通常有三种方法：

#聚合 : 将时间轴缩短，以一段时间内星期/月/年的均值作为数据值。使不同时间段内的值差距缩小。
#平滑： 以一个滑动窗口内的均值代替原来的值，为了使值之间的差距缩小
#多项式过滤：用一个回归模型来拟合现有数据，使得数据更平滑。
#本文主要使用平滑方法

moving_avg = pd.rolling_mean(ts_log,12)
plt.plot(ts_log ,color = 'blue')
plt.plot(moving_avg, color='red')
#可以看出moving_average要比原值平滑许多。

#然后作差：
ts_log_moving_avg_diff = ts_log-moving_avg   #取对数的减去平滑以后的
ts_log_moving_avg_diff.dropna(inplace = True)
test_stationarity(ts_log_moving_avg_diff)    #红色为差值的均值，黑色为差值的标准差

#可以看到，做了处理之后的数据基本上没有了随时间变化的趋势，DFtest的结果告诉我们在95%的置信度下，数据是稳定的。

#上面的方法是将所有的时间平等看待，而在许多情况下，可以认为越近的时刻越重要。
#所以引入指数加权移动平均-- Exponentially-weighted moving average.（pandas中通过ewma()函数提供了此功能。）

# halflife的值决定了衰减因子alpha：  alpha = 1 - exp(log(0.5) / halflife)
expweighted_avg = pd.ewma(ts_log,halflife=12)
ts_log_ewma_diff = ts_log - expweighted_avg
test_stationarity(ts_log_ewma_diff)

#可以看到相比普通的Moving Average，新的数据平均标准差更小了。而且DFtest可以得到结论：数据在99%的置信度上是稳定的。

#检测和去除季节性
#有两种方法：

  #1 差分化： 以特定滞后数目的时刻的值的作差
  #2 分解： 对趋势和季节性分别建模在移除它们

#Differencing--差分
#ts_log_diff = ts_log - ts_log.shift()
#ts_log_diff.dropna(inplace=True)    #删掉有缺失的数据
#test_stationarity(ts_log_diff)

ts_log_diff = ts_log.diff(1)
ts_log_diff.dropna(inplace=True)
test_stationarity(ts_log_diff)
plt.show()
#如图，可以看出相比MA方法，Differencing方法处理后的数据的均值和方差的在时间轴上的振幅明显缩小了。DFtest的结论是在90%的置信度下，数据是稳定的。
#上图看出一阶差分大致已经具有周期性，不妨绘制二阶差分对比：
#diff()函数是求导数和差分的
ts_log_diff1 = ts_log.diff(1)   #一阶差分
ts_log_diff2 = ts_log.diff(2)    #二阶差分
ts_log_diff1.plot()
ts_log_diff2.plot()
plt.show()
#基本已经没有变化。所以使用一阶差分。
#ts_log_diff2 = ts_log.diff(2)
#ts_log_diff2.dropna(inplace=True)
#test_stationarity(ts_log_diff2)
#plt.show()

   #3.Decomposing-分解

#分解(decomposing) 可以用来把时序数据中的趋势和周期性数据都分离出来:

def decompose(timeseries):
    
    # 返回包含三个部分 trend（趋势部分） ， seasonal（季节性部分） 和residual (残留部分)
    decomposition = seasonal_decompose(timeseries)

    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='Seasonarity')
    plt.legend(loc='best')
    plt.subplot(414)
    plt.plot(residual,label='Residual')
    plt.legend(loc='best')
    plt.tight_layout()
    plt.show()
    
    return trend , seasonal, residual

trend , seasonal, residual = decompose(ts_log)

#如图可以明显的看到，将original数据 拆分成了三份。Trend数据具有明显的趋势性，
#Seasonality数据具有明显的周期性，Residuals是剩余的部分，可以认为是去除了趋势和季节性数据之后，稳定的数据，是我们所需要的。

#消除了trend 和seasonal之后，只对residual部分作为想要的时序数据进行处理

residual.dropna(inplace=True)
test_stationarity(residual)

#如图所示，数据的均值和方差趋于常数，几乎无波动(看上去比之前的陡峭，但是要注意他的值域只有[-0.05,0.05]之间)，
#所以直观上可以认为是稳定的数据。另外DFtest的结果显示，Statistic值原小于1%时的Critical value，所以在99%的置信度下，数据是稳定的。

#6. 对时序数据进行预测

#假设经过处理，已经得到了稳定时序数据。接下来，我们使用ARIMA模型

#step1： 通过ACF,PACF进行ARIMA（p，d，q）的p，q参数估计

#由前文Differencing差分部分已知，一阶差分后数据已经稳定，所以d=1。
#所以用一阶差分化的ts_log_diff = ts_log - ts_log.shift() 作为输入。
#等价于
#yt=Yt−Yt−1
#作为输入。

#先画出ACF,PACF的图像,代码如下：
#ACF and PACF plots:
lag_acf = acf(ts_log_diff, nlags=20)   #nlags表示要计算ACF的滞后数
lag_pacf = pacf(ts_log_diff, nlags=20, method='ols')   #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
b44a
.title('Partial Autocorrelation Function')
plt.tight_layout()

#图中，上下两条灰线之间是置信区间，p的值就是ACF第一次穿过上置信区间时的横轴值。q的值就是PACF第一次穿过上置信区间的横轴值。所以从图中可以得到p=2，q=2。

#方法二：信息准则定阶 
#目前选择模型常用如下准则： （其中L为似然函数，k为参数数量，n为观察数） 
#AIC = -2 ln(L) + 2 k 中文名字：赤池信息量 akaike information criterion 
#BIC = -2 ln(L) + ln(n)*k 中文名字：贝叶斯信息量 bayesian information criterion 
#HQ = -2 ln(L) + ln(ln(n))*k hannan-quinn criterion 
#我们常用的是AIC准则，同时需要尽量避免出现过拟合的情况。所以优先考虑的模型应是AIC值最小的那一个模型。
#为了控制计算量，在此限制AR最大阶不超过6，MA最大阶不超过4。 但是这样带来的坏处是可能为局部最优。
#import statsmodels.api as sm
#sm.tsa.arma_order_select_ic(ts,max_ar=6,max_ma=4,ic='aic')['aic_min_order']  # AIC
#sm.tsa.arma_order_select_ic(ts,max_ar=6,max_ma=4,ic='bic')['bic_min_order']  # BIC
#sm.tsa.arma_order_select_ic(ts,max_ar=6,max_ma=4,ic='hqic')['hqic_min_order'] # HQIC

#step2： 得到参数估计值p，d，q之后，生成模型ARIMA（p，d，q）p=2   d=1   q=2
#为了突出差别，用三种参数取值的三个模型作为对比。
#模型1：AR模型(ARIMA(2,1,0))

#from statsmodels.tsa.arima_model import ARIMA
model = ARIMA(ts_log, order=(2, 1, 0))  
results_AR = model.fit(disp=-1)  
plt.plot(ts_log_diff)
plt.plot(results_AR.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_AR.fittedvalues-ts_log_diff)**2))

#图中，蓝线是输入值，红线是模型的拟合值，RSS的累计平方误差。

#模型2：MA模型（ARIMA（0,1,2））

model = ARIMA(ts_log, order=(0, 1, 2))  
results_MA = model.fit(disp=-1)  
plt.plot(ts_log_diff)
plt.plot(results_MA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_MA.fittedvalues-ts_log_diff)**2))

#模型3：ARIMA模型(ARIMA(2,1,2))

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))

#由RSS，可知模型3--ARIMA（2,1,2）的拟合度最好，所以我们确定了最终的预测模型。

#step3: 将模型代入原数据进行预测
#因为上面的模型的拟合值是对原数据进行稳定化之后的输入数据的拟合，所以需要对拟合值进行相应处理的逆操作，使得它回到与原数据一致的尺度。

#ARIMA拟合的其实是一阶差分ts_log_diff，predictions_ARIMA_diff[i]是第i个月与i-1个月的ts_log的差值。
#由于差分化有一阶滞后，所以第一个月的数据是空的，
predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True)
#print (predictions_ARIMA_diff.head)
#累加现有的diff，得到每个值与第一个月的差分（同log底的情况下）。
#即predictions_ARIMA_diff_cumsum[i] 是第i个月与第1个月的ts_log的差值。
predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum()
#先ts_log_diff => ts_log=>ts_log => ts 
#先以ts_log的第一个值作为基数，复制给所有值，然后每个时刻的值累加与第一个月对应的差值(这样就解决了，第一个月diff数据为空的问题了)
#然后得到了predictions_ARIMA_log => predictions_ARIMA
'''

'''
Cumsum ：计算轴向元素累加和，返回由中间结果组成的数组   0是
重点就是返回值是“由中间结果组成的数组”
'''

'''
例题1
arr  = np.array([[[1,2,3],[8,9,12]],[[1,2,4],[2,4,5]]])
print(arr.cumsum(0))
print(arr.cumsum(1))
print(arr.cumsum(2))
例题2
a=np.array([[1,2,3],[4,5,6]])
a.cumsum(0)
a.cumprod(1)
a.cumsum(1)
a.cumprod()
a.cumsum()
 
B = cumsum(A)   这种用法返回数组不同维数的累加和
数组就是一个一行n列的矩阵，向量就是一个n行一列的矩阵
如果A是一个向量， cumsum(A) 返回一个向量，该向量中第m行的元素是A中第1行到第m行的所有元素累加和；
如果A是一个矩阵， cumsum(A) 返回一个和A同行同列的矩阵，矩阵中第m行第n列元素是A中第1行到第m行的所有第n列元素的累加和；
如果A是一个多维数组， cumsum(A)只对A中第一个非奇异维进行计算。
B = cumsum(A,dim)
这种调用格式返回A中由标量dim所指定的维数的累加和。
例如：cumsum(A,1)返回的是沿着第一维（各列）的累加和，cumsum(A,2)返回的是沿着第二维（各行）的累加和。
'''
'''
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 = np.exp(predictions_ARIMA_log)
plt.figure()
plt.plot(ts)
plt.plot(predictions_ARIMA)
plt.title('RMSE: %.4f'% np.sqrt(sum((predictions_ARIMA-ts)**2)/len(ts)))

#ARIMA建模的步骤。
#(1). 获取被观测系统时间序列数据；
#(2). 对数据绘图，观测是否为平稳时间序列；对于非平稳时间序列要先进行d阶差分运算，化为平稳时间序列；
#(3). 经过第二步处理，已经得到平稳时间序列。要对平稳时间序列分别求得其自相关系数ACF 和偏自相关系数PACF，通过对自相关图和偏自相关图的分析，得到最佳的阶层 p 和阶数 q
#(4). 由以上得到的d、q、p，得到ARIMA模型。然后开始对得到的模型进行模型检验。

'''