利用ARMA模型对平稳非白噪声序列进行建模的步骤及代码实现

假如某个观察值序列通过序列预处理可以判定为平稳非白噪声序列&＃xff0c;就可以利用ARMA模型对该序列进行建模。建模的基本步骤如下&＃xff1a;

(1)求出该观察值序列的样本自相关系数(ACF)和样本偏自相关系数(PACF)的值。

(2)根据样本自相关系数和偏自相关系数的性质&＃xff0c;选择适当的ARMA(p&＃xff0c;q)模型进行拟合。

(3)估计模型中位置参数的值。

(4)检验模型的有效性。如果模型不通过检验&＃xff0c;转向步骤(2)&＃xff0c;重新选择模型再拟合。

(5)模型优化。如果拟合模型通过检验&＃xff0c;仍然转向不走(2)&＃xff0c;充分考虑各种情况&＃xff0c;建立多个拟合模型&＃xff0c;从所有通过检验的拟合模型中选择最优模型。

(6)利用拟合模型&＃xff0c;预测序列的将来走势。

二、代码实现

1、绘制时序图&＃xff0c;查看数据的大概分布

trainSeting.head()

Out[36]:

date

2017-10-01 126.4

2017-10-02 82.4

2017-10-03 78.1

2017-10-04 51.1

2017-10-05 90.9

Name: sales, dtype: float64

plt.plot(trainSeting)

2、平稳性检验

&＃39;&＃39;&＃39;进行ADF检验

adf_test的返回值

Test statistic&＃xff1a;代表检验统计量

p-value&＃xff1a;代表p值检验的概率

Lags used&＃xff1a;使用的滞后k&＃xff0c;autolag&＃61;AIC时会自动选择滞后

Number of Observations Used&＃xff1a;样本数量

Critical Value(5%) : 显著性水平为5%的临界值。

(1)假设是存在单位根&＃xff0c;即不平稳&＃xff1b;

(2)显著性水平&＃xff0c;1%&＃xff1a;严格拒绝原假设&＃xff1b;5%&＃xff1a;拒绝原假设&＃xff0c;10%类推。

(3)看P值和显著性水平a的大小&＃xff0c;p值越小&＃xff0c;小于显著性水平的话&＃xff0c;就拒绝原假设&＃xff0c;认为序列是平稳的&＃xff1b;大于的话&＃xff0c;不能拒绝&＃xff0c;认为是不平稳的

(4)看检验统计量和临界值&＃xff0c;检验统计量小于临界值的话&＃xff0c;就拒绝原假设&＃xff0c;认为序列是平稳的&＃xff1b;大于的话&＃xff0c;不能拒绝&＃xff0c;认为是不平稳的

&＃39;&＃39;&＃39;

#滚动统计

def rolling_statistics(timeseries):

#Determing rolling statistics

rolmean &＃61; pd.rolling_mean(timeseries, window&＃61;12)

rolstd &＃61; pd.rolling_std(timeseries, window&＃61;12)

#Plot rolling statistics:

orig &＃61; plt.plot(timeseries, color&＃61;&＃39;blue&＃39;,label&＃61;&＃39;Original&＃39;)

mean &＃61; plt.plot(rolmean, color&＃61;&＃39;red&＃39;, label&＃61;&＃39;Rolling Mean&＃39;)

std &＃61; plt.plot(rolstd, color&＃61;&＃39;black&＃39;, label &＃61; &＃39;Rolling Std&＃39;)

plt.legend(loc&＃61;&＃39;best&＃39;)

plt.title(&＃39;Rolling Mean & Standard Deviation&＃39;)

plt.show(block&＃61;False)

##ADF检验

from statsmodels.tsa.stattools import adfuller

def adf_test(timeseries):

rolling_statistics(timeseries)#绘图

print (&＃39;Results of Augment Dickey-Fuller Test:&＃39;)

dftest &＃61; adfuller(timeseries, autolag&＃61;&＃39;AIC&＃39;)

