python-arima模型statsmodels库实现-有数据集（续）-statsmodels-0.9.0版本

python-arima模型statsmodels库实现-有数据集（续）

这篇博客是上一篇python-arima模型statsmodels库实现的续集，上一篇采用的statsmodels版本应该要高一点，如果使用低版本的statsmodels代码会有bug，这一篇则是针对statsmodels-0.9.0版本的代码。

代码如下：

#coding=gbk
import  numpy  as np
import pandas as pd
import os
from numpy import NaN
from numpy import nan
import matplotlib.pyplot as plt
import statsmodels.api as sm     #acf,pacf图
from statsmodels.tsa.stattools import adfuller  #adf检验
from pandas.plotting import autocorrelation_plot
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.stats.diagnostic import acorr_ljungbox

import statsmodels.api as sm
import matplotlib as mpl
path="E:/data/china_data.xlsx"
# 为了控制计算量，我们限制AR最大阶不超过6，MA最大阶不超过4。





plt.style.use('fivethirtyeight')
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['font.serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
df=pd.read_excel(path)
#print(df)
#help(df)

#for index, row in df.iterrows():

df=df.replace(NaN, "null")
#  print(index, row)
print(df)
def  f(column):
    r=0
    inde1=0
    index2=len(column)-1
    for i in range(len(column)):
     #   print(column[len(column)-i-1])
        if   column[len(column)-i-1] is "null" and r==1:
            index2=i
            return index1,index2
       
        if   column[len(column)-i-1]!= "null" and r==0:
             index1=i
             r=1
    return index1,index2
           
        
#df['时间(年)']=pd.to_datetime(df['时间(年)'])

print(df.columns)
print(df[df.columns[0]])
indexz=df.columns[0]


def adf_test(data):#小于0.05则是平稳序列
   # print("data:",data.values)
    data_z=np.array(list(data.values))
    #print(data_z.reshape(-1,))
    
   
    t = adfuller(data_z.reshape(-1,))
    
    print("p-value:",t[1])
def  box_pierce_test(data):#小于0.05，不是白噪声序列
    print(acorr_ljungbox(data, lags=1)) 

def  stability_judgment(data):
    fig = plt.figure(figsize=(12,8))
    ax1=fig.add_subplot(211)
    fig = sm.graphics.tsa.plot_acf(data,lags=5,ax=ax1)
    ax2 = fig.add_subplot(212)
    fig = sm.graphics.tsa.plot_pacf(data,lags=5,ax=ax2)
    plt.show()


def  model_fit(data,df,index,length,index1,index2):
    data_diff=df[["时间(年)",index]][length-index2:length-index1]
  #  sm.tsa.arma_order_select_ic(data_diff,max_ar=6,max_ma=4,ic='aic')['aic_min_order']  # AIC


    #对模型进行定阶
    pmax = int(len(data) / 10)    #一般阶数不超过 length /10
    qmax = int(len(data) / 10)
    if  pmax>4:pmax=6
    if  qmax>4:qmax=4
    bic_matrix = []
    print("data",data)

   # help(sm.tsa.arima.ARIMA)
    for p in range(pmax +1):
        temp= []
        for q in range(qmax+1):
            try:
         #  ARIMA(train_data, order=(1,1,1))

          # print(sm.tsa.arima.ARIMA(data,order=(p,1,q)).fit())
                  temp.append(sm.tsa.ARIMA(data,order=(p,1,q)).fit().bic)
                #  print(temp)
            except:
                temp.append(None)
           # temp.append(sm.tsa.arima.ARIMA(data,order=(p,1,q)).fit().bic)
            bic_matrix.append(temp)
             
    bic_matrix = pd.DataFrame(bic_matrix)   #将其转换成Dataframe 数据结构
    print("bic_matrix",bic_matrix)
    p,q = bic_matrix.stack().astype(float).idxmin()   #先使用stack 展平， 然后使用 idxmin 找出最小值的位置
    print(u'BIC 最小的p值 和 q 值：%s,%s' %(p,q))  #  BIC 最小的p值 和 q 值：0,1
   
    model = sm.tsa.ARIMA(data, order=(p,1,q)).fit()
    model.summary()        #生成一份模型报告
    predictions_ARIMA_diff = pd.Series(model.fittedvalues, copy=True)
    print(predictions_ARIMA_diff)
    model.forecast(5)   #为未来5天进行预测， 返回预测结果， 标准误差， 和置信区间


    


for index, column in df.iteritems():
    if index==indexz:
        continue
    index1,index2 =f(column)
    length=len(column)
   # print("index1 index2：",index1,index2)

  #  print(column[length-index2-1:length-index1])
    print(index)
    df[index]=df[index].replace( "null",0)
    df[index].astype('float')
    
    df[str(index)+"diff1"]=df[index].diff(1)
    df[str(index)+"diff2"]=df[index+"diff1"].diff(1)
    # 一阶差分还原
    # tmpdata2：原数据
    # pred：一阶差分后的预测数据
    #df_shift = tmpdata2['ecpm_tomorrow'].shift(1)
    #predict = pred.add(df_shift)
    # predict = pred + df_shift

   # print(index2-index1)
    #print(df[["时间(年)",index]][length-index2:length-index1])
    adf_test(df[[index]][length-index2:length-index1])
    box_pierce_test(df[[index]][length-index2:length-index1])
    model_fit(df[[index]][length-index2:length-index1],df,index,length,index1,index2)
   
  ## model_fit(data,p,q)


    stability_judgment(df[[index]][length-index2:length-index1])

    stability_judgment(df[[str(index)+"diff1"]][length-index2:length-index1])
  #  stability_judgment(df[[str(index)+"diff2"]][length-index2:length-index1])

    plt.plot(df[["时间(年)"]][length-index2:length-index1],df[[index]][length-index2:length-index1],label="diff0")
    plt.plot(df[["时间(年)"]][length-index2:length-index1],df[[str(index)+"diff1"]][length-index2:length-index1],label="diff1")
 #   plt.plot(df[["时间(年)"]][length-index2:length-index1],df[[str(index)+"diff2"]][length-index2:length-index1],label="diff2")

   # df[["时间(年)",index]][length-index2:length-index1].plot(x=indexz,y=index,figsize=(9,9))
    plt.xlabel("时间（年）")
    plt.ylabel(index)
    plt.legend()
    plt.show()
 




os.system("pause")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173

运行结果如下：

大家可

大家可以学习一下哈。