python⽩噪声检验-利⽤python实现平稳时间序列的建模⽅式假如某个观察值序列通过序列预处理可以判定为平稳⾮⽩噪声序列,就可以利⽤ARMA模型对该序列进⾏建模。建模的基本步骤如下:
(1)求出该观察值序列的样本⾃相关系数(ACF)和样本偏⾃相关系数(PACF)的值。
miqimiaomiaowu(2)根据样本⾃相关系数和偏⾃相关系数的性质,选择适当的ARMA(p,q)模型进⾏拟合。
(3)估计模型中位置参数的值。
(4)检验模型的有效性。如果模型不通过检验,转向步骤(2),重新选择模型再拟合。
(5)模型优化。如果拟合模型通过检验,仍然转向不⾛(2),充分考虑各种情况,建⽴多个拟合模型,从所有通过检验的拟合模型中选择最优模型。
(6)利⽤拟合模型,预测序列的将来⾛势。
⼆、代码实现
1、绘制时序图,查看数据的⼤概分布
trainSeting.head()
Out[36]:
date
2017-10-01 126.4
2017-10-02 82.4
2017-10-03 78.1
chinglish2017-10-04 51.1
2017-10-05 90.9
Name: sales, dtype: float64
plt.plot(trainSeting)ml是什么意思是
2、平稳性检验
'''进⾏ADF检验
adf_test的返回值
Test statistic:代表检验统计量
p-value:代表p值检验的概率
Lags ud:使⽤的滞后k,autolag=AIC时会⾃动选择滞后rabbit是什么意思
Number of Obrvations Ud:样本数量
Critical Value(5%) : 显著性⽔平为5%的临界值。
(1)假设是存在单位根,即不平稳;
(2)显著性⽔平,1%:严格拒绝原假设;5%:拒绝原假设,10%类推。
(3)看P值和显著性⽔平a的⼤⼩,p值越⼩,⼩于显著性⽔平的话,就拒绝原假设,认为序列是平稳的;⼤于的话,不能拒绝,认为是不平稳的
(4)看检验统计量和临界值,检验统计量⼩于临界值的话,就拒绝原假设,认为序列是平稳的;⼤于的话,不能拒绝,认为是不平稳的
'''
#滚动统计
def rolling_statistics(timeries):
#Determing rolling statistics
rolmean = pd.rolling_mean(timeries, window=12)
rolstd = pd.rolling_std(timeries, window=12)
#Plot rolling statistics:
orig = plt.plot(timeries, 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(block=Fal)
##ADF检验
from statsmodels.tsa.stattools import adfuller
def adf_test(timeries):
rolling_statistics(timeries)#绘图
print ('Results of Augment Dickey-Fuller Test:')
dftest = adfuller(timeries, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Ud','Number of Obrvations Ud'])
for key,value in dftest[4].items():
大学生必看电影
nmec
dfoutput['Critical Value (%s)'%key] = value #增加后⾯的显著性⽔平的临界值
print (dfoutput)
adf_test(trainSeting) #从结果中可以看到p值为0.1097>0.1,不能拒绝H0,认为该序列不是平稳序列
返回结果如下
Results of Augment Dickey-Fuller Test:
lawn mowerTest Statistic? ? -5.718539e+00
p-value? ? ? 7.028398e-07
#Lags Ud? ? ? 0.000000e+00
Number of Obrvations Ud 6.200000e+01
Critical Value (1%)? -3.540523e+00
Critical Value (5%)? -2.909427e+00
Critical Value (10%)? -2.592314e+00
dtype: float64
通过上⾯可以看到,p值⼩于0.05,可以认为该序列为平稳时间序列。
3、⽩噪声检验
'''acorr_ljungbox(x, lags=None, boxpierce=Fal)函数检验⽆⾃相关
lags为延迟期数,如果为整数,则是包含在内的延迟期数,如果是⼀个列表或数组,那么所有时滞都包含在列表中最⼤的时滞中boxpierce为True时表⽰除开返回LB统计量还会返回Box和Pierce的Q统计量
返回值:
lbvalue:测试的统计量
pvalue:基于卡⽅分布的p统计量
bpvalue:((optionsal), float or array) – test statistic for Box-Pierce test
bppvalue:((optional), float or array) – p-value bad for Box-Pierce test on chi-square distribution
'''
from statsmodels.stats.diagnostic import acorr_ljungbox
def test_stochastic(ts,lag):
p_value = acorr_ljungbox(ts, lags=lag) #lags可⾃定义
return p_value
test_stochastic(trainSeting,[6,12])
Out[62]: (array([13.28395274, 14.89281684]), array([0.03874194, 0.24735042]))
从上⾯的分析结果中可以看到,延迟6阶的p值为0.03<0.05,因此可以拒绝原假设,认为该序列不是⽩噪声序列。
4、确定ARMA的阶数
(1)利⽤⾃相关图和偏⾃相关图
####⾃相关图ACF和偏相关图PACF
import statsmodels.api as sm
def acf_pacf_plot(ts_log_diff):
(2)借助AIC、BIC统计量⾃动确定
##借助AIC、BIC统计量⾃动确定
from statsmodels.tsa.arima_model import ARMA
def proper_model(data_ts, maxLag):
init_bic = float("inf")
init_p = 0
init_q = 0
init_properModel = None
for p in np.arange(maxLag):
for q in np.arange(maxLag):
model = ARMA(data_ts, order=(p, q))
try:
results_ARMA = model.fit(disp=-1, method='css')
except:
continue
bic = results_ARMA.bic
if bic < init_bic:
init_p = p
init_q = q童话英文版
init_properModel = results_ARMA
init_bic = bic
return init_bic, init_p, init_q, init_properModel
proper_model(trainSeting,40)
#在statsmodels包⾥还有更直接的函数:
import statsmodels.tsa.stattools as st
order = st.arma_order_lect_ic(ts_log_diff2,max_ar=5,max_ma=5,ic=['aic', 'bic', 'hqic'])
order.bic_min_order
'''
我们常⽤的是AIC准则,AIC⿎励数据拟合的优良性但是尽量避免出现过度拟合(Overfitting)的情况。所以优先考虑的模型应是AIC值最⼩的那⼀个模型。
为了控制计算量,我们限制AR最⼤阶不超过5,MA最⼤阶不超过5。 但是这样带来的坏处是可能为局部最优。
timeries是待输⼊的时间序列,是pandas.Series类型,max_ar、max_ma是p、q值的最⼤备选值。
order.bic_min_order返回以BIC准则确定的阶数,是⼀个tuple类型
返回值如下:
order.bic_min_order
Out[13]: (1, 0)
新东方考研网
5、建模
从上述结果中可以看到,可以选择AR(1)模型pid是什么
>>>>>>##模型>>>>>>>#### AR模型,q=0
#RSS是残差平⽅和
# disp为-1代表不输出收敛过程的信息,True代表输出
from statsmodels.tsa.arima_model import ARIMA
model = ARIMA(trainSeting,order=(1,0,0)) #第⼆个参数代表使⽤了⼆阶差分
results_AR = model.fit(disp=-1)
plt.plot(trainSeting)
plt.plot(results_AR.fittedvalues, color='red') #红⾊线代表预测值
plt.title('RSS:%.4f' % sum((results_AR.fittedvalues-trainSeting)**2))#残差平⽅和