python建立ARIMA模型进行时间序列分析（氵论文）

如何从零开始氵一篇数据处理的论文

我是真的不怎么会，纯粹入门级，画个图交个论文就完事，下面介绍我到底时怎么做的（仅参考）

并且，一直给我报这个错，不过感觉好像问题不大，估计是因为时间间隔不相等引起的

ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.

CSV数据

注意一下日期格式，用excel打开的时候会变成默认日期格式，记得把斜杠换成 - ，有时候会怪怪的

坝顶垂直位移监测是变形监测的一项主要内容，其成果能准确反映水工建筑物垂直向的变形性态。下表是我国某重力坝坝顶5#104点2008年1月~2012年1月坝顶垂直位移，原始数据共50组。

Date,dy
2008-01-11,0.62
2008-02-06,1.01
2008-03-17,1.78
2008-04-13,1.29
2008-05-17,0.11
2008-06-15,-0.35
2008-07-01,-0.44
2008-07-12,-0.3
2008-08-10,-1.11
2008-09-14,-1.78
2008-10-12,-1.39
2008-11-16,-0.94
2008-12-13,-0.36
2009-01-19,1.47
2009-02-16,1.75
2009-03-13,2.04
2009-04-18,1.03
2009-05-16,0.02
2009-06-20,-0.59
2009-07-11,-1.35
2009-08-15,-2.14
2009-09-19,-1.96
2009-10-16,-1.46
2009-11-14,-0.56
2009-12-11,0.04
2010-01-15,0.96
2010-02-20,1.58
2010-03-13,1.43
2010-04-17,0.95
2010-05-15,0.14
2010-06-12,-0.3
2010-07-16,-1.35
2010-08-14,-1.6
2010-09-18,-1.98
2010-10-16,-1.58
2010-11-19,-0.98
2010-12-24,0.56
2011-01-21,1.14
2011-02-18,1.19
2011-03-19,1.18
2011-04-17,0.61
2011-05-15,0.76
2011-06-18,-0.66
2011-07-16,-1.14
2011-08-20,-1.35
2011-09-24,-1.85
2011-10-22,-0.95
2011-11-19,-0.65
2011-12-24,0.44
2012-01-14,1.09

实验过程是用50组数据来拟合模型，选出一个最好的；然后再用前四十组来建立模型，预测下来十组数据，并且画个图来看。其实我自己也觉得不怎么严谨。别问 问就懒得改了

程序

因为我的前后矛盾 所以我写了两程序，一个专门用来预测；我打注释的也是可以跑的，调程序嫌麻烦就直接注释掉了。其实这程序是用来画图的（迫真）

拟合

import warnings  # 忽略警告
warnings.filterwarnings('ignore')


import pandas as pd #表格和数据操作
import numpy as np
import matplotlib.pyplot as plt
from random import randrange
from statsmodels.tsa.stattools import adfuller as ADF
from statsmodels.stats.diagnostic import acorr_ljungbox #白噪声检验
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima_model import ARMA
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.api import tsa
import statsmodels.api as sm
from statsmodels.graphics.api import qqplot
from itertools import product

dta=pd.read_csv('时间序列1 - 副本 - 副本.CSV',header=0 ,index_col=0)



dta = dta.dropna() #去除不完整项数据

print(dta.describe())

#用单位根检验（ADF）进行平稳性检验

dta1= dta
diff=0
adf=ADF(dta1)
if adf[1]>0.05:
    print (u'原始序列经检验不平稳，p值为:%s'%(adf[1]))
else:
    print (u'原始序列经检验平稳，p值为:%s'%(adf[1]))
while adf[1]>=0.05:#adf[1]为p值，p值小于0.05认为是平稳的
    diff=diff+1
    adf=ADF(dta1.diff(diff).dropna())
print (u'原始序列经过%s阶差分后归于平稳，p值为%s'%(diff,adf[1]))

for i in range(1,3):
    dta1 = dta1.diff(i).dropna()
    adf = ADF(dta1)
    print(u'原始序列经过%s阶差分,p值为%s'%(i,adf[1]))

#采用LB统计量的方法进行白噪声检验

dta2 = dta
[[lb],[p]]=acorr_ljungbox(dta2,lags=1)
if p<0.05:
    print (u'原始序列为非白噪声序列，p值为:%s'%p)
else:
    print (u'原始序列为白噪声序列，p值为:%s'%p)
[[lb],[p]]=acorr_ljungbox(dta2.diff(1).dropna(),lags=1)
if p<0.05:
    print (u'一阶差分序列为非白噪声序列，p值为:%s'%p)
else:
    print (u'一阶差分序列为白噪声序列，p值为:%s'%p)


dta.plot(figsize=(8,7))
plt.show()


#绘制自相关和偏向关图
fig = plt.figure(figsize=(8,7))
ax1= fig.add_subplot(211)
ax2= fig.add_subplot(212)
fig = plot_acf(dta,ax=ax1)
fig = plot_pacf(dta,ax=ax2)
fig.show()


#模型识别
#确定最佳p,d,q值

