专题二 四种基础的时间序列模型
专题二 四种基础的时间序列模型
本文讲解了四种最简单的时间序列模型，从定义、确立模型、模型应用、结果分析出发，通过阅读可以迅速上手简易的时间序列模型。
2.1 自回归(AR)模型
AR模型的表达式如下：
其中：
表示时间序列在时间点
的观测值。
是常数项（截距）。
是模型的阶数，表示使用过去
是对应于时间点
的系数，表示在时间点
是白噪声（或误差项），表示模型无法解释的随机波动。
AR模型的预测依赖于过去的观测值和模型中的参数。模型的阶数 决定了需要考虑多少个过去时间点的观测值。
AR模型可以用来探索时间序列数据中的趋势和周期性，帮助分析历史数据并进行未来的预测。为了确定AR模型的阶数 和参数
，常常会使用统计方法如最小二乘法、信息准则（如AIC、BIC）或者其他模型选择技术。
注意：AR模型假设时间序列是平稳的，即其均值和方差在时间上保持不变。如果时间序列数据不满足平稳性假设，需要进行平稳性转换或考虑其他时间序列模型，如自回归移动平均模型（ARMA）或自回归积分移动平均模型（ARIMA）。
# AP(p)
# 导入需要的包
import pandas as pd
import tushare as ts #获取股票数据需要安装的库
import numpy as np
import matplotlib.pyplot as plt #绘图
import seaborn as sns #seaborn画出的图更好看，且代码更简单
sns.set(color_codes=True) #seaborn设置背景
data=pd.read_excel('股价数据.xlsx')
data.set_index('date', inplace=True) #设定日期为索引
r2=np.log(data['close'])-np.log(data['close'].shift(1)) #计算对数收益率
r2=r2.dropna()
r2.head()
date
2020-01-03 -0.051530
2020-01-06 -0.000586
2020-01-07 0.016870
2020-01-08 -0.006484
2020-01-09 0.014713
Name: close, dtype: float64
2.1.1 AR
模型的定阶
一般有两种方法来决定p：
第一种：利用偏相关函数(Partial Auto Correlation Function,PACF)
第二种：利用信息准则函数
1. 偏相关函数判断p
对于偏相关函数(PACF)的介绍，见1.4.2，这里不赘述：
AR§序列的样本偏相关函数是p步截尾的。
所谓截尾，就是快速收敛应该是快速的降到几乎为0或者在置信区间以内。
PACF图形是描述时间序列数据的偏自相关结构的图形。在PACF图中，横轴也表示滞后阶数（lag），纵轴表示偏相关系数。PACF图的截尾点之后的值也接近于零，截尾点前的值通常显示出截尾点之后的指数衰减。
截尾点之后第一个超过显著性水平的偏相关系数通常用来确定AR模型的阶数。
具体我们看下面的例子
from statsmodels.graphics.tsaplots import plot_pacf as PACF #偏自相关图
fig = PACF(r2,lags = 30) #使用对数收益率序列
plt.show()
可以看出，1阶之后基本接近0，定阶1或者6或者26
2. 信息准则— AIC、BIC、HQ
给出的结果!=0，!=1 使用其，否则根据PACF图形的结果进行尝试
现在有以上这么多可供选择的模型，我们通常采用AIC法则。我们知道：增加自由参数的数目提高了拟合的优良性，AIC鼓励数据拟合的优良性但是尽量避免出现过度拟合(Overfitting)的情况。所以优先考虑的模型应是AIC值最小的那一个。赤池信息准则的方法是寻找可以最好地解释数据但包含最少自由参数的模型。不仅仅包括AIC准则，目前选择模型常用如下准则：
中文名字：赤池信息量 akaike information criterion
中文名字：贝叶斯信息量 bayesian information criterion
汉南-奎恩准则 hannan-quinn criterion
下面我们来测试下3种准则下确定的p，仍然用茅台日收益率序列。
from statsmodels.tsa.arima.model import ARIMA
#14.0 statamdels.org
aicList = []
bicList = []
hqicList = []
for i in range(1,10): #从1阶开始算
od = (i,0,0) # 这里使用了ARMA模型，order 代表了模型的(p,q)值，我们令q始终为0，就只考虑了AR情况。
tempModel = ARIMA(r2,order=od).fit()
aicList.append(tempModel.aic)
bicList.append(tempModel.bic)
hqicList.append(tempModel.hqic)
plt.figure(figsize=(15,6))
plt.plot(aicList,'r',label='aic value')
plt.plot(bicList,'b',label='bic value')
plt.plot(hqicList,'k',label='hqic value')
plt.legend(loc=0)
直接输出AIC确定的结果
from scipy import stats
import statsmodels.api as sm # 统计相关的库
print(sm.tsa.arma_order_select_ic(r2,max_ar=6,max_ma=0,ic='aic')['aic_min_order'])
(0, 0)
因此我们建立AR(26)模型
from scipy import stats
import statsmodels.api as sm # 统计相关的库
from statsmodels.tsa.ar_model import AutoReg
temp = np.array(r2) # 载入收益率序列
model =AutoReg(temp,lags=[1,6,25,26])
res = model.fit()
out = 'AIC: {0:0.3f}, HQIC: {1:0.3f}, BIC: {2:0.3f}'
print(out.format(res.aic, res.hqic, res.bic))
print(res.summary())#输出模型结果
plt.rcParams['font.sans-serif'] = ['simhei'] #字体为黑体，正常显示中文字体
plt.rcParams['axes.unicode_minus'] = False #正常显示负号
plt.figure(figsize=(10,4))
plt.plot(temp,'b',label='对数收益率') #核心画图代码 原始收益率序列
plt.plot(res.fittedvalues, 'r',label='AR model')#模型拟合结果
plt.legend()
AIC: -3400.990, HQIC: -3390.432, BIC: -3373.675
AutoReg Model Results
==============================================================================
Dep. Variable: y No. Observations: 727
Model: Restr. AutoReg(26) Log Likelihood 1706.495
Method: Conditional MLE S.D. of innovations 0.021
Date: Wed, 09 Aug 2023 AIC -3400.990
Time: 22:19:45 BIC -3373.675
Sample: 26 HQIC -3390.432
727
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 0.0008 0.001 1.054 0.292 -0.001 0.002
y.L1 0.0029 0.037 0.078 0.938 -0.070 0.076
y.L6 -0.0426 0.037 -1.142 0.253 -0.116 0.031
y.L25 0.0546 0.037 1.465 0.143 -0.018 0.128
y.L26 -0.1368 0.037 -3.680 0.000 -0.210 -0.064
Roots
==============================================================================
Real Imaginary Modulus Frequency
------------------------------------------------------------------------------
AR.1 -1.0603 -0.1318j 1.0685 -0.4803
AR.2 -1.0603 +0.1318j 1.0685 0.4803
AR.3 -0.9941 -0.3837j 1.0656 -0.4414
AR.4 -0.9941 +0.3837j 1.0656 0.4414
AR.5 -0.8740 -0.6099j 1.0657 -0.4030
AR.6 -0.8740 +0.6099j 1.0657 0.4030
AR.7 -0.7048 -0.8059j 1.0706 -0.3644
AR.8 -0.7048 +0.8059j 1.0706 0.3644
AR.9 -0.4882 -0.9566j 1.0740 -0.3251
AR.10 -0.4882 +0.9566j 1.0740 0.3251
AR.11 -0.2412 -1.0461j 1.0735 -0.2861
AR.12 -0.2412 +1.0461j 1.0735 0.2861
AR.13 0.0147 -1.0741j 1.0742 -0.2478
AR.14 0.0147 +1.0741j 1.0742 0.2478
AR.15 0.2711 -1.0459j 1.0805 -0.2096
AR.16 0.2711 +1.0459j 1.0805 0.2096
AR.17 0.5188 -0.9556j 1.0873 -0.1708
AR.18 0.5188 +0.9556j 1.0873 0.1708
AR.19 0.7355 -0.8040j 1.0896 -0.1321
AR.20 0.7355 +0.8040j 1.0896 0.1321
AR.21 0.9045 -0.6082j 1.0900 -0.0942
AR.22 0.9045 +0.6082j 1.0900 0.0942
AR.23 1.0925 -0.1318j 1.1004 -0.0191
AR.24 1.0925 +0.1318j 1.1004 0.0191
AR.25 1.0253 -0.3830j 1.0945 -0.0569
AR.26 1.0253 +0.3830j 1.0945 0.0569
------------------------------------------------------------------------------
2.1.2 模型的检验
根据AR模型定义式，如果模型是充分的，其残差序列应该是白噪声，根据我们第一节里1.6.3介绍的混成检验，可以用来检验残差与白噪声的接近程度。
我们先求出残差序列：
delta = res.fittedvalues - temp[26:] # 残差
plt.figure(figsize=(10,6))
plt.plot(delta,'r',label=' residual error')
plt.legend(loc=0)
然后我们检查它是不是接近白噪声序列
acf,q,p = sm.tsa.acf(delta,nlags=10,qstat=True) ## 计算自相关系数 及p-value
out = np.c_[range(1,11), acf[1:], q, p]
output=pd.DataFrame(out, columns=['lag', "AC", "Q", "P-value"])
output = output.set_index('lag')
output
| AC | Q | P-value | |
|---|---|---|---|
| lag | |||
| 1.0 | 0.000429 | 0.000130 | 0.990913 |
| 2.0 | 0.019465 | 0.267240 | 0.874922 |
| 3.0 | -0.029681 | 0.889214 | 0.828030 |
| 4.0 | 0.050446 | 2.688470 | 0.611234 |
| 5.0 | -0.012696 | 2.802594 | 0.730388 |
| 6.0 | -0.003114 | 2.809472 | 0.832352 |
| 7.0 | 0.039262 | 3.904103 | 0.790754 |
| 8.0 | 0.030937 | 4.584702 | 0.800900 |
| 9.0 | -0.036161 | 5.515918 | 0.787218 |
| 10.0 | 0.030898 | 6.196767 | 0.798470 |
观察p-value可知，该序列可以认为没有相关性，近似得可以认为残差序列接近白噪声。说明模型是显著的。
2.1.3 拟合优度
我们使用下面的统计量来衡量拟合优度：
但是，对于一个给定的数据集，是用参数个数的非降函数，为了克服该缺点，推荐使用调整后的
:
它的值在0-1之间，越接近1，拟合效果越好。
score = 1 - delta.var()/temp[26:].var()
print(score)
0.023891486642034998
可以看出，模型的拟合程度并不好，当然，这并不重要，也许是这个序列并不适合用AR模型拟合。
2.2 滑动平均
模型
滑动平均（Moving Average，MA）模型是时间序列分析中的一种基本模型。它是ARIMA模型的一部分，用于描述时间序列中的随机波动。
MA模型的基本思想是将时间序列数据表示为滞后项（延迟项）的线性组合。具体来说，在MA模型中，当前时刻的观测值被解释为过去时刻的随机噪声和滞后项的线性组合。
MA(q)模型的表示为：
其中，
是常数项，表示序列的均值；
是MA模型的参数，表示滞后项的权重。
MA(q)模型的特点是，随着滞后阶数 的增加，模型中的滞后项增多，从而增强了对过去噪声的依赖程度。MA模型适用于一些在趋势和季节性变化不明显的时间序列数据。
在实际应用中，通常会使用ACF（自相关函数）和PACF（偏自相关函数）来辅助确定MA模型的阶数 。ACF和PACF可以帮助分析滞后项之间的相关性，从而帮助选择适合的MA模型。然后，可以使用最小二乘估计等方法来估计MA模型的参数。
MA模型与ARIMA模型的另外两个部分——自回归（AR）模型和差分（I）模型结合在一起，构成了常见的ARIMA模型，用于时间序列的建模和预测。
2.2.1 MA模型的性质
1. 平稳性
MA模型总是弱平稳的，因为他们是白噪声序列（残差序列）的有限线性组合。因此，根据弱平稳的性质可以得出两个结论：
2. 自相关函数
对q阶的MA模型，其自相关函数ACF总是q步截尾的。 因此MA(q)序列只与其前q个延迟值线性相关，从而它是一个“有限记忆”的模型。
这一点可以用来确定模型的阶次.
3.可逆性
当满足可逆条件的时候，MA(q)模型可以改写为AR(p)模型。这里不进行推导，给出1阶和2阶MA的可逆性条件。
1阶：
2阶：
2.2.2 MA的阶次判定
我们通常利用上面介绍的第二条性质：MA(q)模型的ACF函数q步截尾来判断模型阶次。示例如下：
仍然使用收益率序列来进行，画出序列的ACF：
from statsmodels.graphics.tsaplots import plot_acf as ACF #自相关图
fig = ACF(r2,lags = 30)
plt.show()
从置信区间来看，1阶时候已经在置信区间内部了，从截尾来看应当在26阶。因此判定MA模型阶次为1或者26阶。
信息准则定阶：
print(sm.tsa.arma_order_select_ic(r2,max_ar=0,max_ma=30,ic='aic')['aic_min_order'])
(0, 0)
2.2.3 建模和预测
由于sm.tsa中没有单独的MA模块，我们利用ARIMA模块，只要将其中ARIMA(p,d,q)中的p,d设置为0即可。
from statsmodels.tsa.arima.model import ARIMA
temp = np.array(r2) # 载入收益率序列
model =ARIMA(temp,order=(0, 0, 26))
res = model.fit()
print(res.summary())
plt.rcParams['font.sans-serif'] = ['simhei'] #字体为黑体
plt.rcParams['axes.unicode_minus'] = False #正常显示负号
plt.figure(figsize=(10,4))
plt.plot(temp,'b',label='对数收益率')
plt.plot(res.fittedvalues, 'r',label='MA model')
plt.legend()
SARIMAX Results
==============================================================================
Dep. Variable: y No. Observations: 727
Model: ARIMA(0, 0, 26) Log Likelihood 1781.241
Date: Wed, 09 Aug 2023 AIC -3506.481
Time: 22:43:21 BIC -3377.991
Sample: 0 HQIC -3456.899
- 727
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 0.0007 0.001 0.848 0.397 -0.001 0.002
ma.L1 0.0183 0.037 0.495 0.621 -0.054 0.091
ma.L2 0.0134 0.036 0.374 0.709 -0.057 0.084
ma.L3 -0.0304 0.034 -0.888 0.374 -0.097 0.037
ma.L4 0.0474 0.039 1.209 0.227 -0.029 0.124
ma.L5 -0.0065 0.036 -0.179 0.858 -0.078 0.065
ma.L6 -0.0575 0.034 -1.682 0.093 -0.124 0.009
ma.L7 0.0436 0.036 1.195 0.232 -0.028 0.115
ma.L8 0.0323 0.036 0.886 0.376 -0.039 0.104
ma.L9 -0.0431 0.038 -1.127 0.260 -0.118 0.032
ma.L10 0.0217 0.033 0.665 0.506 -0.042 0.086
ma.L11 -0.0494 0.037 -1.336 0.181 -0.122 0.023
ma.L12 0.0521 0.038 1.367 0.172 -0.023 0.127
ma.L13 -0.0349 0.038 -0.924 0.356 -0.109 0.039
ma.L14 0.0311 0.038 0.809 0.418 -0.044 0.107
ma.L15 -0.0483 0.038 -1.283 0.199 -0.122 0.025
ma.L16 0.0453 0.039 1.173 0.241 -0.030 0.121
ma.L17 -0.0086 0.038 -0.226 0.821 -0.083 0.066
ma.L18 -0.0532 0.041 -1.283 0.200 -0.135 0.028
ma.L19 0.0537 0.041 1.310 0.190 -0.027 0.134
ma.L20 0.0289 0.042 0.694 0.488 -0.053 0.110
ma.L21 -0.0045 0.037 -0.123 0.902 -0.077 0.068
ma.L22 0.0308 0.040 0.764 0.445 -0.048 0.110
ma.L23 0.0174 0.038 0.458 0.647 -0.057 0.092
ma.L24 -0.0262 0.041 -0.639 0.523 -0.107 0.054
ma.L25 0.0366 0.041 0.884 0.377 -0.045 0.118
ma.L26 -0.1432 0.039 -3.633 0.000 -0.220 -0.066
sigma2 0.0004 1.9e-05 22.885 0.000 0.000 0.000
===================================================================================
Ljung-Box (L1) (Q): 0.00 Jarque-Bera (JB): 76.16
Prob(Q): 0.98 Prob(JB): 0.00
Heteroskedasticity (H): 0.97 Skew: -0.15
Prob(H) (two-sided): 0.78 Kurtosis: 4.56
===================================================================================
2.2.4 模型的检验和拟合优度
delta = res.fittedvalues - temp # 残差
plt.figure(figsize=(10,6))
plt.plot(delta,'r',label=' residual error')
plt.legend(loc=0)
acf,q,p = sm.tsa.acf(delta,nlags=10,qstat=True) ## 计算自相关系数 及p-value
out = np.c_[range(1,11), acf[1:], q, p]
output=pd.DataFrame(out, columns=['lag', "AC", "Q", "P-value"])
output = output.set_index('lag')
output
| AC | Q | P-value | |
|---|---|---|---|
| lag | |||
| 1.0 | 0.001044 | 0.000795 | 0.977507 |
| 2.0 | -0.003513 | 0.009819 | 0.995103 |
| 3.0 | 0.004257 | 0.023087 | 0.999073 |
| 4.0 | -0.002255 | 0.026813 | 0.999911 |
| 5.0 | -0.003205 | 0.034351 | 0.999989 |
| 6.0 | 0.002258 | 0.038099 | 0.999999 |
| 7.0 | -0.002213 | 0.041703 | 1.000000 |
| 8.0 | -0.001425 | 0.043201 | 1.000000 |
| 9.0 | 0.000063 | 0.043204 | 1.000000 |
| 10.0 | 0.004830 | 0.060450 | 1.000000 |
score = 1 - delta.var()/temp.var() #计算拟合优度
print(score)
0.05142262438717027
效果也比较一般，说明日收益率序列的预测还要使用其他模型才更准确。
2.3 ARMA模型
ARMA（AutoRegressive Moving Average）模型是一种常用的时间序列分析模型，结合了自回归（AR）模型和滑动平均（MA）模型的特性。ARMA模型用于描述时间序列数据中的自相关性和随机波动。
ARMA(p, q)模型的表示为：
其中，
是AR模型的参数，表示滞后项的权重；
ARMA模型适用于具有自相关性和滞后项之间相关性的时间序列数据。它在拟合时间序列数据方面比单独使用AR或MA模型更加灵活，因为它可以同时考虑多个滞后项的影响。
在实际应用中，选择ARMA模型的合适阶数 和
是很重要的一步。通常可以使用ACF（自相关函数）和PACF（偏自相关函数）来辅助确定ARMA模型的阶数。然后，可以使用最小二乘估计等方法来估计ARMA模型的参数。
ARMA模型是时间序列建模和预测中的重要工具，但对于一些时间序列数据，可能需要更复杂的模型，例如ARIMA（包含差分项）或季节性ARIMA（SARIMA）模型，以更好地捕捉数据的特征。
特点：ARMA模型的应用对象应该为平稳序列！我们下面的步骤都是建立在假设原序列平稳的条件下的~
2.3.1 ARMA模型的定阶
1. PACF、ACF 判断模型阶次
我们通过观察PACF和ACF截尾，分别判断p、q的值。
from statsmodels.graphics.tsaplots import plot_acf as ACF #自相关图
fig = ACF(r2,lags = 30)
plt.show()
from statsmodels.graphics.tsaplots import plot_pacf as PACF #偏自相关图
fig = PACF(r2,lags = 30)
plt.show()
结合前文关于AR、MA模型的定阶情况，ARMA模型的阶次应该定位(26,26)
2. 信息准则定阶（AIC/BIC/HQ）
sm.tsa.arma_order_select_ic(r2,max_ar=6,max_ma=4,ic='aic')['aic_min_order'] # AIC
(0, 0)
2.3.2 模型的建立及预测
from statsmodels.tsa.arima.model import ARIMA
temp = np.array(r2) # 载入收益率序列
model =ARIMA(temp,order=(26, 0, 26))
res = model.fit()
print(res.summary())
plt.rcParams['font.sans-serif'] = ['simhei'] #字体为黑体
plt.rcParams['axes.unicode_minus'] = False #正常显示负号
plt.figure(figsize=(10,4))
plt.plot(temp,'b',label='对数收益率')
plt.plot(res.fittedvalues, 'r',label='ARMA model')
plt.legend()
SARIMAX Results
==============================================================================
Dep. Variable: y No. Observations: 727
Model: ARIMA(26, 0, 26) Log Likelihood 1783.552
Date: Wed, 09 Aug 2023 AIC -3459.105
Time: 22:46:23 BIC -3211.303
Sample: 0 HQIC -3363.482
- 727
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 0.0007 0.001 0.797 0.425 -0.001 0.002
ar.L1 0.0322 0.445 0.072 0.942 -0.841 0.905
ar.L2 -0.0384 0.418 -0.092 0.927 -0.858 0.782
ar.L3 0.0342 0.429 0.080 0.936 -0.806 0.874
ar.L4 -0.0198 0.393 -0.050 0.960 -0.790 0.751
ar.L5 4.234e-05 0.431 9.82e-05 1.000 -0.845 0.845
ar.L6 -0.0192 0.425 -0.045 0.964 -0.851 0.813
ar.L7 0.0036 0.442 0.008 0.994 -0.863 0.870
ar.L8 0.0297 0.378 0.078 0.937 -0.712 0.771
ar.L9 -0.0128 0.374 -0.034 0.973 -0.745 0.720
ar.L10 0.0385 0.390 0.099 0.921 -0.725 0.802
ar.L11 -0.0029 0.410 -0.007 0.994 -0.807 0.801
ar.L12 0.0761 0.365 0.209 0.835 -0.639 0.792
ar.L13 -0.0970 0.346 -0.280 0.779 -0.776 0.582
ar.L14 0.0021 0.375 0.005 0.996 -0.733 0.737
ar.L15 -0.0485 0.368 -0.132 0.895 -0.770 0.673
ar.L16 0.0229 0.336 0.068 0.946 -0.635 0.680
ar.L17 0.0362 0.346 0.105 0.917 -0.642 0.715
ar.L18 -0.0328 0.337 -0.097 0.922 -0.693 0.627
ar.L19 -0.0520 0.370 -0.141 0.888 -0.777 0.673
ar.L20 -0.0141 0.341 -0.041 0.967 -0.683 0.655
ar.L21 0.0619 0.345 0.179 0.858 -0.614 0.738
ar.L22 -0.0590 0.348 -0.170 0.865 -0.741 0.623
ar.L23 0.0952 0.327 0.291 0.771 -0.546 0.737
ar.L24 0.0410 0.317 0.129 0.897 -0.581 0.663
ar.L25 0.0477 0.321 0.149 0.882 -0.582 0.677
ar.L26 -0.0397 0.291 -0.136 0.892 -0.610 0.531
ma.L1 -0.0183 0.439 -0.042 0.967 -0.878 0.841
ma.L2 0.0473 0.421 0.112 0.911 -0.779 0.873
ma.L3 -0.0523 0.424 -0.123 0.902 -0.884 0.779
ma.L4 0.0569 0.398 0.143 0.886 -0.723 0.837
ma.L5 -0.0104 0.436 -0.024 0.981 -0.864 0.843
ma.L6 -0.0355 0.424 -0.084 0.933 -0.867 0.796
ma.L7 0.0372 0.442 0.084 0.933 -0.830 0.904
ma.L8 0.0128 0.382 0.033 0.973 -0.735 0.761
ma.L9 -0.0346 0.377 -0.092 0.927 -0.774 0.705
ma.L10 -0.0146 0.400 -0.036 0.971 -0.799 0.770
ma.L11 -0.0526 0.415 -0.127 0.899 -0.865 0.760
ma.L12 -0.0101 0.366 -0.028 0.978 -0.727 0.707
ma.L13 0.0408 0.358 0.114 0.909 -0.661 0.743
ma.L14 0.0436 0.382 0.114 0.909 -0.704 0.792
ma.L15 -0.0133 0.388 -0.034 0.973 -0.775 0.748
ma.L16 0.0324 0.349 0.093 0.926 -0.652 0.717
ma.L17 -0.0465 0.352 -0.132 0.895 -0.736 0.643
ma.L18 -0.0210 0.340 -0.062 0.951 -0.686 0.644
ma.L19 0.0977 0.369 0.265 0.791 -0.625 0.821
ma.L20 0.0404 0.346 0.117 0.907 -0.637 0.718
ma.L21 -0.0478 0.353 -0.135 0.892 -0.740 0.644
ma.L22 0.0767 0.351 0.219 0.827 -0.611 0.764
ma.L23 -0.0594 0.340 -0.175 0.861 -0.726 0.607
ma.L24 -0.0883 0.330 -0.268 0.789 -0.735 0.559
ma.L25 0.0096 0.332 0.029 0.977 -0.640 0.659
ma.L26 -0.1122 0.315 -0.356 0.722 -0.729 0.505
sigma2 0.0004 2.12e-05 20.383 0.000 0.000 0.000
===================================================================================
Ljung-Box (L1) (Q): 0.03 Jarque-Bera (JB): 74.71
Prob(Q): 0.86 Prob(JB): 0.00
Heteroskedasticity (H): 0.97 Skew: -0.14
Prob(H) (two-sided): 0.81 Kurtosis: 4.54
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
2.3.3模型的检验和拟合优度
delta = res.fittedvalues - temp # 残差
plt.figure(figsize=(10,6))
plt.plot(delta,'r',label=' residual error')
plt.legend(loc=0)
acf,q,p = sm.tsa.acf(delta,nlags=10,qstat=True) ## 计算自相关系数 及p-value
out = np.c_[range(1,11), acf[1:], q, p]
output=pd.DataFrame(out, columns=['lag', "AC", "Q", "P-value"])
output = output.set_index('lag')
output
| AC | Q | P-value | |
|---|---|---|---|
| lag | |||
| 1.0 | 0.006433 | 0.030206 | 0.862024 |
| 2.0 | -0.000137 | 0.030220 | 0.985004 |
| 3.0 | -0.007029 | 0.066384 | 0.995541 |
| 4.0 | 0.009800 | 0.136786 | 0.997765 |
| 5.0 | -0.001881 | 0.139383 | 0.999633 |
| 6.0 | 0.002682 | 0.144672 | 0.999940 |
| 7.0 | -0.000363 | 0.144769 | 0.999992 |
| 8.0 | -0.012086 | 0.252436 | 0.999990 |
| 9.0 | 0.002153 | 0.255856 | 0.999998 |
| 10.0 | 0.007901 | 0.301996 | 0.999999 |
score = 1 - delta.var()/temp.var() #计算拟合优度
print(score)
0.057448712094786636
如果对比之前建立的AR、MA模型，可以发现拟合精度上有所提升，但仍然是不够看的级别
2.4 ARIMA模型
ARIMA（AutoRegressive Integrated Moving Average）模型是一种广泛应用于时间序列分析和预测的经典模型。ARIMA模型结合了自回归（AR）模型、差分（I）模型和滑动平均（MA）模型的特性，用于描述时间序列数据中的自相关性、趋势性和随机波动。
ARIMA模型通常适用于具有一定自相关性和趋势性的时间序列数据。它的基本原理是在数据中引入差分操作，使得时间序列数据在一阶或多阶差分后变得平稳，然后再结合自回归和滑动平均部分来捕捉数据的自相关性和随机波动。
ARIMA(p, d, q)模型的表示为：
其中，
在ARIMA模型中， 表示需要进行的差分次数，用来使得时间序列数据变得平稳。当
时，ARIMA模型就退化为ARMA模型。
ARIMA模型的选择和参数估计通常需要通过观察时间序列数据的ACF（自相关函数）和PACF（偏自相关函数）来确定合适的阶数 、
和
。然后，可以使用最小二乘估计等方法来估计ARIMA模型的参数。
ARIMA模型在时间序列分析和预测中被广泛应用，可以用于对未来的数据进行预测和对过去的数据进行拟合，从而提供对时间序列数据的更深入理解和更准确的预测。
2.4.1 平稳性检验（ADF）
见1.5.2 使用收盘价数据
d2=data['close'] # 贵州茅台收盘价
d2.plot(figsize=(15,4))
from statsmodels.tsa.stattools import adfuller as ADF #平稳性检测
temp = np.array(d2)
t = ADF(temp) # ADF检验
output=pd.DataFrame(index=['Test Statistic Value', "p-value", "Lags Used", "Number of Observations Used","Critical Value(1%)","Critical Value(5%)","Critical Value(10%)"],columns=['value'])
output['value']['Test Statistic Value'] = t[0]
output['value']['p-value'] = t[1]
output['value']['Lags Used'] = t[2]
output['value']['Number of Observations Used'] = t[3]
output['value']['Critical Value(1%)'] = t[4]['1%']
output['value']['Critical Value(5%)'] = t[4]['5%']
output['value']['Critical Value(10%)'] = t[4]['10%']
output
| value | |
|---|---|
| Test Statistic Value | -2.30762 |
| p-value | 0.169522 |
| Lags Used | 0 |
| Number of Observations Used | 727 |
| Critical Value(1%) | -3.439377 |
| Critical Value(5%) | -2.865524 |
| Critical Value(10%) | -2.568891 |
可以看出，p-value为0.1704489，大于显著性水平0.05。原假设：序列具有单位根即非平稳。不能被拒绝。因此收盘价序列为非平稳的。
data2Diff = d2.diff() # 差分
data2Diff.plot(figsize=(15,4))
temp = np.array(data2Diff)[1:] # 差分后第一个值为NaN,舍去
t = ADF(temp) # ADF检验
print("p-value: ",t[1])
p-value: 6.887197234767805e-24
显示已经平稳了，故d的取值可为1阶
d=1
2.4.2 ARIMA(p,d,q)模型阶次确定
上面初步确定了d的取值，接着结合差分后序列的ACF和PACF来判断q、p
temp = np.array(data2Diff)[1:] # 差分后第一个值为NaN,舍去
fig = plt.figure(figsize=(12,8))
ax1=fig.add_subplot(211)
fig = ACF(temp,lags=30,ax=ax1)
ax2 = fig.add_subplot(212)
fig = PACF(temp,lags=30,ax=ax2)
可以看出，模型的阶次为(26,26)。建模计算了太大。我们再看看AIC准则
sm.tsa.arma_order_select_ic(temp,max_ar=6,max_ma=4,ic='aic')['aic_min_order'] # AIC
(3, 3)
2.4.3 ARIMA模型建立及预测
根据确定的阶次 建立ARIMA(26,1,26)以及ARIMA(3,1,3)
from statsmodels.tsa.arima.model import ARIMA
temp = np.array(d2) # 载入收盘价序列
model =ARIMA(temp,order=(3, 1, 3))
res = model.fit()
print(res.summary())
plt.rcParams['font.sans-serif'] = ['simhei'] #字体为黑体
plt.rcParams['axes.unicode_minus'] = False #正常显示负号
plt.figure(figsize=(10,4))
plt.plot(temp,'b',label='收盘价')
plt.plot(res.fittedvalues[1:], 'r',label='ARIMA model')
plt.legend()
SARIMAX Results
==============================================================================
Dep. Variable: y No. Observations: 728
Model: ARIMA(3, 1, 3) Log Likelihood -3638.880
Date: Wed, 09 Aug 2023 AIC 7291.759
Time: 22:53:07 BIC 7323.882
Sample: 0 HQIC 7304.155
- 728
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 -1.8971 0.085 -22.389 0.000 -2.063 -1.731
ar.L2 -1.8304 0.090 -20.333 0.000 -2.007 -1.654
ar.L3 -0.8498 0.074 -11.480 0.000 -0.995 -0.705
ma.L1 1.8984 0.095 19.919 0.000 1.712 2.085
ma.L2 1.8679 0.101 18.492 0.000 1.670 2.066
ma.L3 0.8482 0.091 9.288 0.000 0.669 1.027
sigma2 1300.6892 50.633 25.689 0.000 1201.451 1399.928
===================================================================================
Ljung-Box (L1) (Q): 0.03 Jarque-Bera (JB): 137.57
Prob(Q): 0.87 Prob(JB): 0.00
Heteroskedasticity (H): 1.48 Skew: -0.07
Prob(H) (two-sided): 0.00 Kurtosis: 5.13
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
from statsmodels.tsa.arima.model import ARIMA
temp = np.array(d2) # 载入收盘价序列
model =ARIMA(temp,order=(26, 1, 26))
res = model.fit()
print(res.summary())
plt.rcParams['font.sans-serif'] = ['simhei'] #字体为黑体
plt.rcParams['axes.unicode_minus'] = False #正常显示负号
plt.figure(figsize=(10,4))
plt.plot(temp,'b',label='收盘价')
plt.plot(res.fittedvalues[1:], 'r',label='ARIMA model')
plt.legend()
SARIMAX Results
==============================================================================
Dep. Variable: y No. Observations: 728
Model: ARIMA(26, 1, 26) Log Likelihood -3613.155
Date: Thu, 03 Aug 2023 AIC 7332.310
Time: 12:02:36 BIC 7575.523
Sample: 0 HQIC 7426.161
- 728
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 -0.1248 0.375 -0.332 0.740 -0.861 0.611
ar.L2 -0.0399 0.252 -0.158 0.874 -0.534 0.454
ar.L3 0.0831 0.254 0.327 0.744 -0.416 0.582
ar.L4 -0.2205 0.287 -0.768 0.443 -0.783 0.342
ar.L5 0.0155 0.308 0.050 0.960 -0.588 0.619
ar.L6 0.0376 0.281 0.134 0.894 -0.514 0.589
ar.L7 0.2008 0.199 1.007 0.314 -0.190 0.591
ar.L8 -0.0452 0.249 -0.182 0.856 -0.533 0.442
ar.L9 0.0180 0.249 0.072 0.942 -0.469 0.505
ar.L10 -0.1558 0.210 -0.741 0.459 -0.568 0.256
ar.L11 0.1428 0.256 0.558 0.577 -0.359 0.645
ar.L12 0.2120 0.261 0.813 0.416 -0.299 0.723
ar.L13 -0.1820 0.210 -0.868 0.385 -0.593 0.229
ar.L14 -0.1011 0.215 -0.471 0.637 -0.522 0.319
ar.L15 -0.0802 0.193 -0.415 0.678 -0.459 0.299
ar.L16 0.0364 0.201 0.181 0.857 -0.358 0.431
ar.L17 0.1169 0.160 0.731 0.465 -0.196 0.430
ar.L18 -0.0500 0.182 -0.274 0.784 -0.408 0.308
ar.L19 -0.1491 0.191 -0.779 0.436 -0.524 0.226
ar.L20 -0.1627 0.161 -1.008 0.313 -0.479 0.153
ar.L21 0.2864 0.195 1.467 0.142 -0.096 0.669
ar.L22 0.0443 0.232 0.191 0.848 -0.410 0.499
ar.L23 -0.0411 0.195 -0.211 0.833 -0.424 0.341
ar.L24 0.2800 0.189 1.480 0.139 -0.091 0.651
ar.L25 0.5272 0.176 2.994 0.003 0.182 0.872
ar.L26 0.1797 0.258 0.696 0.487 -0.327 0.686
ma.L1 0.1586 0.369 0.430 0.667 -0.564 0.882
ma.L2 0.0997 0.247 0.404 0.686 -0.384 0.583
ma.L3 -0.1194 0.249 -0.479 0.632 -0.608 0.369
ma.L4 0.2670 0.293 0.912 0.362 -0.307 0.841
ma.L5 -0.0450 0.322 -0.140 0.889 -0.675 0.585
ma.L6 -0.0618 0.293 -0.211 0.833 -0.637 0.513
ma.L7 -0.2217 0.205 -1.079 0.281 -0.624 0.181
ma.L8 0.0909 0.258 0.353 0.724 -0.414 0.596
ma.L9 -0.0613 0.262 -0.234 0.815 -0.575 0.453
ma.L10 0.2255 0.226 0.999 0.318 -0.217 0.668
ma.L11 -0.1986 0.285 -0.696 0.486 -0.758 0.361
ma.L12 -0.1590 0.296 -0.536 0.592 -0.740 0.422
ma.L13 0.0957 0.234 0.409 0.682 -0.363 0.554
ma.L14 0.1560 0.231 0.676 0.499 -0.296 0.608
ma.L15 0.0241 0.244 0.099 0.921 -0.454 0.502
ma.L16 0.0315 0.234 0.135 0.893 -0.427 0.490
ma.L17 -0.1301 0.186 -0.699 0.485 -0.495 0.235
ma.L18 0.0173 0.214 0.081 0.936 -0.403 0.437
ma.L19 0.1058 0.215 0.491 0.623 -0.316 0.528
ma.L20 0.2322 0.169 1.376 0.169 -0.099 0.563
ma.L21 -0.2887 0.231 -1.249 0.212 -0.742 0.164
ma.L22 -0.0226 0.280 -0.081 0.936 -0.572 0.527
ma.L23 0.0682 0.219 0.311 0.756 -0.362 0.498
ma.L24 -0.2839 0.205 -1.385 0.166 -0.686 0.118
ma.L25 -0.5246 0.176 -2.983 0.003 -0.869 -0.180
ma.L26 -0.3573 0.283 -1.261 0.207 -0.913 0.198
sigma2 1296.4820 65.956 19.657 0.000 1167.211 1425.753
===================================================================================
Ljung-Box (L1) (Q): 0.12 Jarque-Bera (JB): 118.88
Prob(Q): 0.73 Prob(JB): 0.00
Heteroskedasticity (H): 1.39 Skew: -0.12
Prob(H) (two-sided): 0.01 Kurtosis: 4.97
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
delta = res.fittedvalues - temp # 残差
score = 1 - delta[1:].var()/temp[1:].var() #计算拟合优度
print(score)
0.9870172288192488
# .predict(10) 进行预测
可以看出 建立的ARIMA（3,1,3）和ARIMA（26,1,26）模型对收盘价序列的拟合效果整体上是不错的，接着我们对其差分序列进行检验。
from statsmodels.tsa.arima.model import ARIMA
temp = np.array(data2Diff) # 载入收盘价一阶差分序列
model =ARIMA(temp,order=(3, 0, 3))
res = model.fit()
print(res.summary())
plt.rcParams['font.sans-serif'] = ['simhei'] #字体为黑体
plt.rcParams['axes.unicode_minus'] = False #正常显示负号
plt.figure(figsize=(10,4))
plt.plot(temp,'b',label='一阶差分close')
plt.plot(res.fittedvalues[1:], 'r',label='ARMA model2')
plt.legend()
SARIMAX Results
==============================================================================
Dep. Variable: y No. Observations: 728
Model: ARIMA(3, 0, 3) Log Likelihood -3638.660
Date: Thu, 03 Aug 2023 AIC 7293.319
Time: 12:02:38 BIC 7330.042
Sample: 0 HQIC 7307.489
- 728
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 0.9330 1.408 0.663 0.508 -1.826 3.692
ar.L1 -1.8995 0.082 -23.293 0.000 -2.059 -1.740
ar.L2 -1.8284 0.087 -20.998 0.000 -1.999 -1.658
ar.L3 -0.8499 0.070 -12.167 0.000 -0.987 -0.713
ma.L1 1.9009 0.092 20.625 0.000 1.720 2.082
ma.L2 1.8684 0.098 19.129 0.000 1.677 2.060
ma.L3 0.8498 0.087 9.745 0.000 0.679 1.021
sigma2 1300.5128 50.576 25.714 0.000 1201.386 1399.639
===================================================================================
Ljung-Box (L1) (Q): 0.02 Jarque-Bera (JB): 135.40
Prob(Q): 0.88 Prob(JB): 0.00
Heteroskedasticity (H): 1.49 Skew: -0.07
Prob(H) (two-sided): 0.00 Kurtosis: 5.11
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
from statsmodels.tsa.arima.model import ARIMA
temp = np.array(data2Diff) # 载入收盘价一阶差分序列
model =ARIMA(temp,order=(26, 0, 26))
res = model.fit()
print(res.summary())
plt.rcParams['font.sans-serif'] = ['simhei'] #字体为黑体
plt.rcParams['axes.unicode_minus'] = False #正常显示负号
plt.figure(figsize=(10,4))
plt.plot(temp,'b',label='一阶差分close')
plt.plot(res.fittedvalues[1:], 'r',label='ARMA model2')
plt.legend()
SARIMAX Results
==============================================================================
Dep. Variable: y No. Observations: 728
Model: ARIMA(26, 0, 26) Log Likelihood -3612.311
Date: Thu, 03 Aug 2023 AIC 7332.622
Time: 12:03:15 BIC 7580.498
Sample: 0 HQIC 7428.267
- 728
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 0.9282 1.416 0.655 0.512 -1.848 3.704
ar.L1 -0.1491 0.418 -0.357 0.721 -0.969 0.670
ar.L2 -0.0441 0.253 -0.174 0.862 -0.541 0.453
ar.L3 0.0976 0.245 0.398 0.691 -0.383 0.578
ar.L4 -0.2318 0.284 -0.816 0.415 -0.789 0.325
ar.L5 0.0079 0.331 0.024 0.981 -0.641 0.657
ar.L6 0.0305 0.284 0.107 0.914 -0.526 0.587
ar.L7 0.2151 0.199 1.082 0.279 -0.174 0.605
ar.L8 -0.0496 0.253 -0.196 0.845 -0.545 0.446
ar.L9 0.0146 0.252 0.058 0.954 -0.479 0.508
ar.L10 -0.1523 0.204 -0.748 0.454 -0.551 0.247
ar.L11 0.1296 0.260 0.499 0.618 -0.379 0.639
ar.L12 0.2153 0.273 0.788 0.431 -0.320 0.751
ar.L13 -0.1653 0.190 -0.868 0.385 -0.539 0.208
ar.L14 -0.1171 0.200 -0.586 0.558 -0.509 0.275
ar.L15 -0.0722 0.185 -0.391 0.696 -0.434 0.289
ar.L16 0.0344 0.188 0.183 0.855 -0.335 0.404
ar.L17 0.1232 0.154 0.801 0.423 -0.178 0.425
ar.L18 -0.0466 0.171 -0.272 0.786 -0.382 0.289
ar.L19 -0.1532 0.186 -0.822 0.411 -0.519 0.212
ar.L20 -0.1639 0.149 -1.100 0.271 -0.456 0.128
ar.L21 0.2835 0.181 1.565 0.117 -0.071 0.638
ar.L22 0.0529 0.242 0.219 0.827 -0.421 0.526
ar.L23 -0.0396 0.181 -0.219 0.826 -0.394 0.315
ar.L24 0.2481 0.178 1.394 0.163 -0.101 0.597
ar.L25 0.5444 0.161 3.390 0.001 0.230 0.859
ar.L26 0.2154 0.274 0.787 0.431 -0.321 0.752
ma.L1 0.1806 0.411 0.439 0.661 -0.626 0.987
ma.L2 0.1009 0.249 0.404 0.686 -0.388 0.590
ma.L3 -0.1344 0.243 -0.553 0.580 -0.610 0.342
ma.L4 0.2777 0.294 0.945 0.345 -0.298 0.853
ma.L5 -0.0458 0.350 -0.131 0.896 -0.733 0.641
ma.L6 -0.0523 0.305 -0.171 0.864 -0.651 0.546
ma.L7 -0.2385 0.211 -1.132 0.257 -0.651 0.174
ma.L8 0.1005 0.267 0.376 0.707 -0.423 0.624
ma.L9 -0.0619 0.278 -0.223 0.823 -0.606 0.482
ma.L10 0.2192 0.232 0.944 0.345 -0.236 0.674
ma.L11 -0.1922 0.301 -0.639 0.523 -0.782 0.398
ma.L12 -0.1593 0.320 -0.498 0.619 -0.786 0.468
ma.L13 0.0693 0.223 0.311 0.756 -0.367 0.506
ma.L14 0.1730 0.220 0.788 0.431 -0.257 0.603
ma.L15 0.0100 0.251 0.040 0.968 -0.483 0.503
ma.L16 0.0384 0.234 0.164 0.869 -0.420 0.497
ma.L17 -0.1403 0.194 -0.724 0.469 -0.520 0.239
ma.L18 0.0180 0.219 0.082 0.934 -0.411 0.447
ma.L19 0.1084 0.217 0.500 0.617 -0.316 0.533
ma.L20 0.2350 0.161 1.464 0.143 -0.080 0.550
ma.L21 -0.2915 0.231 -1.262 0.207 -0.744 0.161
ma.L22 -0.0245 0.304 -0.081 0.936 -0.620 0.571
ma.L23 0.0594 0.216 0.275 0.784 -0.364 0.483
ma.L24 -0.2556 0.200 -1.279 0.201 -0.647 0.136
ma.L25 -0.5529 0.166 -3.324 0.001 -0.879 -0.227
ma.L26 -0.3927 0.310 -1.268 0.205 -1.000 0.214
sigma2 1294.7855 67.427 19.203 0.000 1162.631 1426.940
===================================================================================
Ljung-Box (L1) (Q): 0.10 Jarque-Bera (JB): 115.33
Prob(Q): 0.75 Prob(JB): 0.00
Heteroskedasticity (H): 1.40 Skew: -0.12
Prob(H) (two-sided): 0.01 Kurtosis: 4.94
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
delta = res.fittedvalues - temp # 残差
score = 1 - delta[1:].var()/temp[1:].var()
print(score)
0.099128474863452
差分序列ARMA模型的拟合效果比较一般，若希望准确率更高，可能要建立更加高阶/高级的模型
小结:
本节主要介绍了AR、MA、ARMA、ARIMA四种最常见的时间序列模型进行了介绍，其中核心是ARIMA模型，其余都是ARIMA模型的简化版本，主要难点在于阶数的确定和建立高阶的模型。
案例2：
自选金融资产，对数据进行必要整理
（1）绘制ACF与PACF图形；
（2）检验数据的平稳性；
（3）构建ARIMA模型；
（4）对数据进行5阶外推预测
关注gzh“finance 褪黑素”获取金融、大数据等科目案例代码和更多学习资料。
版权声明：本文为博主作者：finance褪黑素原创文章，版权归属原作者，如果侵权，请联系我们删除！
原文链接：https://blog.csdn.net/celiaweiwei/article/details/133844753
