基本概念

在这里插入图片描述

基本流程

导入模型

import sys
import os
import warnings
warnings.filterwarnings("ignore")
import pandas as pd
import numpy as np
from arch.unitroot import ADF
import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import style
style.use('ggplot')
import statsmodels.api as sm
import statsmodels.formula.api as smf
import statsmodels.tsa.api as smt
from statsmodels.tsa.stattools import adfuller 
from statsmodels.stats.diagnostic import acorr_ljungbox 
from statsmodels.graphics.api import qqplot
pd.set_option('display.float_format', lambda x: '%.5f' % x) 
np.set_printoptions(precision=5, suppress=True) 
"""中文显示问题"""
plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.sans-serif'] = ['SimHei']

在这里插入图片描述

1.平稳性检验

1）什么是序列平稳

在这里插入图片描述

2）检验是否平稳

方法一:观察序列图:若个体值要围绕序列均值上下波动，没有明显的上升或下降趋势，则是平稳序列。

在这里插入图片描述

方法二:观察自相关图:平稳的序列的自相关图（ACF）和偏相关图（PACF）不是拖尾就是截尾。

在这里插入图片描述
这张图可以看出，左图（自相关函数图）呈现出拖尾特征（2阶拖尾），而右图（偏自相关函数图）呈现出截尾特征（2阶截尾）

在这里插入图片描述
这张图可以看出，左图（自相关函数图）呈现出截尾特征（1阶截尾），而右图（偏自相关函数图）呈现出拖尾特征（1阶拖尾）

在这里插入图片描述
这张图可以看出，左图（自相关函数图）呈现出拖尾特征（2阶拖尾），右图（偏自相关函数图）也呈现出拖尾特征（1阶拖尾）

方法三：单位根检验：当一个时间序列的滞后算子多项式方程存在单位根时，我们认为该时间序列是非平稳的；反之，当该方程不存在单位根时，我们认为该时间序列是平稳的。

第一种方法：使用statsmodel的方法为：

from statsmodels.stats.diagnostic import unitroot_adf

unitroot_adf(df.pct_chg)

输出为：
这里包含了检验值、p-value、滞后阶数、自由度等信息。
我们看到了检验统计量为-14.46，远小于1%的临界值-3.47，即p值远小于0.01，因此我们拒绝原假设，认为该时间序列是平稳的。（这里原假设是存在单位根，即时间序列为非平稳的。)

第二种方法：使用arch的方法为：

from arch.unitroot import ADF

ADF(df.pct_chg)

其输出信息与statsmodel基本是一致的。

或

其他案例使用方法二

print("单位根检验:\n")
print(ADF(data.diff1.dropna()))

单位根检验:

   Augmented Dickey-Fuller Results   
=====================================
Test Statistic                 -3.156
P-value                         0.023
Lags                                0
-------------------------------------

Trend: Constant
Critical Values: -3.63 (1%), -2.95 (5%), -2.61 (10%)
Null Hypothesis: The process contains a unit root.
Alternative Hypothesis: The process is weakly stationary.

3）不平稳序列转为平稳序列

方法：差分变换

fig = plt.figure(figsize=(12,8))
ax1= fig.add_subplot(111)
diff1 = dta.diff(1)
diff1.plot(ax=ax1)

在这里插入图片描述
从下图可以看出，预测数列取对数并作一阶差分后的图形

选择其中一个检验方法，检验之后发现：
预测数列取对数并作一阶差分后的图形显示基本消除了长期趋势性的影响，趋于平稳化，满足ARIMA模型建模的基本要求。

2.白噪声检验

在这里插入图片描述

from statsmodels.stats.diagnostic import acorr_ljungbox
acorr_ljungbox(data.diff1.dropna(), lags = [i for i in range(1,12)],boxpierce=True)

(array([11.30402, 13.03896, 13.37637, 14.24184, 14.6937 , 15.33042,
        16.36099, 16.76433, 18.15565, 18.16275, 18.21663]),
 array([0.00077, 0.00147, 0.00389, 0.00656, 0.01175, 0.01784, 0.02202,
        0.03266, 0.03341, 0.05228, 0.07669]),
 array([10.4116 , 11.96391, 12.25693, 12.98574, 13.35437, 13.85704,
        14.64353, 14.94072, 15.92929, 15.93415, 15.9696 ]),
 array([0.00125, 0.00252, 0.00655, 0.01135, 0.02027, 0.03127, 0.04085,
        0.06031, 0.06837, 0.10153, 0.14226]))

3.模型定阶：计算自相关函数ACF 和偏自相关函数PACF

若你检验方法选择方法二，则你已经得到了ACF和PACF，但检验方法最好还是选择单位根检验

其中 L 是该模型下的最大似然, n 是数据数量， k 是模型的变量个数。

在这里插入图片描述

模型	ACF	PACF
AR	拖尾	截尾
MA	截尾	拖尾
ARMA	拖尾	拖尾

例子1:

1) 观察ACF和PACF图，判断是哪种模型

