ARIMA时间序列建模过程——原理及python实现.docx

ARIMA时间序列建模过程原理及python实现ARIMA模型的全称叫做自回归查分移动平均模型，全称是(ARIMA, Autoregressive Integrated Moving Average Model)，是统计模型(statistic model)中最常见的一种用来进行时间序列预测的模型，AR、MA、ARMA模型都可以看作它的特殊形式。1. ARIMA的优缺点优点：模型十分简单，只需要内生变量而不需要借助其他外生变量。缺点：要求时序数据是稳定的（stationary），或者是通过差分化(differencing)后是稳定的;本质上只能捕捉线性关系，而不能捕捉非线性关系。2. ARIMA的参数与数学形式ARIMA模型有三个参数:p,d,q。p-代表预测模型中采用的时序数据本身的滞后数(lags) ,也叫做AR/Auto-Regressive项;d-代表时序数据需要进行几阶差分化，才是稳定的，也叫Integrated项;q-代表预测模型中采用的预测误差的滞后数(lags)，也叫做MA/Moving Average项。差分： 假设y表示t时刻的Y的差分。ifd=0,yt=Yt,ifd=1,yt=YtYt1,ifd=2,yt=(YtYt1)(Yt1Yt2)=Yt2Yt1+Yt2ARIMA的预测模型可以表示为：Y的预测值 = 白噪音+1个或多个时刻的加权+一个或多个时刻的预测误差。假设p，q，d已知，ARIMA用数学形式表示为：yt=+1yt1+.+pytp+1et1+.+qetq其中,表示AR的系数，表示MA的系数3.Python建模#构建初始序列import numpy as npimport matplotlib.pyplot as pltimport statsmodels.api as smfrom statsmodels.graphics.tsaplots import acf,pacf,plot_acf,plot_pacffrom statsmodels.tsa.arima_model import ARMAfrom statsmodels.tsa.arima_model import ARIMA#序列化 time_series_ = pd.Series(151.0, 188.46, 199.38, 219.75, 241.55, 262.58, 328.22, 396.26, 442.04, 517.77, 626.52, 717.08, 824.38, 913.38, 1088.39, 1325.83, 1700.92, 2109.38, 2499.77, 2856.47, 3114.02, 3229.29, 3545.39, 3880.53, 4212.82, 4757.45, 5633.24, 6590.19, 7617.47, 9333.4, 11328.92, 12961.1, 15967.61)time_series_.index = pd.Index(sm.tsa.datetools.dates_from_range(1978,2010)time_series_.plot(figsize=(12,8)plt.show()3.1 异常值及缺失值处理异常值一般采用移动中位数方法：from pandas import rolling_median threshold = 3 #指的是判定一个点为异常的阈值 dfpandas = rolling_median(dfu, window=3, center=True).fillna(method=bfill).fillna(method=ffill)#dfu是原始数据，dfpandas 是求移动中位数后的结果，window指的是移动平均的窗口宽度difference = np.abs(dfu - dfpandas)outlier_idx = difference threshold缺失值一般是用均值代替（若连续缺失，且序列不平稳，求查分时可能出现nan）或直接删除。3.2 判断是时序数据是稳定的方法一般是观察时序图（稳定的数据是没有趋势(trend)，没有周期性(seasonality)的; 即它的均值，在时间轴上拥有常量的振幅，并且它的方差，在时间轴上是趋于同一个稳定的值的。该方式并不严谨）或单位跟检测（ADF是一种常用的单位根检验方法，他的原假设为序列具有单位根，即非平稳，对于一个平稳的时序数据，就需要在给定的置信水平上显著，拒绝原假设。如果ADF统计量比临界值的值小，则可在该显著性水平下，拒绝原序列存在单位根的原假设，即原序列是平稳的。 反之，则接受原假设，是非平稳序列）确定。ADF检验：t=sm.tsa.stattools.adfuller(time_series_log, )#ADF检测output=pd.DataFrame(index=Test Statistic Value, p-value, Lags Used, Number of Observations Used,Critical Value(1%),Critical Value(5%),Critical Value(10%),columns=value)outputvalueTest Statistic Value = t0outputvaluep-value = t1outputvalueLags Used = t2outputvalueNumber of Observations Used = t3outputvalueCritical Value(1%) = t41%outputvalueCritical Value(5%) = t45%outputvalueCritical Value(10%) = t410%print(output)#-result-# valueTest Statistic Value 0.807369p-value 0.991754Lags Used 1Number of Observations Used 31Critical Value(1%) -3.66143Critical Value(5%) -2.96053Critical Value(10%) -2.619323.3 不平稳处理1）对数处理。对数处理可以减小数据的波动；time_series_log = np.log(time_series_)time_series_log.plot(figsize=(8,6)v0 = time_series_log0print(np.exp(time_series_log)plt.show()2）差分。一般来说，非纯随机的时间序列经过一阶差分或二阶差分就会变得平稳。在保证ADF检验的p0.01的情况下，阶数越小越好，否则会带来样本减少、还原序列麻烦、预测困难等问题。time_series = time_series_log.diff(1)time_series = time_series.dropna(how=any)time_series.plot(figsize=(8,6)plt.show()3.4 参数确定关于pdq的选择一般是通过计算自相关系数（acf）和偏相关系数确定（pacf），观察两个相关系数图确定，平稳的序列的自相关图和偏相关图要么拖尾，要么是截尾。截尾就是在某阶之后，系数都为0或趋于0，拖尾就是有一个缓慢衰减的趋势，但是不都为0。自相关系数ACF偏相关系数PACF模型选择拖尾p阶截尾AR(p)q阶截尾拖尾MA(q)p阶拖尾q阶拖尾ARMA(p, q)同时，还可利用aic、bic统计量自动确定。plot_acf(time_series)#acf图plot_pacf(time_series)#pacf图plt.show()#bic最小化选择最优参数，aic同理def proper_model(data_ts, maxLag):#init_bic = sys.maxint/python2init_bic = sys.maxsizeinit_p = 0init_q = 0for p in np.arange(maxLag):for q in np.arange(maxLag):model = ARMA(data_ts, order=(p,q)try:results_ARMA = model.fit()except:continuebic = results_ARMA.bicif bic init_bic:init_p = pinit_q = qinit_bic = bicreturn init_bic, init_p, init_q3.5 模型拟合b,p,q = proper_model(time_series, 10)arima_mod = ARIMA(time_series_,order=(p,0,q)arma_mod = ARMA(time_series_,order=(p,q)result_arma = arma_mod.fit()3.6 模型检测假设检验，判断（诊断）残差序列是否为白噪声序列从而确定模型优劣。Ljung-Box test是对randomness的检验,或者说是对时间序列是否存在滞后相关的一种统计检验。对于滞后相关的检验，我们常常采用的方法还包括计算ACF和PCAF并观察其图像，但是无论是ACF还是PACF都仅仅考虑是否存在某一特定滞后阶数的相关。LB检验则是基于一系列滞后阶数，判断序列总体的相关性或者说随机性是否存在。给定置信值a（0.05、0.01.），如果p值大于a值则接受原假定，即序列为白噪声序列，否认序列存在相关性。from statsmodels.stats.diagnostic import acorr_ljungboxp_value = acorr_ljungbox(result_arima.resid.values.squeeze(), lags=5)1print(p_value)3.6 模型预测result_arma_pre = result_arma.predict(start=2011, end=2011)result_arima = arima_mod.fit()result_arima_pre = result_arima.predict(start=2011, end=2011)print(result_arma_pre)print(result_arima_pre)如果对序列使用了差分，最后则需要对结果进行还原操作：# print(result_arma_pre)# v0 = time_series_log0# realV = # realI = # preV = # preI = # realV.append(v0)# preV.append(v0)# for i in range(len(time_series_log):#realI.append(str(time_series_log.indexi).split( )0.split(-)0)#preI.append(str(time_series_log.indexi).split( )0.split(-)0)# for i in range(len(time_series):#v = v0+time_seriesi#v0 = v#realV.append(v)#preV.append(v)# v1 = realVlen(realV)-1# for i in range(len(result_arma_pre):#v = v1 + result_arma_prei#v1 = v#preV.append(v)#