dfoutput &＃61; pd.Series(dftest[0:4], index&＃61;[&＃39;Test Statistic&＃39;,&＃39;p-value&＃39;,&＃39;#Lags Used&＃39;,&＃39;Number of Observations Used&＃39;])

for key,value in dftest[4].items():

dfoutput[&＃39;Critical Value (%s)&＃39;%key] &＃61; value #增加后面的显著性水平的临界值

print (dfoutput)

adf_test(trainSeting) #从结果中可以看到p值为0.1097>0.1&＃xff0c;不能拒绝H0,认为该序列不是平稳序列

返回结果如下

Results of Augment Dickey-Fuller Test:

Test Statistic    -5.718539e&＃43;00

p-value      7.028398e-07

#Lags Used      0.000000e&＃43;00

Number of Observations Used 6.200000e&＃43;01

Critical Value (1%)  -3.540523e&＃43;00

Critical Value (5%)  -2.909427e&＃43;00

Critical Value (10%)  -2.592314e&＃43;00

dtype: float64

通过上面可以看到&＃xff0c;p值小于0.05&＃xff0c;可以认为该序列为平稳时间序列。

3、白噪声检验

&＃39;&＃39;&＃39;acorr_ljungbox(x, lags&＃61;None, boxpierce&＃61;False)函数检验无自相关

lags为延迟期数&＃xff0c;如果为整数&＃xff0c;则是包含在内的延迟期数&＃xff0c;如果是一个列表或数组&＃xff0c;那么所有时滞都包含在列表中最大的时滞中

boxpierce为True时表示除开返回LB统计量还会返回Box和Pierce的Q统计量

返回值&＃xff1a;

lbvalue:测试的统计量

pvalue:基于卡方分布的p统计量

bpvalue:((optionsal), float or array) – test statistic for Box-Pierce test

bppvalue:((optional), float or array) – p-value based for Box-Pierce test on chi-square distribution

&＃39;&＃39;&＃39;

from statsmodels.stats.diagnostic import acorr_ljungbox

def test_stochastic(ts,lag):

p_value &＃61; acorr_ljungbox(ts, lags&＃61;lag) #lags可自定义

return p_value

test_stochastic(trainSeting,[6,12])

Out[62]: (array([13.28395274, 14.89281684]), array([0.03874194, 0.24735042]))

从上面的分析结果中可以看到&＃xff0c;延迟6阶的p值为0.03<0.05&＃xff0c;因此可以拒绝原假设&＃xff0c;认为该序列不是白噪声序列。

4、确定ARMA的阶数

(1)利用自相关图和偏自相关图

####自相关图ACF和偏相关图PACF

import statsmodels.api as sm

def acf_pacf_plot(ts_log_diff):

sm.graphics.tsa.plot_acf(ts_log_diff,lags&＃61;40) #ARIMA,q

sm.graphics.tsa.plot_pacf(ts_log_diff,lags&＃61;40) #ARIMA,p

acf_pacf_plot(trainSeting) #查看数据的自相关图和偏自相关图

(2)借助AIC、BIC统计量自动确定

##借助AIC、BIC统计量自动确定

from statsmodels.tsa.arima_model import ARMA

def proper_model(data_ts, maxLag):

init_bic &＃61; float("inf")

init_p &＃61; 0

init_q &＃61; 0

init_properModel &＃61; None

for p in np.arange(maxLag):

for q in np.arange(maxLag):

model &＃61; ARMA(data_ts, order&＃61;(p, q))

try:

results_ARMA &＃61; model.fit(disp&＃61;-1, method&＃61;&＃39;css&＃39;)

except:

continue

bic &＃61; results_ARMA.bic

if bic

init_p &＃61; p

init_q &＃61; q

init_properModel &＃61; results_ARMA

init_bic &＃61; bic

return init_bic, init_p, init_q, init_properModel

proper_model(trainSeting,40)

#在statsmodels包里还有更直接的函数&＃xff1a;

import statsmodels.tsa.stattools as st

order &＃61; st.arma_order_select_ic(ts_log_diff2,max_ar&＃61;5,max_ma&＃61;5,ic&＃61;[&＃39;aic&＃39;, &＃39;bic&＃39;, &＃39;hqic&＃39;])

order.bic_min_order

&＃39;&＃39;&＃39;

我们常用的是AIC准则&＃xff0c;AIC鼓励数据拟合的优良性但是尽量避免出现过度拟合(Overfitting)的情况。所以优先考虑的模型应是AIC值最小的那一个模型。

为了控制计算量&＃xff0c;我们限制AR最大阶不超过5&＃xff0c;MA最大阶不超过5。 但是这样带来的坏处是可能为局部最优。

timeseries是待输入的时间序列&＃xff0c;是pandas.Series类型&＃xff0c;max_ar、max_ma是p、q值的最大备选值。

order.bic_min_order返回以BIC准则确定的阶数&＃xff0c;是一个tuple类型

返回值如下&＃xff1a;

order.bic_min_order

Out[13]: (1, 0)

5、建模

从上述结果中可以看到&＃xff0c;可以选择AR(1)模型

################################模型######################################

# AR模型&＃xff0c;q&＃61;0

#RSS是残差平方和

# disp为-1代表不输出收敛过程的信息&＃xff0c;True代表输出

from statsmodels.tsa.arima_model import ARIMA

model &＃61; ARIMA(trainSeting,order&＃61;(1,0,0)) #第二个参数代表使用了二阶差分

results_AR &＃61; model.fit(disp&＃61;-1)

plt.plot(trainSeting)

plt.plot(results_AR.fittedvalues, color&＃61;&＃39;red&＃39;) #红色线代表预测值

plt.title(&＃39;RSS:%.4f&＃39; % sum((results_AR.fittedvalues-trainSeting)**2))#残差平方和

6、预测未来走势

############################预测未来走势##########################################

# forecast方法会自动进行差分还原&＃xff0c;当然仅限于支持的1阶和2阶差分

forecast_n &＃61; 12 #预测未来12个天走势

forecast_AR &＃61; results_AR.forecast(forecast_n)

forecast_AR &＃61; forecast_AR[0]

print (forecast_AR)

print (forecast_ARIMA_log)

[90.49452199 84.05407353 81.92752342 81.22536496 80.99352161 80.91697003

80.89169372 80.88334782 80.88059211 80.87968222 80.87938178 80.87928258]

##将预测的数据和原来的数据绘制在一起&＃xff0c;为了实现这一目的&＃xff0c;我们需要增加数据索引&＃xff0c;使用开源库arrow:

import arrow

def get_date_range(start, limit, level&＃61;&＃39;day&＃39;,format&＃61;&＃39;YYYY-MM-DD&＃39;):

start &＃61; arrow.get(start, format)

result&＃61;(list(map(lambda dt: dt.format(format) , arrow.Arrow.range(level, start,limit&＃61;limit))))

dateparse2 &＃61; lambda dates:pd.datetime.strptime(dates,&＃39;%Y-%m-%d&＃39;)

return map(dateparse2, result)

# 预测从2017-12-03开始&＃xff0c;也就是我们训练数据最后一个数据的后一个日期

new_index &＃61; get_date_range(&＃39;2017-12-03&＃39;, forecast_n)

forecast_ARIMA_log &＃61; pd.Series(forecast_AR, copy&＃61;True, index&＃61;new_index)

print (forecast_ARIMA_log.head())

##绘图如下

plt.plot(trainSeting,label&＃61;&＃39;Original&＃39;,color&＃61;&＃39;blue&＃39;)

plt.plot(forecast_ARIMA_log, label&＃61;&＃39;Forcast&＃39;,color&＃61;&＃39;red&＃39;)

plt.legend(loc&＃61;&＃39;best&＃39;)

plt.title(&＃39;forecast&＃39;)

以上这篇利用python实现平稳时间序列的建模方式就是小编分享给大家的全部内容了&＃xff0c;希望能给大家一个参考&＃xff0c;也希望大家多多支持。