dta= dta.diff(1)#我们已经知道要使用一阶差分的时间序列，之前判断差分的程序可以注释掉
fig = plt.figure(figsize=(12,8))
ax1=fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(dta,lags=40,ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(dta,lags=40,ax=ax2)

在这里插入图片描述

2) AIC BIC选择其中最佳模型

arma_mod20 = sm.tsa.ARMA(dta,(7,0)).fit() # dta是差分后的序列
print(arma_mod20.aic,arma_mod20.bic,arma_mod20.hqic)
arma_mod30 = sm.tsa.ARMA(dta,(0,1)).fit()
print(arma_mod30.aic,arma_mod30.bic,arma_mod30.hqic)
arma_mod40 = sm.tsa.ARMA(dta,(7,1)).fit()
print(arma_mod40.aic,arma_mod40.bic,arma_mod40.hqic)
arma_mod50 = sm.tsa.ARMA(dta,(8,0)).fit()
print(arma_mod50.aic,arma_mod50.bic,arma_mod50.hqic)

在这里插入图片描述

例子2:

1) 观察ACF和PACF图，判断是哪种模型

在这里插入图片描述

自相关图拖尾或一阶截尾
偏自相关图一阶截尾,
所以我们可以建立ARIMA(1,1,0)、ARIMA(1,1,1)、ARIMA(0,1,1)模型。

2) AIC BIC选择其中最佳模型

arma_mod20 = sm.tsa.ARIMA(data["xt"],(1,1,0)).fit() # data["xt"]是没有差分的数据
arma_mod30 = sm.tsa.ARIMA(data["xt"],(0,1,1)).fit()
arma_mod40 = sm.tsa.ARIMA(data["xt"],(1,1,1)).fit()
values = [[arma_mod20.aic,arma_mod20.bic,arma_mod20.hqic],[arma_mod30.aic,arma_mod30.bic,arma_mod30.hqic],[arma_mod40.aic,arma_mod40.bic,arma_mod40.hqic]]
df = pd.DataFrame(values,index=["AR(1,1,0)","MA(0,1,1)","ARMA(1,1,1)"],columns=["AIC","BIC","hqic"])
df

	AIC	BIC	hqic
AR(1,1,0)	253.09159	257.84215	254.74966
MA(0,1,1)	251.97340	256.72396	253.63147
ARMA(1,1,1)	254.09159	258.84535	259.74966

4.参数估计

from statsmodels.tsa.arima_model import ARIMA
model = ARIMA(data["xt"], order=(0,1,1))
result = model.fit()
print(result.summary())

                             ARIMA Model Results                              
==============================================================================
Dep. Variable:                   D.xt   No. Observations:                   36
Model:                 ARIMA(0, 1, 1)   Log Likelihood                -122.987
Method:                       css-mle   S.D. of innovations              7.309
Date:                Tue, 22 Dec 2020   AIC                            251.973
Time:                        09:11:55   BIC                            256.724
Sample:                    01-01-1953   HQIC                           253.631
                         - 01-01-1988                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.9956      2.014      2.481      0.013       1.048       8.943
ma.L1.D.xt     0.6710      0.165      4.071      0.000       0.348       0.994
                                    Roots                                    
=============================================================================
                  Real          Imaginary           Modulus         Frequency
-----------------------------------------------------------------------------
MA.1           -1.4902           +0.0000j            1.4902            0.5000
-----------------------------------------------------------------------------

5.模型检验

在这里插入图片描述

1) 参数的显著性检验

在这里插入图片描述

2）模型的显著性检验

在这里插入图片描述

resid = result.resid#残差
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
fig = qqplot(resid, line='q', ax=ax, fit=True)

在这里插入图片描述

6.模型预测

pred = result.predict('1988', '1990',dynamic=True, typ='levels')
print (pred)

1988-01-01   278.35527
1989-01-01   283.35088
1990-01-01   288.34649
Freq: AS-JAN, dtype: float64

plt.figure(figsize=(12, 8))
plt.xticks(rotation=45)
plt.plot(pred)
plt.plot(data.xt)
plt.show()

在这里插入图片描述

7.预测结果分析

在这里插入图片描述
平均绝对误差：

平均相对误差：

预测均方差：

其中，y_i为序列的实际数据， $\widehat{y}_i$ 为模型预测数据。

例如：
在这里插入图片描述

引用：
https://blog.csdn.net/qixizhuang/article/details/86716149
https://zhidao.baidu.com/question/288273884.html
https://zhuanlan.zhihu.com/p/343617191
https://blog.csdn.net/qq_43613793/article/details/109908418
https://blog.csdn.net/qq_45176548/article/details/116771331
https://blog.csdn.net/weixin_41555408/article/details/117218457
https://blog.csdn.net/weixin_41636030/article/details/89138638
https://www.zhihu.com/question/48447779/answer/112330897
https://blog.csdn.net/u010414589/article/details/49622625
https://wiki.mbalib.com/wiki/ARIMA%E6%A8%A1%E5%9E%8B
https://blog.csdn.net/qq_45176548/article/details/111504846