对于时间序列任意时刻的序列值Xt都是一个随机变量，每一个随机变量都会有均值和方差，记Xt的均值为ut，方差为任取t,s，定义序列的自协方差函数=E[(xt-ut)(Xs-us)]和自相关系数，之所以称它们为自协方差函数和自相关系数，是因为它们衡量的是同一个事件在两个不同时期（时刻t和时刻s）之间的相关程度，形象地讲就是度量子集过去的行为对自己现在的影响。
如果事件序列在某一常数附近波动且波动范围有限，即有常数均值和常数方差，并且延迟k期的序列变量的自协方差和自相关系数是相等的或者说延迟k期的序列变量之间的影响程度是一样的，则称为平稳序列。
平稳性的检验：
对序列平稳性的检验有两种检验方法：一种是根据时序图和自相关图的特征做出判断的图检验，该方法操作简单、应用广泛；另一种是构造检验统计量进行的方法，目前最常用的方法是单位根检验。
l 时序图检验：根据平稳时间序列的均值和方差都为常数的性质，平稳序列的时序图显示该序列值始终在一个常数附近随机波动，而且波动的范围有界；如果有明显的趋势性或周期性那它通常不是平稳序列。
l 自相关图检验：平稳序列具有短期相关性，这个性质表明对平稳序列而言通常只有近期的序列值对现时值的影响比较明显，间隔越远的过去值对现时值的影响越小。随着延迟期k的增加，平稳序列的自相关系数会比较快的衰减趋向于零，并在零附近随机波动，而非平稳序列的自相关系数衰减的速度比较慢，这就是利用自相关图进行平稳性检验的标准。
l 单位根检验：单位根检验是指检验序列中是否存在单位根，因为存在单位根就是非平稳性序列了。
from statsmodels.tsa.stattools import adfuller def check_stationary(ts,thresh=0.05): results = adfuller(ts) return results[1]<=thresh check_stationary(train_ts['co2'])
纯随机性检验：
如果一个序列是纯随机序列，那么它的序列值之间应该没有任何关系，即满足 ，这是一种理论上才会出现的理想状态，实际上纯随机序列的样本自相关系数不会绝对为零，但是很接近为零，并在零附近随机波动。
纯随机性检验也称为白噪声检验，一般是构造检验统计量来检验序列的纯随机性，常用的检验统计量有Q统计量、LB统计量，由样本各延迟期数的自相关系数可以计算得到检验统计量，然后计算出对应的p值，如果p值显著大于显著性水平，则表示该序列不能拒绝纯随机性的原假设，可以停止对该序列的分析。
from statsmodels.stats.diagnostic import acorr_ljungbox def check_stochastic(ts): p_value = acorr_ljungbox(ts,lags=1)[1] return p_value check_stochastic(train_ts['co2'])
差分运算：
差分运算具有强大的确定性信息提取能力，许多非平稳序列差分后会显示出平稳序列的性质，这时称这个非平稳序列为差分平稳序列。对差分平稳序列可以使用ARMA模型进行拟合。下图显示了差分前后的情况。
train_ts['diff1'] = train_ts['co2'].diff(1) train_ts.plot(figsize=(15,6))
对差分后的序列同样进行平稳性与随机性验证。
check_stationary(train_ts['diff1'].dropna()) check_stochastic(train_ts['diff1'].dropna())
验证通过后需要进行模型的定阶，确定p、q：
第一种方法是人为识别的方法：将pacf、acf图画出，依据截尾还是拖尾进行查找
import statsmodels.api as sm fig = plt.figure(figsize=(12,8)) ax1 = fig.add_subplot(211) fig = sm.graphics.tsa.plot_acf(train_ts['diff1'].dropna(),lags=20,ax=ax1) ax2 = fig.add_subplot(212) fig = sm.graphics.tsa.plot_pacf(train_ts['diff1'].dropna(),lags=20,ax=ax2)
很多时候，并不能看出p、q的阶数
第二种方法是相对最优模型识别：将p、d、q的可能值进行遍历，一次带入模型求最小的aic、bic、hqic
from statsmodels.tsa.arima_model import ARIMA ts = train_ts['co2'] import itertools p = d =q =range(0,3) pdq = list(itertools.product(p,d,q)) param_list = [] aic_list = [] bic_list = [] results_df = pd.DataFrame() for param in pdq: try: model = ARIMA(ts,order=param) results = model.fit() print 'ARIMA{} - AIC:{} - BIC:{} '.format(param,results.aic,results.bic) param_list.append(param) aic_list.append(results.aic) bic_list.append(results.bic) except: continue results_df['pdq'] = param_list results_df['aic'] = aic_list results_df['bic'] = bic_list results_df.loc[results_df['aic'].idxmin()]
取得p、q、d之后，带入模型
best_result = ARIMA(ts,order=(2,1,2)).fit()
如果预测指定长期序列
pred_ts = best_result.predict(start='2000-01-01',end='2001-12-01',typ='levels')
但这样的话，后期序列值会震荡开。这时可以采用滚动预测：
先创建要预测时间段的时间索引
test_range = pd.date_range(start='2000-01-01', end='2001-12-01', freq='MS')
利用已创建的时间索引，创建序列
pred_ts = pd.Series(index=test_range)
遍历时间索引，对不断迭代更新的序列，进行训练、建模、预测
for cur_month in test_range: prev_month = cur_month - pd.Timedelta('1M') train_ts = proc_data_df.loc[:prev_month]['co2'] model = ARIMA(train_ts, order=(2, 1, 2)).fit() #pred = model.predict(start=cur_month, typ='levels') pred = model.forecast() pred_ts.loc[cur_month] = pred[0]
然后对测试序列、预测序列进行图形上的比较
test_ts.plot(figsize=(15,8))
pred_ts.plot()