###定阶
##pmax = int(len(dta)/10) #一般阶数不超过length/10
##qmax = int(len(dta)/10) #一般阶数不超过length/10
##bic_matrix = [] #bic矩阵
##for p in range(pmax+1):
##    tmp = []
##    for q in range(qmax+1):
##        try: #存在部分报错，所以用try来跳过报错。
##            tmp.append(ARIMA(dta, (p,1,q)).fit().bic)
##        except:
##            tmp.append(None)
##    bic_matrix.append(tmp)
##
##bic_matrix = pd.DataFrame(bic_matrix) #从中可以找出最小值
##
##p,q = bic_matrix.stack().astype('float64').idxmin() #先用stack展平，然后用idxmin找出最小值位置。
##print (u'BIC最小的p值和q值为：%s、%s' %(p,q))

p=2  #这里是已经出结果了，212模型最好，所以就直接赋值了，免得调程序调半天
q=2

#建立ARIMA(p,1,q)模型
arima = ARIMA(dta.dropna(), (p,1,q)).fit()
print("最优模型", arima.summary())

#模型检验：模型确立后，检验其残差序列是否为白噪声

dta3 = dta.drop(axis = 0, index = '2008/1/11')#删除首项，对应差分缺失
dta_pred = arima.predict(typ = 'levels') #按模型预测
print(dta_pred)#手动操作。。
#绘拟合图
fig3 = plt.figure(figsize=(12,8))
plt.plot(dta3)
plt.plot(dta_pred, color='red')
fig3.show()

#修正残差序列格式
dta_pred2=pd.read_csv('预测.CSV',header=0 ,index_col=0)#这里就是把上面的dta_pred打印出来复制到一个csv里面  因为直接用的话有bug不知道咋回事，好像是运算出来列的名字丢掉了
dta_pred2 = dta_pred2.dropna() #去除不完整项数据

pred_error = (dta_pred2 - dta3).dropna()#计算残差
lb,p_l= acorr_ljungbox(pred_error, lags = 16)#LB检验
print(p_l)
h = (p_l < 0.05).sum() #p值小于0.05，认为是非白噪声。
if h > 0:
    print (u'模型ARIMA(%s,1,%s)不符合白噪声检验'%(p,q))
else:
    print (u'模型ARIMA(%s,1,%s)符合白噪声检验' %(p,q))


#残差的自相关图
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = plot_acf(pred_error,ax=ax1)
ax2 = fig.add_subplot(212)
fig = plot_pacf(pred_error, ax=ax2)
fig.show()

#D-W检验
print(sm.stats.durbin_watson(pred_error))

#绘制qq图
fig = plt.figure(figsize=(12,8))
fig = qqplot(pred_error, line='q', fit=True)
fig.show()

###不同差分次数的精度
##print(ARIMA(dta, (p,0,q)).fit().bic)
##print(ARIMA(dta, (p,1,q)).fit().bic)
##
##print(ARIMA(dta, (p,0,q)).fit().aic)
##print(ARIMA(dta, (p,1,q)).fit().aic)
###print(ARIMA(dta, (p,2,q)).fit().aic) #MA系数不可逆
##
###差分比较
##fig1 = plt.figure(figsize=(8,7))
##ax1= fig1.add_subplot(211)
##diff1 = dta.diff(1)
##diff1.plot(ax=ax1)
##ax2= fig1.add_subplot(212)
##diff2 = dta.diff(2)
##diff2.plot(ax=ax2)
##fig1.show()
##
###差分后的自相关图
##dta= dta.diff(1)#1阶差分
##dta = dta.dropna()
##
##fig2 = plt.figure(figsize=(8,7))
##ax1=fig2.add_subplot(211)
##fig2 = plot_acf(dta,ax=ax1)
##ax2 = fig2.add_subplot(212)
##fig2 = plot_pacf(dta,ax=ax2)
##fig2.show()

预测

import warnings  # 忽略警告
warnings.filterwarnings('ignore')

import pandas as pd #表格和数据操作
import numpy as np
import matplotlib.pyplot as plt

from statsmodels.tsa.arima_model import ARMA
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.api import tsa
import statsmodels.api as sm

from itertools import product


dta0=pd.read_csv('时间序列_40.CSV',header=0)
dta0.index = pd.to_datetime(dta0['Date'])

data = pd.DataFrame()
data['date'] = ['2008/1/11','2008/2/6','2008/3/17','2008/4/13','2008/5/17','2008/6/15','2008/7/1','2008/7/12','2008/8/10','2008/9/14','2008/10/12','2008/11/16','2008/12/13','2009/1/19','2009/2/16','2009/3/13','2009/4/18','2009/5/16','2009/6/20','2009/7/11','2009/8/15','2009/9/19','2009/10/16','2009/11/14','2009/12/11','2010/1/15','2010/2/20','2010/3/13','2010/4/17','2010/5/15','2010/6/12','2010/7/16','2010/8/14','2010/9/18','2010/10/16','2010/11/19','2010/12/24','2011/1/21','2011/2/18','2011/3/19','2011/4/17','2011/5/15','2011/6/18','2011/7/16','2011/8/20','2011/9/24','2011/10/22','2011/11/19','2011/12/24','2012/1/14']
data['dy'] = [0.62,1.01,1.78,1.29,0.11,-0.35,-0.44,-0.3,-1.11,-1.78,-1.39,-0.94,-0.36,1.47,1.75,2.04,1.03,0.02,-0.59,-1.35,-2.14,-1.96,-1.46,-0.56,0.04,0.96,1.58,1.43,0.95,0.14,-0.3,-1.35,-1.6,-1.98,-1.58,-0.98,0.56,1.14,1.19,1.18,0.61,0.76,-0.66,-1.14,-1.35,-1.85,-0.95,-0.65,0.44,1.09]
data.index = pd.to_datetime(data['date'])


