时间序列问题研究Task2
import warnings
warnings.simplefilter('ignore')1 时间序列建模方法汇总演示
1.1 ARIMA建模
ARMA建模的前提是输入的时间序列必须是平稳序列，因此要对输入数据进行平稳检验，如果不平稳则要进行平稳化处理。并且平稳化处理之后检查序列是否为白噪声序列，如果序列为白噪声序列则没有建模意义。
1.1.1 读取数据
首先读取数据，采用statsmodels中自带的美国夏威夷莫纳罗亚天文台大气二氧化碳数据，并转化为月度间隔数据
import numpy as np
import pandas as pd
from statsmodels.datasets import co2
from statsmodels.tsa.stattools import adfuller
from matplotlib import pyplot as plt
# 加载数据
data = co2.load(as_pandas=True).data
# 月度数据
data = data['co2']
data = data.resample('M').mean().ffill()
# 可视化
data.plot(figsize=(12,4))
plt.title('co2 data')
plt.show()
png
从上图可以看出二氧化碳数据有明显的趋势，不为平稳序列，为了相对客观的判别，使用ADF检验。
1.1.2 平稳检验
res = adfuller(data)
print('p value:', res[1])p value: 0.9989453312516823adf检验p值为0.99，明显不平稳，对原始序列进行处理，使其平稳。
1.1.3 差分平稳化
由于是二氧化碳的月度数据，其自然会有季节趋势在内，所以使用季节差分。
# # 一阶差分
# data_diff1 = data.diff()
# # 二阶差分
# data_diff2 = data_diff1.diff()
# 季节差分
data_diff = data.diff(12).dropna()
data_diff.plot(figsize=(12,4))
plt.title('co2 - seasonal difference')
plt.show()
# ADF检验
res = adfuller(data_diff)
print('p value:', res[1])
png
p value: 0.000778539367442881季节差分后，adf检验p值小于0.05，且时序图可以看出其基本平稳。
1.1.4 白噪声检验
检验差分后的序列是否为白噪声序列
from statsmodels.stats.diagnostic import acorr_ljungbox
res = acorr_ljungbox(data_diff, lags=[6,12,24], return_df=True)
print(res) lb_stat lb_pvalue
6 1249.517352 9.173589e-267
12 1534.064761 0.000000e+00
24 1547.525124 0.000000e+00使用lb检验，滞后6，12，24期p值均小于0.05，说明序列不为白噪声序列
1.1.5 模型构建
通过之前的序列平稳化操作，我们让序列平稳化并且不为白噪声，接下来通过对数据绘制其ACF和PACF图来进行人工模型定阶。
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
#fig, ax = plt.subplots(1, 2,figsize = (10,10))
#fig.subplots_adjust(hspace=1,wspace=0.2)
plot_acf(data_diff )
#ax[0][0].set_title('ACF(white_noise)')
plot_pacf(data_diff)
#ax[0][1].set_title('PACF(white_noise)')
plt.show()
png
png
绘制出差分后的模型的ACF和PACF图，通过如下的表格来人工确定模型阶数。
| 模型 | ACF | PACF |
|---|---|---|
| AR(p) | 拖尾 | p阶截尾 |
| MA(q) | q阶截尾 | 拖尾 |
| ARMA(p,q) | 拖尾 | 拖尾 |
从上面的ACF和PACF图可以手动定阶为ARMA(8,3),这里是人工定阶，难度较大，需要对模型原理有一定认识，并且会存在数据ACF图和PACF图不好定阶的情形，所以建议使用信息准则自动给模型定阶。
from statsmodels.tsa.arima_model import ARMA
model1 = ARMA(data_diff,order=(8,3)).fit(disp = 1)
print(model1.summary()) ARMA Model Results
==============================================================================
Dep. Variable: co2 No. Observations: 514
Model: ARMA(8, 3) Log Likelihood -234.507
Method: css-mle S.D. of innovations 0.379
Date: Mon, 09 May 2022 AIC 495.015
Time: 17:31:19 BIC 550.163
Sample: 03-31-1959 HQIC 516.629
- 12-31-2001
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 1.2968 0.103 12.621 0.000 1.095 1.498
ar.L1.co2 1.2418 0.052 23.657 0.000 1.139 1.345
ar.L2.co2 0.4291 0.071 6.026 0.000 0.290 0.569
ar.L3.co2 -1.3402 0.079 -17.051 0.000 -1.494 -1.186
ar.L4.co2 0.5042 0.094 5.354 0.000 0.320 0.689
ar.L5.co2 0.1648 0.094 1.757 0.079 -0.019 0.349
ar.L6.co2 -0.1359 0.077 -1.765 0.078 -0.287 0.015
ar.L7.co2 0.1521 0.073 2.097 0.036 0.010 0.294
ar.L8.co2 -0.1034 0.048 -2.136 0.033 -0.198 -0.009
ma.L1.co2 -0.6607 0.040 -16.435 0.000 -0.740 -0.582
ma.L2.co2 -0.7286 0.044 -16.395 0.000 -0.816 -0.641
ma.L3.co2 0.9322 0.034 27.545 0.000 0.866 0.998
Roots
=============================================================================
Real Imaginary Modulus Frequency
-----------------------------------------------------------------------------
AR.1 0.9041 -0.7182j 1.1547 -0.1068
AR.2 0.9041 +0.7182j 1.1547 0.1068
AR.3 1.1989 -0.1749j 1.2116 -0.0231
AR.4 1.1989 +0.1749j 1.2116 0.0231
AR.5 -1.1947 -0.1553j 1.2047 -0.4794
AR.6 -1.1947 +0.1553j 1.2047 0.4794
AR.7 -0.1731 -1.8370j 1.8451 -0.2650
AR.8 -0.1731 +1.8370j 1.8451 0.2650
MA.1 -1.0000 -0.0000j 1.0000 -0.5000
MA.2 0.8908 -0.5285j 1.0357 -0.0852
MA.3 0.8908 +0.5285j 1.0357 0.0852
-----------------------------------------------------------------------------上面输出了我们建立的ARMA(8,3)模型的细节，可以看出有部分参数的p值比0.05稍微大一点，这样的模型也可以使用，但是为了提升模型效果，我们进一步使用信息准则来自动化寻找模型参数。
1.1.6 使用信息准则自动化寻找模型参数
这里使用最小化BIC准则确定p,q。这个计算耗时略大为了控制计算量，这里限制AR阶数不超过10,MA最大阶数不超过4。
from statsmodels.tsa.stattools import arma_order_select_ic
#bic_min_order = arma_order_select_ic(data_diff,max_ar=8,max_ma=10,ic='bic')['bic_min_order']
#print(bic_min_order)这里将模型最大阶数分别设置到(10,10),在本地运行了22分钟才出结果。最终模型自动定阶为ARMA(8,10)
model2 = ARMA(data_diff,order=(8,10)).fit(disp=1)
print(model2.summary()) ARMA Model Results
==============================================================================
Dep. Variable: co2 No. Observations: 514
Model: ARMA(8, 10) Log Likelihood -173.335
Method: css-mle S.D. of innovations 0.331
Date: Mon, 09 May 2022 AIC 386.670
Time: 17:32:52 BIC 471.514
Sample: 03-31-1959 HQIC 419.923
- 12-31-2001
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 1.2928 0.110 11.715 0.000 1.077 1.509
ar.L1.co2 0.5929 0.084 7.090 0.000 0.429 0.757
ar.L2.co2 -0.5868 0.098 -5.968 0.000 -0.780 -0.394
ar.L3.co2 0.1013 0.098 1.031 0.303 -0.091 0.294
ar.L4.co2 0.1455 0.074 1.963 0.050 0.000 0.291
ar.L5.co2 -0.3106 0.058 -5.329 0.000 -0.425 -0.196
ar.L6.co2 0.4224 0.060 7.025 0.000 0.305 0.540
ar.L7.co2 -0.2488 0.059 -4.228 0.000 -0.364 -0.133
ar.L8.co2 0.2189 0.054 4.090 0.000 0.114 0.324
ma.L1.co2 -0.0515 0.077 -0.670 0.503 -0.202 0.099
ma.L2.co2 0.8900 0.060 14.815 0.000 0.772 1.008
ma.L3.co2 0.3326 0.099 3.352 0.001 0.138 0.527
ma.L4.co2 0.2951 0.055 5.341 0.000 0.187 0.403
ma.L5.co2 0.8529 0.052 16.461 0.000 0.751 0.954
ma.L6.co2 -0.1583 0.054 -2.959 0.003 -0.263 -0.053
ma.L7.co2 0.9139 0.050 18.110 0.000 0.815 1.013
ma.L8.co2 -0.0148 0.092 -0.161 0.872 -0.194 0.165
ma.L9.co2 0.5002 0.057 8.771 0.000 0.388 0.612
ma.L10.co2 0.5369 0.073 7.380 0.000 0.394 0.680
Roots
==============================================================================
Real Imaginary Modulus Frequency
------------------------------------------------------------------------------
AR.1 -1.1070 -0.0000j 1.1070 -0.5000
AR.2 -0.4916 -1.1114j 1.2153 -0.3163
AR.3 -0.4916 +1.1114j 1.2153 0.3163
AR.4 1.2232 -0.0000j 1.2232 -0.0000
AR.5 0.8046 -0.9110j 1.2154 -0.1349
AR.6 0.8046 +0.9110j 1.2154 0.1349
AR.7 0.1973 -1.2278j 1.2436 -0.2246
AR.8 0.1973 +1.2278j 1.2436 0.2246
MA.1 0.8962 -0.5241j 1.0382 -0.0842
MA.2 0.8962 +0.5241j 1.0382 0.0842
MA.3 0.4976 -0.8737j 1.0054 -0.1676
MA.4 0.4976 +0.8737j 1.0054 0.1676
MA.5 -0.0106 -1.0000j 1.0001 -0.2517
MA.6 -0.0106 +1.0000j 1.0001 0.2517
MA.7 -0.4948 -0.8692j 1.0001 -0.3324
MA.8 -0.4948 +0.8692j 1.0001 0.3324
MA.9 -1.0001 -0.0000j 1.0001 -0.5000
MA.10 -1.7083 -0.0000j 1.7083 -0.5000
------------------------------------------------------------------------------使用最小化BIC准则建立模型，模型的部分参数的p值比较大，但是整体的BIC最小，所以这两个是比较难以平衡的，对建立的模型进行进一步可视化探索。
可视化残差直方图、QQ图，看是否正态分布，可视化ACF看是否仍存有自相关性。
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf
fig, axs = plt.subplots(2, 2,figsize=(8,8))
fig.subplots_adjust(hspace=0.3)
model2.resid.plot(ax=axs[0][0])
axs[0][0].set_title('residual')
model2.resid.plot(kind='hist', ax=axs[0][1])
axs[0][1].set_title('histogram')
sm.qqplot(model2.resid, line='45', fit=True, ax=axs[1][0])
axs[1][0].set_title('Normal Q-Q')
plot_acf(model2.resid, ax=axs[1][1])
plt.show()
png
上图可以看出残差基本服从正态分布，残差的ACF图也基本都趋于0,所以最终将模型定为ARMA(8,10)
1.1.7 模型预测
使用建立好的ARMA(8,10)对数据进行预测。预测六期数据。
preds = model2.predict(0,len(data_diff)+6)
plt.figure(figsize=(12,4))
data_diff.plot(color='g', label='data_diff')
preds.plot(color='r', label='predict')
plt.legend()
plt.show()
png
1.1.8 差分还原
我们的建模是对原始数据进行了季节差分之后的建模，所以在预测之后==千万== 要记得将数据还原。
df1 = pd.DataFrame(data)
df2 = pd.DataFrame(preds, columns=['predict'])
df = pd.concat([df1, df2], axis=1)
df['result'] = df['predict'] + df['co2'].shift(12)
plt.figure(figsize=(12, 4))
df['result'].plot(color='r', label='arma result')
df['co2'].plot(color='g', label='co2 data')
plt.legend()
plt.show()
png
还原之后绘制出图片看出模型预测效果不错。但是需要量化评判指标。时间序列模型可以类似回归模型使用均方误差进行评价。
print("模型均方误差为：",(model2.resid**2).sum()/len(model2.resid))模型均方误差为： 0.117768180699071771.1.9 ARMA建模总结
ARMA是传统的时序建模，其中AR模型和MA模型基本都涵盖其中，是传统统计时序建模的基础。但是也存在模型建立原理较为复杂，模型对数据要求较高：平稳，非白噪声。如果不属于的话还要进行平稳化操作。具体建模流程如下图。
ARMA模型的建立流程如上图所示。目前感觉可以通过程序来实现自动化建模，但是有可能输出的模型不是最优的模型，也可以让用户自己进行手动建模，但是应该提醒用户每一步的意义和操作。
1.2 GARCH建模
ARCH模型（Autoregressive conditional heteroskedasticity model）全称“自回归条件异方差模型”，解决了传统的计量经济学对时间序列变量的第二个假设（方差恒定）所引起的问题。GARCH模型称为广义ARCH模型，是ARCH模型的拓展，由Bollerslev(1986)发展起来的。
简单来说，GARCH模型就是ARMA模型应用在时序的方差上，它包含一个自回归项和一个移动平均项. 如果时间序列数据可以表示为：\[x_t = \sigma_t\omega_t\] 其中\(\{\omega_t\}\)是白噪声序列，均值为0，方差为1，这里\(\sigma_t\)为:\[\sigma_{t}^{2}=\alpha_{0}+\sum_{i=1}^{q} \alpha_{i} y_{t-i}^{2}+\sum_{j=1}^{p} \beta_{j} \sigma_{t-j}^{2}\]
一个简单的GARCH(1,1)模型： \[ \begin{gathered} y_{t}=\sigma_{t} w_{t} \\ \sigma_{t}^{2}=\alpha_{0}+\alpha_{1} y_{t-1}^{2}+\beta_{1} \sigma_{t-1}^{2} \end{gathered} \]
下面是建立ARCH模型的一般步骤:
（1）检验收益率序列是否平稳，根据自相关性建立合适的均值方程，如ARMA模型，描述收益率如何随时间变化，根据拟合的模型和实际值，得到残差序列。
（2）对拟合的均值方程得到的残差序列进行ARCH效应检验，即检验收益率围绕均值的偏差是否时大时小。检验序列是否具有ARCH效应的方法有两种：Ljung-Box检验和LM检验。
（3）若ARCH效应在统计上显著，则需要再设定一个波动率模型来刻画波动率的动态变化。
（4）对均值方差和波动率方差进行联合估计，即假设实际数据服从前面设定的均值方差和波动率方差后，对均值方差和波动率方差中的参数进行估计，并得到估计的误差。
（5）对拟合的模型进行检验。如果估计结果（残差项）不满足模型本身的假设，则模型的可用性较差。
简单来说，一般的金融数据由于现实的复杂性，很难满足平稳条件的方差不变性，当序列数据通过ARMA等模型建立之后的残差的平方具有自相关性，说明高阶的时间序列特征还没有捕捉到，这时需要GARCH模型建模。GARCH模型的核心思想是残差项的条件方差依赖于它的前期的值的大小。
1.2.1 ARCH效应检验
一般来说，要进行GARCH模型建立首先需要检验是否存在ARCH效应，具体可以用已经建立好的ARMA模型，与原始序列拟合的残差，然后对残差序列$ {_t ^ 2}\(做ARCH效应检验，这里介绍相对简单Ljung-Box白噪声检验。如果对\) {_t ^ 2}$进行白噪声检验，检验不显时候没有ARCH效应，检验显著时有ARCH效应。
1.2.2 ARCH模型建立
这里模拟一个GARCH(1,1)的数据
import statsmodels.tsa.api as smt
import pandas as pd
import numpy as np
import statsmodels.tsa.api as smt
import scipy.stats as scs #tsa为Time Series analysis缩写import statsmodels.api as smimport scipy.
%matplotlib inline
#正常显示画图时出现的中文和负号from pylab import mpl
import statsmodels.api as sm
import scipy.stats as scs
from arch import arch_model#画图
import matplotlib.pyplot as plt
import matplotlib as mpl
def tsplot(y, lags=None, figsize=(16, 10), style='bmh'):
if not isinstance(y, pd.Series):
y = pd.Series(y)
with plt.style.context(style):
fig = plt.figure(figsize=figsize)
# mpl.rcParams['font.family'] = 'Ubuntu Mono'
layout = (3, 2)
ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)
acf_ax = plt.subplot2grid(layout, (1, 0))
pacf_ax = plt.subplot2grid(layout, (1, 1))
qq_ax = plt.subplot2grid(layout, (2, 0))
pp_ax = plt.subplot2grid(layout, (2, 1))
y.plot(ax=ts_ax)
ts_ax.set_title('Time Series Analysis Plots')
smt.graphics.plot_acf(y, lags=lags, ax=acf_ax, alpha=0.5)
smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax, alpha=0.5)
sm.qqplot(y, line='s', ax=qq_ax)
qq_ax.set_title('QQ Plot')
scs.probplot(y, sparams=(y.mean(), y.std()), plot=pp_ax)
plt.tight_layout()
return# Simulating a GARCH(1, 1) process
np.random.seed(2)
a0 = 0.2
a1 = 0.5
b1 = 0.3
n = 10000
w = np.random.normal(size=n)
y = np.zeros_like(w)
sigsq = np.zeros_like(w)
for i in range(1, n):
sigsq[i] = a0 + a1*(y[i-1]**2) + b1*sigsq[i-1]
y[i] = w[i] * np.sqrt(sigsq[i])
_ = tsplot(y, lags=30)
png
模拟出的时间序列很类似于一个白噪声序列，但是将其去平方之后再次观察。
tsplot(y*y, lags=30)
png
从上图可以看出在进行平方之后的序列出现了波动聚集效应，对平方后的序列进行白噪声检验。
from statsmodels.stats.diagnostic import acorr_ljungbox as lb_test
print('LB检验P值:',lb_test(y**2)[1])LB检验P值: [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]由输出结果可以看到，，p值都很小，所以该序列拒绝原假设，即认为该数据不是是纯随机数据（白噪声）。 接着拟合一个GARCH(1,1)模型
from arch import arch_model
am = arch_model(y,p=1,q=1)
res = am.fit(update_freq=5)
print(res.summary())Iteration: 5, Func. Count: 34, Neg. LLF: 13129.789119756882
Iteration: 10, Func. Count: 62, Neg. LLF: 12237.312782262137
Optimization terminated successfully (Exit mode 0)
Current function value: 12237.303267318166
Iterations: 13
Function evaluations: 76
Gradient evaluations: 13
Constant Mean - GARCH Model Results
==============================================================================
Dep. Variable: y R-squared: 0.000
Mean Model: Constant Mean Adj. R-squared: 0.000
Vol Model: GARCH Log-Likelihood: -12237.3
Distribution: Normal AIC: 24482.6
Method: Maximum Likelihood BIC: 24511.4
No. Observations: 10000
Date: Mon, May 09 2022 Df Residuals: 9999
Time: 17:33:05 Df Model: 1
Mean Model
==============================================================================
coef std err t P>|t| 95.0% Conf. Int.
------------------------------------------------------------------------------
mu -6.7225e-03 6.735e-03 -0.998 0.318 [-1.992e-02,6.478e-03]
Volatility Model
========================================================================
coef std err t P>|t| 95.0% Conf. Int.
------------------------------------------------------------------------
omega 0.2021 1.043e-02 19.383 1.084e-83 [ 0.182, 0.223]
alpha[1] 0.5162 2.016e-02 25.611 1.144e-144 [ 0.477, 0.556]
beta[1] 0.2879 1.870e-02 15.395 1.781e-53 [ 0.251, 0.325]
========================================================================
Covariance estimator: robust然后再看GARCH模型的残差的相关图片。
_ = tsplot(res.resid,lags=30)
png
可以看出构建GARCH(1,1)模型之后的残差近似于白噪声序列，进一步进行LB检验。
lb_test(res.resid,lags=[6,12],return_df=True)| lb_stat | lb_pvalue | |
|---|---|---|
| 6 | 17.352678 | 0.008071 |
| 12 | 21.605643 | 0.042185 |
从LB检验的P值看出残差序列的LB检验P值小于0.05，可以认为是白噪声序列。