p=2
q=2

arima = ARIMA(dta0['dy'].dropna(), (p,1,q)).fit()

print(arima.summary())

dta_pred = arima.predict(typ = 'levels') #预测

###拟合
##fig1 = plt.figure(figsize=(12,8))
##plt.plot(dta0['dy'], color='green')
##plt.plot(dta_pred, color='yellow')
##fig1.show()

#模型预测
forecast_ts = arima.forecast(10)
fore = pd.DataFrame()
fore['date'] = ['2011-4-17','2011-5-15','2011-6-18','2011-7-16','2011-8-20','2011-9-24','2011-10-22','2011-11-19','2011-12-24','2012-1-14']
fore['result'] = pd.DataFrame(forecast_ts[0])
fore.index = pd.to_datetime(fore['date'])

print(fore['result'])

#绘制成果表
dta0['dy'].plot(color='blue', label='Original',figsize=(12,8))
dta_pred.plot(color='red', label='Predict',figsize=(12,8))
data.dy.plot(color='green', label='Original',figsize=(12,8))
fore.result.plot(color='yellow', label='Forecast',figsize=(12,8))

plt.show()

结果

上图

原始数据时序图

一次差分的时序图

2次差分的时序图

平稳性检验和白噪声检验

差分阶数	ADF检验（小于0.05可视作平稳）	LB检验（小于0.05可视作非白噪声序列）
0	1.1067394073215614e-08	4.4372756913599096e-09
1	2.3032333702927282e-11	0.0002755346487733437
2	1.2900060307415488e-08	1.1951415276363802e-07

自相关图

偏相关图

这个表和上面的图是用SPSS跑出来的，因为我用程序画的偏向关图很离谱

我用程序画的自相关偏向关图：

程序算出来p=q=2
然后算一下差分阶数

差分阶数	AIC	BIC
0	53.22929865693882	64.7014366895077
1	61.743467139654456	73.09438892831821
2	MA系数不可逆	MA系数不可逆

虽然0次差分最小，但是202模型算出来的残差不符合白噪声检验，这也很离谱，所以最后用的是212

然后是模型检验

先搞张图氵页数

在指数平滑模型下，观察ARIMA模型的残差是否是平均值为0且方差为常数的正态分布（服从零均值、方差不变的正态分布），同时也要观察连续残差是否（自）相关
先看看残差的自相关特性，嗯嗯 相当不错，很欣慰

LB检验
模型确立后，还需检验其残差序列是否为白噪声，如果不是白噪声，说明残差中还存在有用的信息，需要修改模型或进一步提取。利用模型和原始数据项做差得到残差项，对该残差序列作滞后值为16的LB检验的取得结果如下序列：[0.30226802 0.31656683 0.26285311 0.39561196 0.14861177 0.20026693 0.28606258 0.32825364 0.35109379 0.28906182 0.33058463 0.36657827 0.4070982 0.47689337 0.53658179 0.5805578 ]；可见残差中有用的信息较少，模型的拟合结果较好，可作为预测模型。

D-W检验
当DW值显著的接近于O或４时，则存在自相关性，而接近于２时，则不存在（一阶）自相关性。检验结果时2.26339224，说明残差的自相关性很小，可视作不存在自相关性

Q-Q图可以直观验证一组数据是否来自某个分布，或者验证某两组数据是否来自同一（族）分布。如果所有点基本上在一条直线上，我们可以说这两个分布是同一分布。所以根据qq图，我们得出结论目标数据大部分分布在一条直线的周围，可认为残差基本上服从正态分布（才怪）

到了预测环节，用前40组数据再建一个模型。为啥是前40个 ，因为前30个建不起来（吐血）

里面那个0.9的残差就很离谱，去掉之后感觉还可以

日期	实际值	预测值	残差	残差分析
2011/4/17	0.61	0.599297	0.010703	count 10
2011/5/15	0.76	-0.20582	0.96582	mean 0.241534
2011/6/18	-0.66	-1.01002	0.35002	std 0.358969
2011/7/16	-1.14	-1.597523	0.457523
2011/8/20	-1.35	-1.820037	0.470037	去极值后残差
2011/9/24	-1.85	-1.63406	-0.21594	count 9
2011/10/22	-0.95	-1.10831	0.15831	mean 0.161057
2011/11/19	-0.65	-0.400514	-0.24949	std 0.268526
2011/12/24	0.44	0.2889	0.1511
2012/1/14	1.09	0.77275	0.31725

实际观测值（绿色），拟合值（红色），和预测值（黄色）