时间序列预测之ARIMA
python 3.8
statsmodels 0.13.5
https://www.biaodianfu.com/arima.html
时间序列与平稳性
在数学上，随机过程被定义为一族时间随机变量，即{x(t),t∈T},其中T表示时间t的变动范围。当T={0,±1,±2,…}时，此类随机过程x(t)是离散时间t的随机函数，称为时间序列。时间序列的构成要素有：
长期趋势(T long-term trend)是在较长时期内受某种根本性因素作用而形成的总的变动趋势
季节变动(S Seasonal variation)是在一年内随着季节的变化而发生的有规律的周期性变动。它是诸如气候条件、生产条件、节假日或人们的风俗习惯等各种因素影响的结果。
循环变动(C Cyclical change)是时间序列呈现出得非固定长度的周期性变动。循环波动的周期可能会持续一段时间，但与趋势不同，它不是朝着单一方向的持续变动，而是涨落相同的交替波动。
不规则变动(I Irregular changes)是时间序列中除去趋势、季节变动和周期波动之后的随机波动。不规则波动通常总是夹杂在时间序列中，致使时间序列产生一种波浪形或震荡式的变动。只含有随机波动的序列也称为平稳序列。
除非你的时间序列是平稳的，否则不能建立一个时间序列模型。在百度词条中是这样粗略的讲的：平稳时间序列粗略地讲，一个时间序列，如果均值没有系统的变化（无趋势）、方差没有系统变化，且严格消除了周期性变化，就称之是平稳的。
平稳性的判断标准
判断一个序列是不是平稳序列有三个评判标准：
均值
均值，是与时间t无关的常数。下图（左）满足平稳序列的条件，下图（右）很明显具有时间依赖。
方差
方差，是与时间t无关的常数。这个特性叫做方差齐性。下图显示了什么是方差对齐，什么不是方差对齐。（注意右手边图中的不同分布）
协方差
协方差，只与时期间隔k有关，与时间t无关的常数。如下图（右），可以注意到随着时间的增加，曲线变得越来越近。因此红色序列的协方差并不是恒定的。
平稳性分类
时间序列的平稳性，和其它随机过程一样，分为严平稳和宽平稳。
严平稳
如果对所有的时刻t，任意正整数k和任意k个正整数 ( t 1 , t 2 , … , t k ) (t_1,t_2,…,t_k) (t1,t2,…,tk), ( r t 1 , r t 2 , … , r t k ) (r_{t1},r_{t2},…,r_{tk}) (rt1,rt2,…,rtk)的联合分布与 ( r t 1 + t , r t 2 + t , … , r t k + t ) (r_{t1+t},r_{t2+t},…,r_{tk+t}) (rt1+t,rt2+t,…,rtk+t)的联合分布相同，我们称时间序列 r t r_t rt是严平稳的。也就是， ( r t 1 , r t 2 , … , r t k ) (r_{t1},r_{t2},…,r_{tk}) (rt1,rt2,…,rtk)的联合分布在时间的平移变换下保持不变，这是个很强的条件。
宽平稳
若时间序列{ r t r_t rt}满足下面两个条件： E ( r t ) = μ E(r_t)=\mu E(rt)=μ， μ \mu μ是常数， C o v ( r t , r t − l ) = γ l Cov(r_t,r_{t−l})=\gamma_l Cov(rt,rt−l)=γl， γ l \gamma_l γl只依赖于 l l l则时间序列{ r t r_t rt}为弱平稳的。即该序列的均值， r t r_t rt与 r t − l r_{t−l} rt−l的协方差不随时间而改变， l l l为任意整数。
平稳性检测方法
在数学中，平稳随机过程（Stationary random process）或者严平稳随机过程（Strictly-sense stationary random process），又称狭义平稳过程，是在固定时间和位置的概率分布与所有时间和位置的概率分布相同的随机过程：即随机过程的统计特性不随时间的推移而变化。这样，数学期望和方差这些参数也不随时间和位置变化。平稳在理论上有严平稳和宽平稳两种，在实际应用上宽平稳使用较多。宽平稳的数学定义为：对于时间序列
y
t
y_t
yt，若对任意的
t
t
t,
k
k
k,
m
m
m，满足：
E
(
y
t
)
=
E
(
y
t
+
m
)
E(y_t)=E(y_{t+m})
E(yt)=E(yt+m)
c
o
v
(
y
t
,
y
t
+
k
)
=
c
o
v
(
y
t
+
k
,
y
t
+
k
+
m
)
cov(y_t,y_{t+k})=cov(y_{t+k},y_{t+k+m})
cov(yt,yt+k)=cov(yt+k,yt+k+m)
则称时间序列
y
t
y_t
yt是宽平稳的。
单位根检验
单位根检验是指检验序列中是否存在单位根，因为存在单位根就是非平稳时间序列了。Dickey-Fuller test和Augmented Dickey Fuller test（ADF）可以测试一个自回归模型是否存在单位根（unit root）。
为了保持单位和平均数相似，这里采用了标准差来代替方差：
import matplotlib.pylab as plt
import numpy as np
import pandas as pd
data = pd.read_csv(
"AirPassengers.csv", parse_dates=["TravelDate"], index_col="TravelDate"
)
print(data.head())
Passengers
TravelDate
1949-01-01 112
1949-02-01 118
1949-03-01 132
1949-04-01 129
1949-05-01 121
from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries, window=12):
rolmean = timeseries.rolling(window=window, center=False).mean()
rolstd = timeseries.rolling(window=window, center=False).std()
# 设置原始图，移动平均图和标准差图的式样
orig = plt.plot(timeseries, color="blue", label="Original")
mean = plt.plot(rolmean, color="red", label="Rolling Mean")
std = plt.plot(rolstd, color="black", label="Rolling Std")
# 使用自动最佳的图例显示位置
plt.legend(loc="best")
plt.title("Rolling Mean & Standard Deviation")
# 供肉眼观察是否平稳
plt.show()
print("ADF检验结果：")
# 使用减小AIC的办法估算ADF测试所需的滞后数
dftest = adfuller(timeseries, autolag="AIC")
print(dftest)
# 将ADF测试结果、显著性概率、所用的滞后数和所用的观测数打印出来
dfoutput = pd.Series(
dftest[0:4],
index=["Test Statistic", "p-value", "Num Lags Used", "Num Observations Used"],
)
for key, value in dftest[4].items():
dfoutput["Critical Value (%s)" % key] = value
print(dfoutput)
ts_week = data["Passengers"]
test_stationarity(ts_week)
ADF检验结果：
(0.8153688792060482, 0.991880243437641, 13, 130, {'1%': -3.4816817173418295, '5%': -2.8840418343195267, '10%': -2.578770059171598}, 996.692930839019)
Test Statistic 0.815369
p-value 0.991880
Num Lags Used 13.000000
Num Observations Used 130.000000
Critical Value (1%) -3.481682
Critical Value (5%) -2.884042
Critical Value (10%) -2.578770
dtype: float64
ADF检验测试结果由测试统计量和一些置信区间的临界值组成。如果“测试统计量”少于“临界值”，我们可以拒绝无效假设，并认为序列是稳定的。
Test Statistic的值如果比Critical Value (5%)小则满足稳定性需求
p-value越低（理论上需要低于0.5）证明序列越稳定。
白噪声检测
对于纯随机序列，又称白噪声序列，序列的各项数值之间没有任何相关关系，序列在进行完全无序的随机波动，可以终止对该序列的分析。白噪声序列是没有信息可提取的平稳序列。
# 生成一个伪随机白噪声用于测试acorr_ljungbox的可靠性
# 导入白噪声检验函数
from random import gauss
from statsmodels.stats.diagnostic import acorr_ljungbox
# 创建一个高斯分布的白噪声
whitenoise = pd.Series([gauss(0.0, 1.0) for i in range(1000)])
# 检验结果：平稳度，p-value。p-value>0.05为白噪声
print("白噪声检验结果为：\n", acorr_ljungbox(whitenoise, lags=1))
白噪声检验结果为：
lb_stat lb_pvalue
1 0.39984 0.527172
时间序列转为平稳序列
将时间序列转为平稳序列，常见的方法如：
Deflation by CPI
Logarithmic（取对数）
First Difference（一阶差分）
Seasonal Difference（季节差分）
Seasonal Adjustment
总结起来，影响时间序列不稳定的主要有2个因素：
趋势 Trend
季节性 Seasonality
测试数据：AirPassengers。这个数据集是1949-1960年每个月国际航空的乘客数量的数据。
加载数据
import matplotlib.pylab as plt
import numpy as np
import pandas as pd
data = pd.read_csv(
"AirPassengers.csv", parse_dates=["TravelDate"], index_col="TravelDate"
)
print(data.head())
ts = data["Passengers"]
ts.head(10)
plt.plot(ts)
Passengers
TravelDate
1949-01-01 112
1949-02-01 118
1949-03-01 132
1949-04-01 129
1949-05-01 121
[<matplotlib.lines.Line2D at 0x186d3f9e6d0>]
检测稳定性
test_stationarity(ts)
ADF检验结果：
(0.8153688792060482, 0.991880243437641, 13, 130, {'1%': -3.4816817173418295, '5%': -2.8840418343195267, '10%': -2.578770059171598}, 996.692930839019)
Test Statistic 0.815369
p-value 0.991880
Num Lags Used 13.000000
Num Observations Used 130.000000
Critical Value (1%) -3.481682
Critical Value (5%) -2.884042
Critical Value (10%) -2.578770
dtype: float64
这里这个结果表明这些序列很不稳定，所以接下来考虑如何处理数据，使得序列相对稳定。
消除趋势
消除趋势的第一个方法是转换。在本例中我们可以清楚地看到有一个显著的趋势。我们可以通过变换，惩罚较高值而不是较小值。这可以采用取对数log,平方根,立方根等。
在《Analysis of Financial Time Series》提到，对原始数据取对数有两个用处：
将指数增长转为线性增长
可以平稳序列的方差
ts_log = np.log(ts)
plt.plot(ts_log)
[<matplotlib.lines.Line2D at 0x186d40a66a0>]
这里我们可以明显看到一个上升的趋势，但是混杂在噪音当中，所以需要去除杂音。这里简单平滑一下数据。平滑的窗口取值12（一年有12个月）。
在这个例子中，很容易看到数据趋势，但是它表现的不是很直观。我们可以使用一些技术来估计或对这个趋势建模，然后将它从序列中删除。常用的方法有：
聚合-取一段时间的平均值（月/周平均值）
平滑-取滚动平均数
多项式回归分析-适合的回归模型
这里主要讨论平滑取滚动平均数方法。平滑是指采取滚动估计，即考虑过去的几个实例：
移动平均数
在这个方法中，根据时间序列的频率采用“K”连续值的平均数。我们可以采用过去一年的平均数，即过去12个月的平均数。
moving_avg = ts_log.rolling(12).mean()
plt.plot(ts_log)
plt.plot(moving_avg, color="red")
[<matplotlib.lines.Line2D at 0x186d4039310>]
红色表示了滚动平均数。让我们从原始序列中减去这个平均数。注意，从我们采用过去12个月的值开始，滚动平均法还没有对前11个月的值定义：
ts_log_moving_avg_diff = ts_log - moving_avg
ts_log_moving_avg_diff.head(20)
TravelDate
1949-01-01 NaN
1949-02-01 NaN
1949-03-01 NaN
1949-04-01 NaN
1949-05-01 NaN
1949-06-01 NaN
1949-07-01 NaN
1949-08-01 NaN
1949-09-01 NaN
1949-10-01 NaN
1949-11-01 NaN
1949-12-01 -0.065494
1950-01-01 -0.093449
1950-02-01 -0.007566
1950-03-01 0.099416
1950-04-01 0.052142
1950-05-01 -0.027529
1950-06-01 0.139881
1950-07-01 0.260184
1950-08-01 0.248635
Name: Passengers, dtype: float64
注意前11个月是NaN值，现在让我们对这11个月排除后测试稳定性。
ts_log_moving_avg_diff.dropna(inplace=True)
test_stationarity(ts_log_moving_avg_diff)
ADF检验结果：
(-3.162907991300869, 0.022234630001243157, 13, 119, {'1%': -3.4865346059036564, '5%': -2.8861509858476264, '10%': -2.579896092790057}, -436.63269481747034)
Test Statistic -3.162908
p-value 0.022235
Num Lags Used 13.000000
Num Observations Used 119.000000
Critical Value (1%) -3.486535
Critical Value (5%) -2.886151
Critical Value (10%) -2.579896
dtype: float64
从上面的测试结果看，已经得到了一个稳定的序列。但是，这个方法有一个缺陷：需要先设定平滑的窗口期，实际在应用的时候像股票这样的数据很难去设定窗口期，所以，我们采取“加权移动平均法”可以对最近的值赋予更高的权重。
加权移动平均
指数加权移动平均法是很受欢迎的方法，所有的权重被指定给先前的值连同衰减系数。这可以通过Pandas实现：
expwighted_avg = ts_log.ewm(halflife=12).mean()
plt.plot(ts_log)
plt.plot(expwighted_avg, color="red")
[<matplotlib.lines.Line2D at 0x186d62b29d0>]
注意,这里使用了参数“halflife”来定义指数衰减量。这只是一个假设，很大程度上取决于业务领域。其他参数，如跨度和质心也可以用来定义衰减。让我们再来检测下新得到的序列的稳定性：
ts_log_ewma_diff = ts_log - expwighted_avg
test_stationarity(ts_log_ewma_diff)
ADF检验结果：
(-3.601262420161663, 0.0057369388015121505, 13, 130, {'1%': -3.4816817173418295, '5%': -2.8840418343195267, '10%': -2.578770059171598}, -462.8241632483199)
Test Statistic -3.601262
p-value 0.005737
Num Lags Used 13.000000
Num Observations Used 130.000000
Critical Value (1%) -3.481682
Critical Value (5%) -2.884042
Critical Value (10%) -2.578770
dtype: float64
检测结果比移动平均的效果更好（Test Statistic的值比1%的临界值还小），另外不会出现前面11个月数据遗漏问题。
消除季节性
两种消除趋势和季节性的方法：
差分：采用一个特定时间差的差值
分解：建立有关趋势和季节性的模型和从模型中删除它们
差分
处理趋势和季节性的最常见的方法之一就是差分法。可以使用Pandas来实现差分：
ts_log_diff = ts_log.diff() # ts_log_diff = ts_log - ts_log.shift()
plt.plot(ts_log_diff)
[<matplotlib.lines.Line2D at 0x186d7428fd0>]
图中可以看出很大程度上减少了趋势。让我们再来检测下：
ts_log_diff.dropna(inplace=True)
test_stationarity(ts_log_diff)
ADF检验结果：
(-2.717130598388133, 0.07112054815085875, 14, 128, {'1%': -3.4825006939887997, '5%': -2.884397984161377, '10%': -2.578960197753906}, -440.35846985568105)
Test Statistic -2.717131
p-value 0.071121
Num Lags Used 14.000000
Num Observations Used 128.000000
Critical Value (1%) -3.482501
Critical Value (5%) -2.884398
Critical Value (10%) -2.578960
dtype: float64
可以看到平均数和标准差随着时间有小的变化。同时，DF检验统计量小于10% 的临界值，因此该时间序列在90%的置信区间上是稳定的。我们也可以采取二阶或三阶差分在具体应用中获得更好的结果。
分解
statsmodels.tsa.seasonal.seasonal_decompose可以将时间序列的数据分解为趋势(trend),季节性(seasonality)和残差(residual)三部分：
from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(ts_log)
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
plt.subplot(411)
plt.plot(ts_log, label="Original")
plt.legend(loc="best")
plt.subplot(412)
plt.plot(trend, label="Trend")
plt.legend(loc="best")
plt.subplot(413)
plt.plot(seasonal, label="Seasonality")
plt.legend(loc="best")
plt.subplot(414)
plt.plot(residual, label="Residuals")
plt.legend(loc="best")
plt.tight_layout()
趋势和季节性，还有残差值都被分解出来，然后我们计算残差值的稳定性。
ts_log_decompose = residual
ts_log_decompose.dropna(inplace=True)
test_stationarity(ts_log_decompose)
ADF检验结果：
(-6.33238679764488, 2.8850589033231607e-08, 9, 122, {'1%': -3.4851223522012855, '5%': -2.88553750045158, '10%': -2.5795685622144586}, -517.2642969403794)
Test Statistic -6.332387e+00
p-value 2.885059e-08
Num Lags Used 9.000000e+00
Num Observations Used 1.220000e+02
Critical Value (1%) -3.485122e+00
Critical Value (5%) -2.885538e+00
Critical Value (10%) -2.579569e+00
dtype: float64
DF测试统计量明显低于1%的临界值，这样时间序列是非常接近稳定。
时间序列建模
相关系数与自相关函数
相关系数
对于两个向量，我们希望定义它们是不是相关。一个很自然的想法，用向量与向量的夹角作为距离的定义，夹角小，就距离小，夹角大，就距离大。我们通常使用余弦公式来计算角度：
c
o
s
<
a
⃗
,
b
⃗
>
=
a
⃗
⋅
b
⃗
∣
a
⃗
∣
∣
b
⃗
∣
cos<\vec{a},\vec{b}>=\frac{\vec{a} \cdot \vec{b}}{|\vec{a}||\vec{b}|}
cos<a,b>=∣a∣∣b∣a⋅b
对于 a ⃗ ⋅ b ⃗ \vec{a} \cdot \vec{b} a⋅b我们叫做内积，例如：
(
x
1
,
y
1
)
⋅
(
x
2
,
y
2
)
=
x
1
x
2
+
y
1
y
2
(x_1,y_1)\cdot(x_2,y_2)=x_1x_2+y_1y_2
(x1,y1)⋅(x2,y2)=x1x2+y1y2
相关系数的定义公式，X和Y的相关系数为：
ρ
x
y
=
C
o
v
(
X
,
Y
)
V
a
r
(
X
)
V
a
r
(
Y
)
\rho_{xy}=\frac{Cov(X,Y)}{\sqrt{Var(X)Var(Y)}}
ρxy=Var(X)Var(Y)Cov(X,Y)
根据样本的估计计算公式为：
ρ x y = ∑ t = 1 T ( x t − x ‾ ) ( y t − y ‾ ) ∑ t = 1 T ( x t − x ‾ ) 2 ∑ t = 1 T ( y t − y ‾ ) 2 \rho_{xy}=\frac{\sum_{t=1}^T{(x_t-\overline{x})(y_t-\overline{y})}} {\sqrt{\sum_{t=1}^T{(x_t-\overline{x})^2}\sum_{t=1}^T{(y_t-\overline{y})^2}}} ρxy=∑t=1T(xt−x)2∑t=1T(yt−y)2∑t=1T(xt−x)(yt−y)
= ( X − x ‾ ) ⟶ ⋅ ( Y − y ‾ ) ⟶ ∣ ( X − x ‾ ) ⟶ ∣ ∣ ( Y − y ‾ ) ⟶ ∣ =\frac{\mathop{(X-\overline{x})} \limits^\longrightarrow \cdot \mathop{(Y-\overline{y})} \limits^\longrightarrow}{|\mathop{(X-\overline{x})} \limits^\longrightarrow||\mathop{(Y-\overline{y})} \limits^\longrightarrow|} =∣(X−x)⟶∣∣(Y−y)⟶∣(X−x)⟶⋅(Y−y)⟶
相关系数实际上就是计算了向量空间中两个向量的夹角，协方差是去均值后两个向量的内积。如果两个向量平行，相关系数等于1或者-1，同向的时候是1，反向的时候就是-1。如果两个向量垂直，则夹角的余弦就等于0，说明二者不相关。两个向量夹角越小，相关系数绝对值越接近1，相关性越高。只不过这里计算的时候对向量做了去均值处理，即中心化操作。而不是直接用向量X,Y计算。
自相关函数(Autocorrelation Function, ACF)
相关系数度量了两个向量的线性相关性，而在平稳时间序列 r t {r_t} rt中，我们有时候很想知道， r t r_t rt与它的过去值 r t − i r_{t−i} rt−i的线性相关性。这时候我们把相关系数的概念推广到自相关系数。
r t r_t rt与 r t − l r_{t−l} rt−l的相关系数称为 r t r_t rt的间隔为 l l l的自相关系数，通常记为 ρ l \rho_l ρl：
ρ l = C o v ( r t , r t − l ) V a r ( r t ) V a r ( r t − l ) = C o v ( r t , t t − l ) V a r ( r t ) \rho_l=\frac{Cov(r_t,r_{t-l})}{\sqrt{Var(r_t)Var(r_{t-l})}} =\frac{Cov(r_t, t_{t-l})}{Var(r_t)} ρl=Var(rt)Var(rt−l)Cov(rt,rt−l)=Var(rt)Cov(rt,tt−l)
这里用到了弱平稳序列的性质：
V
a
r
(
r
t
)
=
V
a
r
(
r
t
−
l
Var(r_t)=Var(r_{t-l}
Var(rt)=Var(rt−l
对一个平稳时间序列的样本 r t {r_t} rt， 1 ≤ t ≤ T 1\le t \le T 1≤t≤T,则间隔为 l l l的样本自相关系数的估计为：
ρ
^
l
=
∑
t
=
l
+
1
T
(
r
t
−
r
‾
)
(
r
t
−
l
−
r
‾
)
∑
t
=
1
T
(
r
t
−
r
‾
)
2
,
0
≤
l
≤
T
−
1
\hat{\rho}_l= \frac{\sum_{t=l+1}^T{(r_t-\overline{r})(r_{t-l}-\overline{r})}} {\sum_{t=1}^T{(r_t-\overline{r})^2}}, 0 \le l \le T-1
ρ^l=∑t=1T(rt−r)2∑t=l+1T(rt−r)(rt−l−r),0≤l≤T−1
则函数
ρ
^
1
,
ρ
^
2
,
ρ
^
3
.
.
.
.
\hat\rho_1,\hat\rho_2,\hat\rho_3 ....
ρ^1,ρ^2,ρ^3....称为
r
t
r_t
rt的样本自相关函数(ACF)
当自相关函数中所有的值都为0时，我们认为该序列是完全不相关的；因此，我们经常需要检验多个自相关系数是否为0。
混成检验
原假设：
H
0
:
ρ
1
=
…
=
ρ
m
=
0
H0:\rho_1=…=\rho_m=0
H0:ρ1=…=ρm=0
备择假设：
H
1
:
∃
i
∈
{
1
,
…
,
m
}
,
ρ
i
≠
0
H1:\exists i\in \{1,…,m\},\rho_i \neq 0
H1:∃i∈{1,…,m},ρi=0
混成检验统计量：
Q
(
m
)
=
T
(
T
+
2
)
∑
l
=
1
m
ρ
^
l
2
T
−
l
Q(m)=T(T+2)\sum_{l=1}^m{\frac{\hat \rho_l^2}{T-l}}
Q(m)=T(T+2)l=1∑mT−lρ^l2
Q
(
m
)
Q(m)
Q(m)渐进服从自由度为
m
m
m的
χ
2
\chi^2
χ2分布
决策规则：
Q
(
m
)
>
χ
α
2
,
拒绝
H
0
Q(m)>\chi^2_\alpha,拒绝H0
Q(m)>χα2,拒绝H0
即,
Q
(
m
)
Q(m)
Q(m)的值大于自由度为
m
m
m的卡方分布
100
(
1
−
α
)
100(1−\alpha)
100(1−α)分位点时，我们拒绝
H
0
H0
H0。
大部分软件会给出 Q ( m ) Q(m) Q(m)的p-value，则当p-value小于等于显著性水平 α \alpha α时拒绝 H 0 H0 H0。
示例：
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import tushare as ts
# tushare需要注册，qq邮箱
pro = ts.pro_api("7b318e083bd4a29877d393a1bf53597423c967278e71056f9c37e18a")
df = pro.daily(ts_code="600519.SH", start_date="20180101", end_date="20220101") # 茅台股价
df = df[["trade_date", "close"]]
df["trade_date"] = pd.to_datetime(df["trade_date"])
data = df.set_index("trade_date")
data["change_pct"] = ""
for i in range(0, len(data["close"]) - 1):
data.iloc[i, 1] = (
100 * (data["close"][i + 1] - data["close"][i]) / data["close"][i]
) # 计算每日收益率（涨跌幅）
data["change_pct"] = data["change_pct"][0:-1]
data["close_diff_1"] = data["close"].diff(1)
data = data[0:-1] # 由于changepct值最后一行缺失，因此去除最后一行
data["change_pct"] = data["change_pct"].astype("float64") # 使changepct数据类型和其他三列保持一致
data.plot(subplots=True, figsize=(12, 12))
array([<AxesSubplot: xlabel='trade_date'>,
<AxesSubplot: xlabel='trade_date'>,
<AxesSubplot: xlabel='trade_date'>], dtype=object)
import statsmodels.api as sm
close_data = data["close"] # 上证指数每日收盘价
m = 10 # 检验10个自相关系数
acf, q, p = sm.tsa.acf(close_data, nlags=m, qstat=True)
# np.r_是按列连接两个矩阵，就是把两矩阵上下相加，要求列数相等，类似于pandas中的concat()。
# np.c_是按行连接两个矩阵，就是把两矩阵左右相加，要求行数相等，类似于pandas中的merge()。
out = np.c_[range(1, 11), acf[1:], q, p]
output = pd.DataFrame(out, columns=["lag", "ACF", "Q", "P-value"])
output = output.set_index("lag")
print(output)
ACF Q P-value
lag
1.0 0.996748 968.671575 1.158897e-212
2.0 0.993510 1932.051494 0.000000e+00
3.0 0.990234 2890.076133 0.000000e+00
4.0 0.986840 3842.528221 0.000000e+00
5.0 0.983285 4789.108616 0.000000e+00
6.0 0.979778 5729.921467 0.000000e+00
7.0 0.976524 6665.464367 0.000000e+00
8.0 0.973601 7596.379456 0.000000e+00
9.0 0.970615 8522.554417 0.000000e+00
10.0 0.967580 9443.903276 0.000000e+00
我们取显著性水平为0.05，可以看出，所有的p-value都小于0.05；则我们拒绝原假设 H 0 H0 H0。因此，我们认为该序列是序列相关的。我们再来看看同期的日收益率序列：
change_data = data["change_pct"]
m = 10 # 检验10个自相关系数
acf, q, p = sm.tsa.acf(change_data, nlags=m, 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")
print(output)
AC Q P-value
lag
1.0 -0.009819 0.093996 0.759158
2.0 0.010291 0.197361 0.906032
3.0 -0.015774 0.440475 0.931760
4.0 0.009815 0.534696 0.970037
5.0 -0.014032 0.727460 0.981428
6.0 -0.055745 3.773007 0.707363
7.0 -0.003956 3.788358 0.803806
8.0 -0.026423 4.474037 0.812025
9.0 0.035592 5.719417 0.767639
10.0 0.003816 5.733751 0.837114
可以看出，p-value均大于显著性水平0.05。我们选择假设H0，即收益率序列没有显著的相关性。
自相关图与偏相关图详解
数据集：澳大利亚墨尔本10年最低气温
import pandas as pd
from matplotlib import pyplot as plt
data = pd.read_csv(
"daily-minimum-temperatures-in-me.csv", parse_dates=["Date"], index_col="Date"
)
data.head()
data.plot(figsize=(12, 6))
<AxesSubplot: xlabel='Date'>
相关与自相关
统计中的相关性说的是两个变量间的相关程度。我们假设如果每个变量都符合正态分布，可以使用皮尔逊相关系数（Pearson correlation coefficient）来统计变量间的相关性。皮尔逊相关系数是-1和1之间的数字分别描述负相关或正相关。值为零表示无相关。
自相关（Autocorrelation），也叫序列相关，是一个信号于其自身在不同时间点的互相关。非正式地来说，它就是两次观察之间的相似度对它们之间的时间差的函数。它是找出重复模式（如被噪声掩盖的周期信号），或识别隐含在信号谐波频率中消失的基频的数学工具。自相关函数在不同的领域，定义不完全等效。在某些领域，自相关函数等同于自协方差。
由于时间序列的相关性与之前的相同系列的值进行了计算，这被称为序列相关或自相关。一个时间序列的自相关系数被称为自相关函数，或简称ACF。这个图被称为相关图或自相关图。
以下是利用statsmodels库中使用plot_acf()函数计算和绘制“每日最低气温”自相关图：
import statsmodels.api as sm
sm.graphics.tsa.plot_acf(data, lags=150)
plt.figure(figsize=(12, 6))
plt.show()
<Figure size 1200x600 with 0 Axes>
这个示例创建了一个2D的平面图，显示沿x轴的延迟值以及y轴上的相关性(-1到1之间)。置信区间被画成圆锥形。默认情况下，置信区间这被设置为95%。
默认情况下，图中显示的是所有的延迟值，这让真实情况掩藏了，我们可以将x轴上的延迟值限制为50，让图更容易看懂。
sm.graphics.tsa.plot_acf(data, lags=50)
plt.figure(figsize=(12, 6))
plt.show()
<Figure size 1200x600 with 0 Axes>
偏自相关函数（PACF）
偏自相关函数用来度量暂时调整所有其他较短滞后的项 ( y t – 1 , y t – 2 , … , y t – k – 1 ) (y_{t–1},y_{t–2},…,y_{t–k–1}) (yt–1,yt–2,…,yt–k–1) 之后，时间序列中以 k k k 个时间单位（ y t y_t yt和 y t – k y_{t–k} yt–k）分隔的观测值之间的相关。
偏自相关是剔除干扰后时间序列观察与先前时间步长时间序列观察之间关系的总结。在滞后 k k k处的偏自相关是在消除由于较短滞后条件导致的任何相关性的影响之后产生的相关性。一项观察的自相关和在先验时间步上的观测包括直接相关和间接相关。这些间接相关是线性函数观察（这个观察在两个时间步长之间）的相关。偏自相关函数试图移除这些间接相关。
下面的示例使用statsmodels库中的plot_pacf()来计算和绘制最低每日温度数据集里的前50个滞后的偏自相关函数：
sm.graphics.tsa.plot_pacf(data, lags=50)
plt.figure(figsize=(12, 6))
plt.show()
D:\soft\tool\Anaconda3\envs\py38\lib\site-packages\statsmodels\graphics\tsaplots.py:348: FutureWarning: The default method 'yw' can produce PACF values outside of the [-1,1] interval. After 0.13, the default will change tounadjusted Yule-Walker ('ywm'). You can use this method now by setting method='ywm'.
warnings.warn(
<Figure size 1200x600 with 0 Axes>
ACF和PACF图的直观认识
时间序列的自相关函数和偏自相关函数的平面图描述了完全不同的情形。我们可以使用ACF和PACF图来探索一些信息。
自回归直观认识(intuition)
考虑由自回归(AR)过程产生的滞后时间为 k k k的时间序列。我们知道，ACF描述了一个观测值与另一个观测值之间的自相关，包括直接和间接的相关性信息。这意味着我们可以预期AR(k)时间序列的ACF强大到(如同使用了) k k k的滞后，并且这种关系的惯性将继续到之后的滞后值，随着效应被削弱而在某个点上缩小到没有。我们知道，PACF只描述观测值与其滞后(lag)之间的直接关系。这表明，超过 k k k的滞后值(lag value)不会再有相关性。这正是ACF和PACF图对AR(k)过程的预期。
滑动平均直观认识(Moving Average Intuition)
考虑由滑动平均(MA)过程产生的滞后(lag)时间为 k k k的时间序列。请记住，滑动平均过程是先前预测的残留偏差的时间序列的自回归模型。考虑滑动平均模型的另一种方法是根据最近预测的错误修正未来的预测。我们期望MA(k)过程的ACF与最近的lag值之间的关系显示出强烈的相关性，然后急剧下降到低或者无相关性。根据定义，这解释了整个过程是如何产生的。对于PACF，我们预计图会显示与滞后(lag)的关系，以及滞后(lag)之前的相关。再次强调，这正是MAF(k)过程的ACF和PACF图的预期。
自相关和偏自相关用于测量当前序列值和过去序列值之间的相关性，并指示预测将来值时最有用的过去序列值。了解了此内容，您就可以确定ARIMA模型中过程的顺序。更具体来说，
自相关函数 (ACF)。延迟为
k
k
k时，这是相距
k
k
k 个时间间隔的序列值之间的相关性。
偏自相关函数 (PACF)。延迟为
k
k
k 时，这是相距
k
k
k 个时间间隔的序列值之间的相关性，同时考虑了间隔之间的值。
截尾与拖尾是什么鬼？
截尾是指时间序列的自相关函数（ACF）或偏自相关函数（PACF）在某阶后均为0的性质（比如AR的PACF）；
拖尾是ACF或PACF并不在某阶后均为0的性质（比如AR的ACF）。
截尾：在大于某个常数
k
k
k后快速趋于0为
k
k
k阶截尾
拖尾：始终有非零取值，不会在
k
k
k大于某个常数后就恒等于零(或在0附近随机波动)
自回归移动平均模型(ARMA(p，q))是时间序列中最为重要的模型之一，它主要由两部分组成： AR代表p阶自回归过程，MA代表q阶移动平均过程，其公式如下：
z
t
=
φ
1
z
t
−
1
+
φ
2
z
t
−
2
+
.
.
.
+
φ
p
z
t
−
p
+
α
t
−
θ
1
α
t
−
1
−
.
.
.
−
θ
p
α
t
−
p
z_t=\varphi_1z_{t-1}+\varphi_2z_{t-2}+...+\varphi_pz_{t-p}+\alpha_t-\theta_1\alpha_{t-1}-...-\theta_p\alpha_{t-p}
zt=φ1zt−1+φ2zt−2+...+φpzt−p+αt−θ1αt−1−...−θpαt−p
模型的简化形式为：
φ
p
(
B
)
z
t
=
θ
q
(
B
)
α
(
t
)
\varphi_p(B)z_t=\theta_q(B)\alpha(t)
φp(B)zt=θq(B)α(t)
其中：
φ
p
(
B
)
=
1
−
φ
1
B
−
φ
2
B
2
−
.
.
.
−
φ
p
B
p
\varphi_p(B)=1-\varphi_1B-\varphi_2B^2-...-\varphi_pB^p
φp(B)=1−φ1B−φ2B2−...−φpBp
θ q ( B ) = 1 − θ 1 B − θ 2 B 2 − . . . − θ q B q \theta_q(B)=1-\theta_1B-\theta_2B^2-...-\theta_qB^q θq(B)=1−θ1B−θ2B2−...−θqBq
依据模型的形式、特性及自相关和偏自相关函数的特征，总结如下：
在时间序列中，ARIMA模型是在ARMA模型的基础上多了差分的操作。
白噪声序列与线性时间序列
白噪声序列
随机变量 X ( t ) ( t = 1 , 2 , 3 … ) X(t) (t=1,2,3…) X(t)(t=1,2,3…)，如果是由一个不相关的随机变量的序列构成的，即对于所有 S S S不等于 T T T，随机变量 X t X_t Xt和 X s X_s Xs的协方差为零，则称其为纯随机过程。对于一个纯随机过程来说，若其期望和方差均为常数，则称之为白噪声过程。白噪声过程的样本实称成为白噪声序列，简称白噪声。之所以称为白噪声，是因为他和白光的特性类似，白光的光谱在各个频率上有相同的强度，白噪声的谱密度在各个频率上的值相同。
线性时间序列
时间序列
{
r
t
}
\{r_t\}
{rt}，如果能写成：
r
t
=
μ
+
∑
i
=
0
∞
Ψ
i
a
t
−
i
r_t=\mu+\sum_{i=0}^\infty{\Psi_ia_{t-i}}
rt=μ+i=0∑∞Ψiat−i
μ
\mu
μ为
r
t
r_t
rt的均值，
Ψ
0
=
1
\Psi_0=1
Ψ0=1,
{
a
t
}
\{a_t\}
{at}为白噪声序列，则我们称
{
r
t
}
\{r_t\}
{rt}为线性序列。其中
a
t
a_t
at称为在
t
t
t时刻的新息(innovation)或扰动(shock)。
很多时间序列具有线性性，即是线性时间序列，相应的有很多线性时间序列模型，例如接下来要介绍的AR、MA、ARMA，都是线性模型，但并不是所有的金融时间序列都是线性的。对于弱平稳序列，我们利用白噪声的性质很容易得到 r t r_t rt的均值和方差：
E
(
r
t
)
=
μ
,
V
a
r
(
r
t
)
=
σ
a
2
∑
i
=
0
∞
Ψ
i
2
E(r_t)=\mu,Var(r_t)=\sigma^2_a\sum_{i=0}^\infty{\Psi_i^2}
E(rt)=μ,Var(rt)=σa2i=0∑∞Ψi2
σ
a
2
\sigma^2_a
σa2为
a
t
a_t
at的方差。因为
V
a
r
(
r
t
)
Var(r_t)
Var(rt)一定小于正无穷，因此
{
Ψ
i
2
}
\{\Psi^2_i\}
{Ψi2}必须是收敛序列，因此满足
i
→
∞
i\rightarrow\infty
i→∞时，
Ψ
i
2
→
0
\Psi_i^2\rightarrow 0
Ψi2→0。即，随着
i
i
i的增大，远处的扰动
a
t
−
i
a_{t−i}
at−i对
r
t
r_t
rt的影响会逐渐消失。
自回归模型(Auto Regressive, AR)
AR 模型的理解
让我们用下面的例子来理解 AR 模型：一个国家当前的GDP称
x
(
t
)
x(t)
x(t)取决于去年的GDP，即
x
(
t
−
1
)
x(t-1)
x(t−1)。一个国家在一个财政年度的GDP取决于前一年建立的制造厂/服务，以及本年度新成立工厂/服务。但是GDP的主要组成部分是前一部分，因此，我们可以把 GDP 的等式写成:
x
(
t
)
=
α
∗
x
(
t
−
1
)
+
e
r
r
o
r
r
(
t
)
x(t)=\alpha * x(t-1) + errorr(t)
x(t)=α∗x(t−1)+errorr(t)
这个方程称为AR公式。公式表示下一个点完全依赖与前面一个点。 α \alpha α是我们寻求的系数，以便最小化误差函数。注意， x ( t − 1 ) x(t-1) x(t−1)同样依赖 x ( t − 2 ) x(t-2) x(t−2)。 因此，波动在未来会逐渐消退。
例如，假设
x
(
t
)
x(t)
x(t)是一个城市在某一天出售的果汁瓶的数量。 在冬天，很少有供应商购买果汁瓶。突然，在一个特定的日子里，温度上升，果汁的需求猛增到1000个。然而，几天后，气候又变得寒冷起来。但是，知道人们在炎热的日子里已经习惯了喝果汁，在寒冷的日子里仍有50% 的人还在喝果汁。接下来的几天，这个比例下降到25% (50%的50%)，然后在数天之后逐渐变小。下面的图表说明了AR系列的惯性特性：
在上节中，我们计算了部分数据端的ACF，看表可知间隔为1时自相关系数是显著的。这说明在 t − 1 t−1 t−1时刻的数据 r t − 1 r_{t−1} rt−1，在预测 t t t时刻 r t r_t rt时可能是有用的！
根据这点我们可以建立下面的模型：
r
t
=
ϕ
0
+
ϕ
1
r
t
−
1
+
a
t
r_t=\phi_0+\phi_1 r_{t-1}+a_t
rt=ϕ0+ϕ1rt−1+at
其中
{
a
t
}
\{a_t\}
{at}是白噪声序列，这个模型与简单线性回归模型有相同的形式，这个模型也叫做一阶自回归(AR)模型，简称AR(1)模型。从AR(1)很容易推广到AR§模型：
r
t
=
ϕ
0
+
ϕ
1
r
t
−
1
+
.
.
.
+
ϕ
p
r
t
−
p
+
a
t
r_t=\phi_0+\phi_1 r_{t-1}+...+\phi_p r_{t-p}+a_t
rt=ϕ0+ϕ1rt−1+...+ϕprt−p+at
AR§模型的特征根及平稳性检验
我们先假定序列是弱平稳的，则有：
E
(
r
t
)
=
μ
,
V
a
r
(
r
t
)
=
γ
0
,
C
o
v
(
r
t
,
t
t
−
j
)
=
γ
j
E(r_t)=\mu , Var(r_t)=\gamma_0 , Cov(r_t,t_{t-j})=\gamma_j
E(rt)=μ,Var(rt)=γ0,Cov(rt,tt−j)=γj
其中
μ
,
γ
0
\mu,\gamma_0
μ,γ0是常数。因为
{
a
t
}
\{a_t\}
{at}是白噪声序列，因此有：
E
(
a
t
)
=
0
,
V
a
r
(
a
t
)
=
σ
a
2
E(a_t)=0,Var(a_t)=\sigma_a^2
E(at)=0,Var(at)=σa2
所以有：
E
(
r
t
)
=
ϕ
0
+
ϕ
1
E
(
r
t
−
1
)
+
ϕ
2
E
(
r
t
−
2
)
+
.
.
.
+
ϕ
p
E
(
r
t
−
p
)
E(r_t)=\phi_0+\phi_1E(r_{t-1})+\phi_2E(r_{t-2})+...+\phi_pE(r_{t-p})
E(rt)=ϕ0+ϕ1E(rt−1)+ϕ2E(rt−2)+...+ϕpE(rt−p)
根据平稳性的性质，又有
E
(
r
t
)
=
E
(
r
t
−
1
)
=
.
.
.
=
μ
E(r_t)=E(r_{t−1})=...=\mu
E(rt)=E(rt−1)=...=μ，从而：
μ
=
ϕ
0
+
ϕ
1
+
μ
+
ϕ
2
μ
+
.
.
.
+
ϕ
p
μ
\mu=\phi_0+\phi_1+\mu+\phi_2\mu+...+\phi_p\mu
μ=ϕ0+ϕ1+μ+ϕ2μ+...+ϕpμ
E
(
r
t
)
=
μ
=
ϕ
0
1
−
ϕ
1
−
ϕ
2
−
.
.
.
−
ϕ
p
E(r_t)=\mu=\frac{\phi_0}{1-\phi_1-\phi_2-...-\phi_p}
E(rt)=μ=1−ϕ1−ϕ2−...−ϕpϕ0
公式中假定分母不为0， 我们将下面的方程称为特征方程:
1
−
ϕ
1
x
−
ϕ
2
x
2
−
.
.
.
−
ϕ
p
x
p
=
0
1-\phi_1x-\phi_2x^2-...-\phi_px^p=0
1−ϕ1x−ϕ2x2−...−ϕpxp=0
该方程所有解的倒数称为该模型的特征根，如果所有的特征根的模都小于1，则该AR§序列是平稳的。
下面我们就用该方法，检验日收益率序列的平稳性：
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import tushare as ts
# tushare需要注册，qq邮箱
pro = ts.pro_api("7b318e083bd4a29877d393a1bf53597423c967278e71056f9c37e18a")
df = pro.daily(ts_code="600519.SH", start_date="20180101", end_date="20190101") # 茅台股价
df = df[["trade_date", "close"]]
df["trade_date"] = pd.to_datetime(df["trade_date"])
data = df.set_index("trade_date")
data["change_pct"] = ""
for i in range(0, len(data["close"]) - 1):
data.iloc[i, 1] = (
100 * (data["close"][i + 1] - data["close"][i]) / data["close"][i]
) # 计算每日收益率（涨跌幅）
data["change_pct"] = data["change_pct"][0:-1]
data["close_diff_1"] = data["close"].diff(1)
data = data[0:-1] # 由于changepct值最后一行缺失，因此去除最后一行
data["change_pct"] = data["change_pct"].astype("float64") # 使changepct数据类型和其他三列保持一致
import statsmodels.api as sm
from statsmodels.tsa.ar_model import AutoReg as AR
change_data = data["change_pct"] # 载入收益率序列
change_data = np.array(change_data, dtype=np.float)
model = AR(change_data, 20)
results_AR = model.fit()
plt.figure(figsize=(12, 6))
plt.plot(change_data, "b", label="changepct")
plt.plot(results_AR.fittedvalues, "r", label="AR model")
plt.legend()
C:\Users\suny\AppData\Local\Temp\ipykernel_1880\1463261479.py:5: DeprecationWarning: `np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
change_data = np.array(change_data, dtype=np.float)
<matplotlib.legend.Legend at 0x186da0eeeb0>
print(len(results_AR.roots)) # 模型阶数
20
画出模型的特征根，来检验平稳性：
r1 = 1
theta = np.linspace(0, 2 * np.pi, 360)
x1 = r1 * np.cos(theta)
y1 = r1 * np.sin(theta)
plt.figure(figsize=(6, 6))
plt.plot(x1, y1, "k") # 画单位圆
roots = 1 / results_AR.roots # 计算results_AR.roots后取倒数
for i in range(len(roots)):
plt.plot(roots[i].real, roots[i].imag, ".r", markersize=8) # 画特征根
plt.show()
可以看出，所有特征根都在单位圆内，则序列为平稳的。
AR模型的阶数确定
对于线性回归模型的阶数，一般认为随着阶数的升高，会提高模型的精度。但是过高的阶数会造成模型过于复杂。下面介绍两种确定AR模型阶数的方法：偏相关函数（pACF）和信息准则函数（AIC、BIC、HQIC）。
偏相关函数(pACF)
首先引入阶数依次升高的AR模型：
r
t
=
ϕ
0
,
1
+
ϕ
1
,
1
r
t
−
1
+
ε
1
t
r_t=\phi_{0,1}+\phi_{1,1}r_{t-1}+\varepsilon_{1t}
rt=ϕ0,1+ϕ1,1rt−1+ε1t
r
t
=
ϕ
0
,
2
+
ϕ
1
,
2
r
t
−
1
+
ϕ
2
,
2
r
t
−
2
+
ε
2
t
r_t=\phi_{0,2}+\phi_{1,2}r_{t-1}+\phi_{2,2}r_{t-2}+\varepsilon_{2t}
rt=ϕ0,2+ϕ1,2rt−1+ϕ2,2rt−2+ε2t
r
t
=
ϕ
0
,
3
+
ϕ
1
,
3
r
t
−
1
+
ϕ
2
,
3
r
t
−
2
+
ϕ
3
,
4
r
t
−
3
+
ε
3
t
r_t=\phi_{0,3}+\phi_{1,3}r_{t-1}+\phi_{2,3}r_{t-2}+\phi_{3,4}r_{t-3}+\varepsilon_{3t}
rt=ϕ0,3+ϕ1,3rt−1+ϕ2,3rt−2+ϕ3,4rt−3+ε3t
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
...................
...................
等式右边第一项为常数项，最后一项为误差项， ϕ ^ i , i \hat\phi_{i,i} ϕ^i,i为 { r t } \{r_t\} {rt}的间隔为 i i i的样本偏自相关函数。 ϕ ^ i , i \hat\phi_{i,i} ϕ^i,i表示的是在第 i − 1 i-1 i−1行的基础上添加的一个阶数对 { r t } \{r_t\} {rt}的贡献。对于AR模型，当超过某个值 p p p时，偏自相关函数应逐渐趋近于0，这既符合我们的直观认识，也有相应的理论证明可以找到。
对于偏相关函数的介绍，这里不详细展开，重点介绍一个性质：AR§序列的样本偏相关函数是 p p p步截尾的。所谓截尾，就是快速收敛应该是快速的降到几乎为0或者在置信区间以内。
下面对2018年茅台股价涨跌幅作出偏自相关函数图，在5%的显著水平下，取置信区间为
2
n
=
2
243
\frac{2}{\sqrt{n}}=\frac{2}{\sqrt{243}}
n2=2432
查看pACF图：
fig = plt.figure(figsize=(12, 6))
ax1 = fig.add_subplot(111)
fig = sm.graphics.tsa.plot_pacf(change_data, lags=40, ax=ax1)
D:\soft\tool\Anaconda3\envs\py38\lib\site-packages\statsmodels\graphics\tsaplots.py:348: FutureWarning: The default method 'yw' can produce PACF values outside of the [-1,1] interval. After 0.13, the default will change tounadjusted Yule-Walker ('ywm'). You can use this method now by setting method='ywm'.
warnings.warn(
很遗憾，从上图中并不太容易确定一个明显的阶数 p p p。不过确实可以看出，随着时间的增加，偏相关函数逐渐趋近于0。这个结果给了我们两点思路：首先偏相关函数趋近于0，说明运用AR模型模拟该序列是有一定的合理性；其次上图很难确定阶数 p p p或阶数 p p p很大，说明模型可能有MA（滑动平均，将在后面介绍）的部分。它至少给了我们一些感性的认识。
信息准则函数
目前比较常用的信息准则函数包括赤池信息量（AIC），贝叶斯信息量（BIC）和汉南-奎因信息量（HQIC）。他们都基于似然函数L和参数的数量k，其中n表示样本大小。由于Python中有相应的函数，所以公式仅作了解：
A
I
C
=
−
2
l
n
L
+
2
k
AIC=−2lnL+2k
AIC=−2lnL+2k
B
I
C
=
−
2
l
n
L
+
l
n
n
⋅
k
BIC=−2lnL+lnn⋅k
BIC=−2lnL+lnn⋅k
H
Q
I
C
=
−
2
l
n
L
+
l
n
(
l
n
n
)
⋅
k
HQIC=−2lnL+ln(lnn)⋅k
HQIC=−2lnL+ln(lnn)⋅k
我们的目标是，选择阶数
k
k
k使得信息量的值最小。以上这么多可供选择的模型，我们通常采用AIC法则。我们知道：增加自由参数的数目提高了拟合的优良性，AIC鼓励数据拟合的优良性但是尽量避免出现过度拟合(Overfitting)的情况。所以优先考虑的模型应是AIC值最小的那一个。
下面我们来测试下3种准则下确定的 p p p，仍然用日收益率序列。为了减少计算量，我们只计算间隔前7个看看效果。
aicList = []
bicList = []
hqicList = []
for i in range(1, 8): # 从1阶开始算
order = (i, 0, 0) # 使用ARMA模型，其中MA阶数为0，只考虑了AR阶数。
chgdataModel = sm.tsa.arima.ARIMA(change_data, order=order).fit()
aicList.append(chgdataModel.aic)
bicList.append(chgdataModel.bic)
hqicList.append(chgdataModel.hqic)
plt.figure(figsize=(10, 6))
plt.plot(aicList, "ro--", label="aic value")
plt.plot(bicList, "bo--", label="bic value")
plt.plot(hqicList, "ko--", label="hqic value")
plt.legend(loc=0)
<matplotlib.legend.Legend at 0x186da35feb0>
计算信息准则函数的时候，我们发现随着阶数的增加，运算速度越来越慢，所以只取了前7阶，可以看出AIC给出的阶数为0，BIC和HQIC给出的阶数为0。当然，前面提到只用AR模型预测该序列是不准确的，且只计算前7阶的结果也不足以判断最佳阶数。当然，利用上面的方法逐个计算是很耗时间的，实际上，有函数可以直接按照准则计算出合适的order，这个是针对ARMA模型的，我们后续再讨论。
AR模型的检验
如果模型是充分的，其残差序列应该是白噪声，根据我们第一章里介绍的混成检验，可以用来检验残差与白噪声的接近程度。先求出残差序列：
delta = results_AR.fittedvalues - change_data[20:] # 取得误差项
plt.figure(figsize=(10, 6))
# plt.plot(change_data[17:],label='original value')
# plt.plot(results_AR.fittedvalues,label='fitted value')
plt.plot(delta, "r", label=" residual error")
plt.legend(loc=0)
<matplotlib.legend.Legend at 0x186da3bc370>
然后我们检查它是不是接近白噪声序列：
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", "ACF", "Q", "P-value"])
output = output.set_index("lag")
print(output)
ACF Q P-value
lag
1.0 -0.008982 0.018154 0.892821
2.0 0.001738 0.018837 0.990626
3.0 -0.020685 0.115988 0.989852
4.0 0.018844 0.196987 0.995457
5.0 0.006205 0.205810 0.999050
6.0 -0.004983 0.211526 0.999818
7.0 -0.018759 0.292917 0.999908
8.0 -0.030360 0.507104 0.999859
9.0 0.021833 0.618396 0.999925
10.0 -0.019082 0.703808 0.999966
观察p-value可知，该序列可以认为没有相关性，近似得可以认为残差序列接近白噪声。
拟合优度
调整后的拟合优度检验：
A
d
j
R
2
=
1
−
σ
^
ϵ
2
σ
^
r
t
2
Ad_jR^2=1-\frac{\hat\sigma^2_\epsilon}{\hat\sigma^2_{r_t}}
AdjR2=1−σ^rt2σ^ϵ2
右边分子代表残差的方差，分母代表
r
t
r_t
rt的方差。
A
d
j
R
2
Ad_jR^2
AdjR2值在0到1之间，越大说明拟合效果越好。计算该序列的
A
d
j
R
²
Ad_jR^²
AdjR²值，发现效果并不好，这也符合我们的预期。
下面我们计算之前对日收益率的AR模型的拟合优度：
adjR = 1 - delta.var() / change_data[20:].var()
print(adjR)
0.06941307925186757
可以看出，模型的拟合程度并不好，当然，这并不重要，也许是这个序列并不适合用AR模型拟合。
AR模型的预测
我们首先得把原来的样本分为训练集和测试集，再来看预测效果，还是以之前的数据为例：
train = change_data[:-10]
test = change_data[-10:]
output = sm.tsa.arima.ARIMA(train, order=(20, 0, 0)).fit()
output.predict()
predicts = output.predict(216, 225, dynamic=True) # 共226个数据
compare = pd.DataFrame()
compare["original"] = change_data[-10:]
compare["predict"] = predicts
print(compare)
compare.plot()
original predict
0 1.608146 0.440614
1 0.388352 0.531311
2 -1.726237 0.742988
3 1.406797 0.654107
4 -0.406002 -0.171440
5 -3.883607 0.552274
6 -1.830801 -0.278129
7 -0.174712 0.727402
8 -2.877610 -0.167412
9 -1.677702 -0.213557
<AxesSubplot: >
该模型的预测结果是很差。那么我们是不是可以通过其他模型获得更好的结果呢？ 我们将在下一部分继续介绍。
滑动平均模型（Moving Average, MA）
MA模型的理解
让我们以另一个例子来理解移动平均时间序列模型：一个公司生产了一批类型在市场上很容易买到的包，在竞争激烈的市场。这个包的销量在很多天都是0，于是，在某一天，这家公司做了一个实验，生产了一个不同类型的包，这种包在市场上任何地方找不到，因此，他可以一下子卖掉了1000个。这个包的需求特别高，很快库存快要完了。这天结束了还有100个包没卖掉。我们把这个误差成为时间点误差。接下来的几天仍有几个客户购买这种包。随着时间的推移，这个袋子已经失去了它的吸引力。下面通过一个简单的公式来描述这个场景：
x
(
t
)
=
b
e
t
a
∗
e
r
r
o
r
(
t
−
1
)
+
e
r
r
o
r
(
t
)
x(t)=beta * error(t-1) + error(t)
x(t)=beta∗error(t−1)+error(t)
如果我们尝试绘制这个图表，它会看起来像这样：
注意到 MA 和 AR 模型的区别了吗？ 在 MA 模型中，噪声/冲击随着时间的推移而迅速消失。 AR模型对冲击具有持久的影响。
这里我们直接给出MA(q)模型的形式：
r
t
=
c
0
+
a
t
−
θ
1
a
t
−
1
−
.
.
.
−
θ
q
a
t
−
q
r_t=c_0+a_t-\theta_1a_{t-1}-...-\theta_qa_{t-q}
rt=c0+at−θ1at−1−...−θqat−q
c
0
c_0
c0为一个常数项。这里的
a
t
a_t
at,是AR模型t时刻的扰动或者说新息，则可以发现，该模型，使用了过去
q
q
q个时期的随机干扰或预测误差来线性表达当前的预测值。
MA模型的性质
平稳性
MA模型总是弱平稳的，因为他们是白噪声序列（残差序列）的有限线性组合。因此，根据弱平稳的性质可以得出两个结论：
E
(
r
t
)
=
c
0
E(r_t)=c_0
E(rt)=c0
V
a
r
(
r
t
)
=
(
1
+
θ
1
2
+
θ
2
2
+
.
.
.
+
θ
q
2
)
σ
a
2
Var(r_t)=(1+\theta_1^2+\theta_2^2+...+\theta_q^2)\sigma^2_a
Var(rt)=(1+θ12+θ22+...+θq2)σa2
自相关函数
对 q q q阶的MA模型，其自相关函数ACF总是 q q q步截尾的。 因此MA(q)序列只与其前q个延迟值线性相关，从而它是一个“有限记忆”的模型。这一点可以用来确定模型的阶次，后面会介绍。
可逆性
当满足可逆条件的时候，MA(q)模型可以改写为AR§模型。这里不进行推导，给出1阶和2阶MA的可逆性条件。
1阶： ∣ θ 1 ∣ < 1 |\theta_1|<1 ∣θ1∣<1
2阶： ∣ θ 2 ∣ < 1 , θ 1 + θ 2 < 1 |\theta_2|<1,\theta_1+\theta_2<1 ∣θ2∣<1,θ1+θ2<1
MA模型的的阶数确定
之前我们知道，通过计算偏自相关函数（pACF）可以粗略确定AR模型的阶数。而对于MA模型，需要采用自相关函数（ACF）进行定阶。
fig2 = plt.figure(figsize=(12, 6))
ax2 = fig2.add_subplot(111)
fig2 = sm.graphics.tsa.plot_acf(change_data, lags=40, ax=ax2)
可以发现还是不太明显，这里以6阶为列，作出MA模型与原日涨跌幅的图像：
results_MA = sm.tsa.arima.ARIMA(
change_data, order=(0, 0, 6)
).fit() # 取ARMA模型中AR的阶数为0，MA的阶数为6
plt.figure(figsize=(12, 6))
plt.plot(change_data, "b", label="changepct")
plt.plot(results_MA.fittedvalues, "r", label="MA model")
plt.legend()
<matplotlib.legend.Legend at 0x186dadd0b80>
MA模型的预测
因此这里就建立6阶的模型作为示例。我们用最后10个数据作为out-sample的样本，用来对比预测值。
train = change_data[:-10]
test = change_data[-10:]
chgdataModel_MA = sm.tsa.arima.ARIMA(
train, order=(0, 0, 6)
).fit() # 取ARMA模型中AR的阶数为0，MA的阶数为6
delta = chgdataModel_MA.fittedvalues - train
score = 1 - delta.var() / train.var()
print(score)
0.026373926167318595
远小于1，拟合效果不好。然后我们用建立的模型进行预测最后10个数据：
predicts_MA = chgdataModel_MA.predict(216, 225, dynamic=True)
compare_MA = pd.DataFrame()
compare_MA["original"] = test
compare_MA["predict_MA"] = predicts_MA
print(compare_MA)
original predict_MA
0 1.608146 0.661439
1 0.388352 0.550726
2 -1.726237 0.248529
3 1.406797 0.805170
4 -0.406002 -0.292676
5 -3.883607 0.564484
6 -1.830801 0.144742
7 -0.174712 0.144742
8 -2.877610 0.144742
9 -1.677702 0.144742
可以看出，建立的模型效果很差，就算只看涨跌方向，正确率也不足。该模型不适用于原数据。没关系，如果真的这么简单能预测股价日涨跌才奇怪了~关于MA的内容只做了简单介绍，下面主要介绍ARMA模型。
自回归-滑动平均模型（ARMA）
AR模型与MA模型的不同
AR和MA 模型之间的主要区别是基于不同时间点的时间序列对象之间的相关性。 x ( t ) x(t) x(t)与 x ( t − n ) x(t-n) x(t−n)的相关性总为0，这直接源于 x ( t ) x(t) x(t)和 x ( t − n ) x(t-n) x(t−n)之间的协方差对于MA模型来说是零。然而，AM模型中 x ( t ) x(t) x(t)与 x ( t − 1 ) x(t-1) x(t−1)的相关性随着时间的推移变得越来越小。不管是否有 AR 模型或 MA 模型，这种差异都可以被利用。
在有些应用中，我们需要高阶的AR或MA模型才能充分地描述数据的动态结构，这样问题会变得很繁琐。为了克服这个困难，提出了自回归滑动平均(ARMA)模型。基本思想是把AR和MA模型结合在一起，使所使用的参数个数保持很小。模型的形式为：
r
t
=
ϕ
0
+
∑
i
=
1
p
ϕ
i
r
t
−
i
+
a
t
+
∑
i
=
1
q
θ
i
a
t
−
i
r_t=\phi_0+\sum_{i=1}^{p}{\phi _ir_{t-i}}+a_t+\sum_{i=1}^q{\theta_ia_{t-i}}
rt=ϕ0+i=1∑pϕirt−i+at+i=1∑qθiat−i
其中， { a t } \{a_t\} {at}为白噪声序列， p p p和 q q q都是非负整数。AR和MA模型都是ARMA(p,q)的特殊形式。
利用向后推移算子B（后移算子B，即上一时刻），上述模型可写成:
(
1
−
ϕ
1
B
−
.
.
.
−
ϕ
p
B
p
)
r
t
=
θ
0
+
(
1
−
θ
1
B
−
.
.
.
−
θ
q
B
q
)
a
t
(1-\phi_1B-...-\phi_pB^p)r_t=\theta_0+(1-\theta_1B-...-\theta_qB^q)a_t
(1−ϕ1B−...−ϕpBp)rt=θ0+(1−θ1B−...−θqBq)at
这时候我们求
r
t
r_t
rt的期望，得到：
E
(
r
t
)
=
ϕ
0
1
−
ϕ
0
−
.
.
.
−
ϕ
p
E(r_t)=\frac{\phi_0}{1-\phi_0-...-\phi_p}
E(rt)=1−ϕ0−...−ϕpϕ0
和上期我们的AR模型一样。因此有着相同的特征方程：
1
−
ϕ
1
x
−
ϕ
2
x
2
−
.
.
.
−
ϕ
p
x
p
=
0
1-\phi_1x-\phi_2x^2-...-\phi_px^p=0
1−ϕ1x−ϕ2x2−...−ϕpxp=0
该方程所有解的倒数称为该模型的特征根，如果所有的特征根的模都小于1，则该ARMA模型是平稳的。有一点很关键：ARMA模型的应用对象应该为平稳序列！下面的步骤都是建立在假设原序列平稳的条件下的。
一旦我们得到平稳时间序列，我们必须回答：是AR 还是MA 过程？
ACF是不同滞后函数之间的总关联图。例如，在GDP问题上，时间点T的GDP是 x ( t ) x(t) x(t)。我们对 x ( t ) x(t) x(t)与 x ( t − 1 ) x(t-1) x(t−1)、 x ( t − 2 ) x(t-2) x(t−2)等的相关性很感兴趣。现在让我们来反思一下上面所学的内容。在移动平均序列的滞后 n n n中， x ( t ) x(t) x(t)与 x ( t − n − 1 ) x(t-n-1) x(t−n−1)之间不存在任何相关。因此，总相关图在第 n n n 次滞后处截止。因此，很容易找到一个MA系列的滞后。对于AR 序列，这种相关性将逐渐下降而没有任何截止值。
如果发现每个滞后的部分相关，则在AR 级数的程度之后会被切断。例如，如果我们有AR系列，如果我们排除第一滞后 x ( t − 1 ) x(t-1) x(t−1)的影响，则我们的第二滞后 x ( t − 2 ) x(t-2) x(t−2)与 x ( t ) x(t) x(t)无关。因此，偏相关函数（PACF）在第一滞后急剧下降。下面是一些例子来澄清你对这个概念的任何怀疑：
上面的蓝线显示出明显不同于零的值。显然，上面的图在第二滞后之后在PACF 曲线上有一个截止，这意味着这主要是AR过程。
显然，上面的图在第二滞后之后在ACF 曲线上有一个截止，这意味着这主要是MA过程。
ARMA模型的阶数确定
PACF、ACF 判断模型阶数
我们通过观察PACF和ACF截尾，分别判断p、q的值。(限定滞后阶数40)
fig = plt.figure(figsize=(12, 12))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(change_data, lags=30, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(change_data, lags=30, ax=ax2)
D:\soft\tool\Anaconda3\envs\py38\lib\site-packages\statsmodels\graphics\tsaplots.py:348: FutureWarning: The default method 'yw' can produce PACF values outside of the [-1,1] interval. After 0.13, the default will change tounadjusted Yule-Walker ('ywm'). You can use this method now by setting method='ywm'.
warnings.warn(
可以看出，模型的阶次应该为(6,24)。然而，这么高的阶次建模的计算量是巨大的。为什么不再限制滞后阶数小一些？如果lags设置为25，20或者更小时，阶数为(0,0)，显然不是我们想要的结果。
综合来看，由于计算量太大，在这里就不使用(6,24)建模了。采用另外一种方法确定阶数。
信息准则定阶
关于信息准则，已有过一些介绍。下面，我们分别应用以上3种法则，为我们的模型定阶，数据仍然是日涨跌幅序列。为了控制计算量，我们限制AR最大阶不超过6，MA最大阶不超过4。 但是这样带来的坏处是可能为局部最优。
import warnings
warnings.filterwarnings("ignore")
sm.tsa.arma_order_select_ic(change_data, max_ar=6, max_ma=4, ic="aic")[
"aic_min_order"
] # AIC
(2, 2)
sm.tsa.arma_order_select_ic(change_data, max_ar=6, max_ma=4, ic="bic")[
"bic_min_order"
] # BIC
# (0, 0)
(0, 0)
sm.tsa.arma_order_select_ic(change_data, max_ar=6, max_ma=4, ic="hqic")[
"hqic_min_order"
] # HQIC
# (0, 0)
(0, 0)
可以看出，AIC求解的模型阶次为(5,4)。我们这里就以AIC准则为准~ 至于到底哪种准则更好，可以分别建模进行对比~
ARMA模型的预测和检验
预测的流程和MA流程一致，这里不再重复。预测结果的拟合优度得分仍然很低，但是从预测图上来看，好像比之前强了一点，至少在中间部分效果还不错。
可能的原因有：
信息准则定阶时，为减少计算量，我们设定了AR和MA模型的最大阶数，这样得到的阶数只是局部最优解。
ARMA模型不足以预测股价的变化趋势（这也正常，不然人人都可以赚大钱了）
差分自回归滑动平均模型（Auto-Regressive Integrated Moving Averages，ARIMA）
到目前为止，我们研究的序列都集中在平稳序列。即ARMA模型研究的对象为平稳序列。如果序列是非平稳的，就可以考虑使用ARIMA模型。ARIMA比ARMA仅多了个”I”，代表着其比ARMA多一层内涵：也就是差分。一个非平稳序列经过 d d d次差分后，可以转化为平稳时间序列。 d d d具体的取值，我们得分被对差分1次后的序列进行平稳性检验，若果是非平稳的，则继续差分。直到 d d d次后检验为平稳序列。
平稳时间序列的ARIMA预测的只不过是一个线性方程(如线性回归)。模型有三个主要参数：
Number of AR (Auto-Regressive) terms §： 现在点使用多少个过往数据计算。AR条件仅仅是因变量的滞后。如：如果 p p p等于5，那么预测 x ( t ) x(t) x(t)将是 x ( t − 1 ) . . . x ( t − 5 ) x(t-1)...x(t-5) x(t−1)...x(t−5)。
Number of MA (Moving Average) terms (q)：使用多少个过往的残余错误值。MA条件是预测方程的滞后预测错误。如：如果 q q q等于5，预测 x ( t ) x(t) x(t)将是 e ( t − 1 ) . . . e ( t − 5 ) e(t-1)...e(t-5) e(t−1)...e(t−5)， e ( i ) e(i) e(i)是移动平均数在第 i i i个瞬间和实际值的差值。
Number of Differences (d)：为时间序列成为平稳时所做的差分次数。有非季节性的差值，这种情况下我们采用一阶差分。
ARIMA模型分析框架
下图的框架展示了如何一步一步的“做一个时间序列分析”
前三步我们在前文已经讨论了。尽管如此，这里还是需要简单说明一下：
第一步：时间序列可视化
在构建任何类型的时间序列模型之前，分析其趋势是至关重要的。我们感兴趣的细节与系列中的任何趋势、季节性或随机行为有关。
第二步：时间序列平稳
一旦我们知道了模式、趋势、周期。我们就可以检查序列是否平稳。Dicky-Fuller是一种很流行的检验方式。如果发现序列是非平稳序列怎么办？
这里有三种比较常用的技术来让一个时间序列平稳。
1、消除趋势：这里我们简单的删除时间序列中的趋势成分。例如，我的时间序列的方程是：
x
(
t
)
=
(
m
e
a
n
+
t
r
e
n
d
∗
t
)
+
e
r
r
o
r
x(t)=(mean+trend*t)+error
x(t)=(mean+trend∗t)+error
我们只需删除括号中的部分，并为其余部分建立模型。
2、差分：这是常用的去除非平稳性的技术。这里我们是对序列的差分的结果建立模型而不是真正的序列。例如：
x
(
t
)
–
x
(
t
−
1
)
=
A
R
M
A
(
p
,
q
)
x(t)–x(t−1)=ARMA(p,q)
x(t)–x(t−1)=ARMA(p,q)
这种差异称为 AR(I)MA 中的积分部分。 现在，我们有三个参数：
p:AR
d:I
q:MA
3、季节性：季节性直接被纳入ARIMA模型中。下面的应用部分我们再讨论这个。
第三步：找到最优参数
使用ACF 和PACF 图可以找到参数p、d、q。这种方法的一个补充是，如果ACF 和PACF 逐渐减小，则表明我们需要使时间序列平稳并将值引入到“d”。
第四步：建立ARIMA模型
找到了这些参数，我们现在就可以尝试建立ARIMA模型了。从上一步找到的值可能只是一个近似估计的值，我们需要探索更多(p,d,q）的组合。最小的BIC和AIC的模型参数才是我们要的。我们也可以尝试一些具有季节性成分的模型。在这里，在ACF/PACF图中我们会注意到一些季节性的东西。
第五步：预测
一旦我们有了最后的ARIMA 模型，我们现在就可以预测未来的时间点。如果模型工作正常，我们也可以用可视化趋势来验证。
ARIMA模型参数p、d、q
看完上面的介绍你一定对于如何确定ARIMA模型的参数p、d、q还有疑问。接下来讲解何确定p、d、q参数。
读取并观察数据
import matplotlib.pyplot as plt
import pandas as pd
import statsmodels.api as sm
data = pd.read_csv("arima-demo.csv", parse_dates=["date"], index_col="date")
print(data.head())
data.plot(figsize=(12, 6))
ct
date
2010-01-01 10930
2010-02-01 10318
2010-03-01 10595
2010-04-01 10972
2010-05-01 7706
<AxesSubplot: xlabel='date'>
从上图可知，存在一定的增长趋势。
时间序列的差分d
ARIMA 模型对时间序列的要求是平稳型。因此，当你得到一个非平稳的时间序列时，首先要做的即是做时间序列的差分，直到得到一个平稳时间序列。如果你对时间序列做d次差分才能得到一个平稳序列，那么可以使用ARIMA(p,d,q)模型，其中d是差分次数。
1阶差分：
diff1 = data.diff(1)
diff1.plot(figsize=(12, 6))
<AxesSubplot: xlabel='date'>
目测已经平稳，再来看看2阶差分的效果：
diff2 = data.diff(2)
diff2.plot(figsize=(12, 6))
<AxesSubplot: xlabel='date'>
可以看到二阶差分侯差异不大，所以这里d设置为1即可。
阶数 p 和阶数 q
如何确定差分阶数和常数项：
假如时间序列在很高的lag（10以上）上有正的ACF，则需要很高阶的差分
假如lag-1的ACF是0或负数或者非常小，则不需要很高的差分，假如lag-1的ACF是-0.5或更小，序列可能overdifferenced。
最优的差分阶数一般在最优阶数下标准差最小（但并不是总数如此）
模型不差分意味着原先序列是平稳的；1阶差分意味着原先序列有个固定的平均趋势；二阶差分意味着原先序列有个随时间变化的趋势
模型没有差分一般都有常数项；有两阶差分一般没常数项；假如1阶差分模型非零的平均趋势，则有常数项
如何确定AR和MA的阶数：
假如PACF显示截尾或者lag-1的ACF是正的（此时序列仍然有点underdifferenced），则需要考虑AR项；PACF的截尾项表明AR的阶数
假如ACF显示截尾或者lag-1的ACF是负的（此时序列有点overdifferenced），则需要加MA项，ACF的截尾项表明AR的阶数
AR和MA可以相互抵消对方的影响，所以假如用AR-MA模型去拟合数据，仍然需要考虑加些AR或MA项。尤其在原先模型需要超过10次的迭代去converge。
假如在AR部分有个单位根（AR系数和大约为1），此时应该减少一项AR增加一次差分
假如在MA部分有个单位根（MA系数和大约为1），此时应该减少一项AR减少一次差分
假如长期预测出现不稳定，则可能AR、MA系数有单位根
现在我们已经得到一个平稳的时间序列，接来下就是选择合适的ARIMA模型，即ARIMA模型中合适的p,q。
第一步我们要先检查平稳时间序列的自相关图和偏自相关图。
diff1.dropna(inplace=True)
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(diff1, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(diff1, lags=40, ax=ax2)
其中lags 表示滞后的阶数，以上分别得到acf 图和pacf 图。通过两图观察得到：
自相关图显示滞后有3（4）个阶超出了置信边界
偏相关图显示在滞后1至7阶（lags 1,2,…，7）时的偏自相关系数超出了置信边界，从lag 7之后偏自相关系数值缩小至0
则有以下模型可以供选择：
ARMA(0,1)模型：即自相关图在滞后1阶之后缩小为0，且偏自相关缩小至0，则是一个阶数q=1的移动平均模型；
ARMA(7,0)模型：即偏自相关图在滞后7阶之后缩小为0，且自相关缩小至0，则是一个阶层p=3的自回归模型；
ARMA(7,1)模型：即使得自相关和偏自相关都缩小至零。则是一个混合模型。
为了确定哪个模型最合适，可以采用信息准测AIC、BIC、HQIC判定：
arma_mod70 = sm.tsa.arima.ARIMA(diff1, order=(7, 0, 0)).fit()
print("arma_mod70:", arma_mod70.aic, arma_mod70.bic, arma_mod70.hqic)
arma_mod01 = sm.tsa.arima.ARIMA(diff1, order=(0, 0, 1)).fit()
print("arma_mod01:", arma_mod01.aic, arma_mod01.bic, arma_mod01.hqic)
arma_mod71 = sm.tsa.arima.ARIMA(diff1, order=(7, 0, 1)).fit()
print("arma_mod71:", arma_mod71.aic, arma_mod71.bic, arma_mod71.hqic)
arma_mod70: 1580.221131275376 1602.6188586029652 1589.2490124074332
arma_mod01: 1632.354930566871 1639.8208396760674 1635.3642242775568
arma_mod71: 1581.592994772257 1606.4793584695785 1591.6239738078764
可以看到ARMA(7,0)的aic，bic，hqic均最小，因此是最佳模型。
ARIMA模型的校验
在指数平滑模型下，观察ARIMA模型的残差是否是平均值为0且方差为常数的正态分布（服从零均值、方差不变的正态分布），同时也要观察连续残差是否（自）相关。
残差的自相关与偏自相关
对ARMA(7,0)模型所产生的残差做自相关图：
resid = arma_mod70.resid # 残差
print(resid)
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(resid.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(resid, lags=40, ax=ax2)
date
2010-02-01 -605.990101
2010-03-01 201.791928
2010-04-01 284.361145
2010-05-01 -3221.618824
2010-06-01 -1778.204030
...
2017-02-01 459.916646
2017-03-01 -321.067931
2017-04-01 266.599962
2017-05-01 -1745.634414
2017-06-01 677.536689
Length: 89, dtype: float64
看一看到大部分都在置信空间内，部分超出也只超出一点点。
D－W检验
Durbin-Watson检验，简称D-W检验，是目前检验自相关性最常用的方法，但它只使用于检验一阶自相关性。当DW值显著的接近于O或４时，则存在自相关性，而接近于２时，则不存在（一阶）自相关性。
print(sm.stats.durbin_watson(resid))
2.0349725121809685
观察是否符合正态分布
这里使用QQ图，它用于直观验证一组数据是否来自某个分布，或者验证某两组数据是否来自同一（族）分布。在教学和软件中常用的是检验数据是否来自于正态分布。
from statsmodels.graphics.api import qqplot
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
fig = qqplot(resid, line="q", ax=ax, fit=True)
Ljung-Box检验
Ljung-Box test是对randomness的检验,或者说是对时间序列是否存在滞后相关的一种统计检验。对于滞后相关的检验，我们常常采用的方法还包括计算ACF和PCAF并观察其图像，但是无论是ACF还是PACF都仅仅考虑是否存在某一特定滞后阶数的相关。LB检验则是基于一系列滞后阶数，判断序列总体的相关性或者说随机性是否存在。 时间序列中一个最基本的模型就是高斯白噪声序列。而对于ARIMA模型，其残差被假定为高斯白噪声序列，所以当我们用ARIMA模型去拟合数据时，拟合后我们要对残差的估计序列进行LB检验，判断其是否是高斯白噪声，如果不是，那么就说明ARIMA模型也许并不是一个适合样本的模型。
import numpy as np
r, q, p = sm.tsa.acf(resid.values.squeeze(), qstat=True)
r
array([ 1. , -0.0214771 , -0.05468141, 0.09198236, 0.0423921 ,
0.15310885, -0.01902258, -0.23662027, 0.07314343, -0.01001086,
-0.25444093, -0.08334865, -0.06207171, -0.09326855, 0.18525471,
-0.22255685, 0.01269276, 0.17604155, -0.05585513, -0.058878 ])
q
array([ 0.04245218, 0.32080234, 1.11758907, 1.28882003, 3.54905028,
3.58435982, 9.11431289, 9.64924293, 9.65938867, 16.29648685,
17.01781677, 17.42307146, 18.3500897 , 22.05611738, 27.47714479,
27.49501874, 30.98103391, 31.33691031, 31.73799806])
p
array([0.83676024, 0.851802 , 0.7728298 , 0.8632695 , 0.61597904,
0.73271426, 0.24455275, 0.29051231, 0.37875655, 0.09145386,
0.10734967, 0.13436378, 0.14467344, 0.07746864, 0.02508078,
0.03629829, 0.02007863, 0.02630744, 0.03344411])
data = np.c_[range(1, 20), r[1:], q, p]
table = pd.DataFrame(data, columns=["lag", "AC", "Q", "Prob(>Q)"])
print(table.set_index("lag"))
AC Q Prob(>Q)
lag
1.0 -0.021477 0.042452 0.836760
2.0 -0.054681 0.320802 0.851802
3.0 0.091982 1.117589 0.772830
4.0 0.042392 1.288820 0.863270
5.0 0.153109 3.549050 0.615979
6.0 -0.019023 3.584360 0.732714
7.0 -0.236620 9.114313 0.244553
8.0 0.073143 9.649243 0.290512
9.0 -0.010011 9.659389 0.378757
10.0 -0.254441 16.296487 0.091454
11.0 -0.083349 17.017817 0.107350
12.0 -0.062072 17.423071 0.134364
13.0 -0.093269 18.350090 0.144673
14.0 0.185255 22.056117 0.077469
15.0 -0.222557 27.477145 0.025081
16.0 0.012693 27.495019 0.036298
17.0 0.176042 30.981034 0.020079
18.0 -0.055855 31.336910 0.026307
19.0 -0.058878 31.737998 0.033444
检验的结果就是看最后一列前十二行的检验概率（一般观察滞后1~12阶），如果检验概率小于给定的显著性水平，比如0.05就拒绝原假设，其原假设是相关系数为零。就结果来看，如果取显著性水平大于0.05，那么相关系数与零没有显著差异，即为白噪声序列。
ARIMA模型实战，通过ACF和PCF确定参数
测试数据使用上面平稳序列使用的AirPassengers数据。
绘制ACF和PACF图
时间序列的自回归函数和偏自回归函数可以在差分后绘制为：
import matplotlib.pylab as plt
import numpy as np
import pandas as pd
data = pd.read_csv(
"AirPassengers.csv", parse_dates=["TravelDate"], index_col="TravelDate"
)
ts = data["Passengers"]
ts_log = np.log(ts)
ts_log_diff = ts_log.diff()
ts_log_diff.dropna(inplace=True)
print(ts_log_diff)
TravelDate
1949-02-01 0.052186
1949-03-01 0.112117
1949-04-01 -0.022990
1949-05-01 -0.064022
1949-06-01 0.109484
...
1960-08-01 -0.026060
1960-09-01 -0.176399
1960-10-01 -0.097083
1960-11-01 -0.167251
1960-12-01 0.102279
Name: Passengers, Length: 143, dtype: float64
# ACF and PACF plots:
from statsmodels.tsa.stattools import acf, pacf
lag_acf = acf(ts_log_diff, nlags=20)
lag_pacf = pacf(ts_log_diff, nlags=20, method="ols")
# Plot ACF:
plt.subplot(121)
plt.plot(lag_acf)
plt.axhline(y=0, linestyle="--", color="gray")
plt.axhline(y=-1.96 / np.sqrt(len(ts_log_diff)), linestyle="--", color="gray")
plt.axhline(y=1.96 / np.sqrt(len(ts_log_diff)), linestyle="--", color="gray")
plt.title("Autocorrelation Function")
# Plot PACF:
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0, linestyle="--", color="gray")
plt.axhline(y=-1.96 / np.sqrt(len(ts_log_diff)), linestyle="--", color="gray")
plt.axhline(y=1.96 / np.sqrt(len(ts_log_diff)), linestyle="--", color="gray")
plt.title("Partial Autocorrelation Function")
plt.tight_layout()
在这个例子中，0刻度线上线的2条虚线为置信区间，用来确定“p”和“q”的值。
p：部分自相关函数表第一次截断的上层置信区间是滞后值。如果你仔细看，该值是p=0
q：自相关函数表第一次截断的上层置信区间是滞后值。如果你仔细看，该值是q=1
现在开始搭建ARIMA的三个模型，及计算各自的RSS（这里的RSS是指残差值,而不是实际序列）
自回归（AR）模型
model = sm.tsa.arima.ARIMA(ts_log, order=(0, 1, 0))
results_AR = model.fit()
plt.plot(ts_log_diff)
plt.plot(results_AR.fittedvalues, color="red")
plt.title("RSS: %.4f" % sum((results_AR.fittedvalues[1:] - ts_log_diff) ** 2))
Text(0.5, 1.0, 'RSS: 4402.7940')
移动平均数（MA）模型
model = sm.tsa.arima.ARIMA(ts_log, order=(0, 1, 1))
results_MA = model.fit()
plt.plot(ts_log_diff)
plt.plot(results_MA.fittedvalues, color="red")
plt.title("RSS: %.4f" % sum((results_MA.fittedvalues[1:] - ts_log_diff) ** 2))
Text(0.5, 1.0, 'RSS: 4405.7765')
ARMA组合
在本例中，由于AR的p值为0，所以组合的结果与MA一致：
model = sm.tsa.arima.ARIMA(ts_log, order=(3, 1, 4))
results_ARIMA = model.fit()
plt.plot(ts_log_diff)
plt.plot(results_ARIMA.fittedvalues, color="red")
plt.title("RSS: %.4f" % sum((results_ARIMA.fittedvalues[1:] - ts_log_diff) ** 2))
Text(0.5, 1.0, 'RSS: 4379.8297')
转化为原数据空间
前面模型采用的数据都是转化后的数据，为了得到最终的结果，需要将数据转化为原数据空间。第一步将预测结果保存为序列。
predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True)
print(predictions_ARIMA_diff.head())
TravelDate
1949-01-01 0.000000
1949-02-01 4.718499
1949-03-01 4.783835
1949-04-01 4.904676
1949-05-01 4.820349
dtype: float64
注意这些是从‘1949-02-01’开始而不是第一个月。为什么?这是因为我们将第一个月份取为滞后值，一月前面没有可以减去的元素。一个简单的解决方法是先确定索引的累积和，然后将它添加到基数。累积和的计算方式如下：
predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum()
print(predictions_ARIMA_diff_cumsum.head())
TravelDate
1949-01-01 0.000000
1949-02-01 4.718499
1949-03-01 9.502334
1949-04-01 14.407010
1949-05-01 19.227359
dtype: float64
将差分转换为对数尺度的方法是这些差值连续地添加到基数。（第一个元素是基数本身）
predictions_ARIMA_log = pd.Series(ts_log.iloc[0], index=ts_log.index)
predictions_ARIMA_log = predictions_ARIMA_log.add(
predictions_ARIMA_diff_cumsum, fill_value=0
)
predictions_ARIMA_log
TravelDate
1949-01-01 4.718499
1949-02-01 9.436998
1949-03-01 14.220833
1949-04-01 19.125509
1949-05-01 23.945858
...
1960-08-01 770.316474
1960-09-01 776.592811
1960-10-01 782.703510
1960-11-01 788.803931
1960-12-01 794.800786
Length: 144, dtype: float64
最后一步是将指数与原序列比较。
predictions_ARIMA = np.exp(predictions_ARIMA_log)
plt.plot(ts)
plt.plot(predictions_ARIMA)
plt.title("RMSE: %.4f" % np.sqrt(sum((predictions_ARIMA - ts) ** 2) / len(ts)))
Text(0.5, 1.0, 'RMSE: inf')
predictions_ARIMA
TravelDate
1949-01-01 1.120000e+02
1949-02-01 1.254400e+04
1949-03-01 1.499786e+06
1949-04-01 2.023499e+08
1949-05-01 2.509308e+10
...
1960-08-01 inf
1960-09-01 inf
1960-10-01 inf
1960-11-01 inf
1960-12-01 inf
Length: 144, dtype: float64
ARIMA模型实战，基于“分解”进行未来预测
前面讲到了，分解出来的df测试值要优于移动平均和移动加权平均，这里就使用“分解”后的数据进行进一步的预测。
对残差值进行预测
residual.dropna(inplace=True)
model_ARIMA = sm.tsa.arima.ARIMA(residual, order=(0, 1, 1)).fit()
predictions_ARIMA = model_ARIMA.predict(start="1950-01-01", end="1962-04-01")
plt.plot(residual)
plt.plot(predictions_ARIMA)
[<matplotlib.lines.Line2D at 0x186e09c6190>]
这里的(p,d,q)与就使用上一篇文章中的(0,1,1)。预测后的数据绘制出来为：
在残差预测的基础上加上趋势和季节数据
predictions_ARIMA = predictions_ARIMA.add(trend, fill_value=0).add(
seasonal, fill_value=0
)
predictions_ARIMA = np.exp(predictions_ARIMA)
plt.plot(data["Passengers"], color="blue")
plt.plot(predictions_ARIMA, color="red")
[<matplotlib.lines.Line2D at 0x186e0a886d0>]
可以看到在中间匹配的还可以，而开头和最后预测结果比较异常，通过查看trend、seasonal的数据可发现其存在数据缺失。trend在1960-06之后便没有数据了，seasonal在1960-12之后也结束了，所以会导致最终的数据在1960-06之后便为0了。
print(trend)
print(seasonal)
TravelDate
1949-01-01 NaN
1949-02-01 NaN
1949-03-01 NaN
1949-04-01 NaN
1949-05-01 NaN
..
1960-08-01 NaN
1960-09-01 NaN
1960-10-01 NaN
1960-11-01 NaN
1960-12-01 NaN
Name: trend, Length: 144, dtype: float64
TravelDate
1949-01-01 -0.085815
1949-02-01 -0.114413
1949-03-01 0.018113
1949-04-01 -0.013046
1949-05-01 -0.008966
...
1960-08-01 0.204512
1960-09-01 0.064836
1960-10-01 -0.075271
1960-11-01 -0.215846
1960-12-01 -0.100315
Name: seasonal, Length: 144, dtype: float64
对trend和seasonal进行预测
我们看看trend的形状：
plt.plot(trend)
[<matplotlib.lines.Line2D at 0x186e0af8bb0>]
上面趋势近似线性，我们可以用线性拟合来预测1960-06之后的数据：
trend.dropna(inplace=True)
from sklearn.linear_model import LinearRegression
x = pd.Series(range(trend.size), index=trend.index)
x = x.to_frame()
linreg = LinearRegression()
linereg = linreg.fit(x, trend)
x = pd.Series(range(0, 154), index=(pd.period_range("1949-07", periods=154, freq="M")))
x = x.to_frame()
res_predict = linereg.predict(x)
trend2 = pd.Series(res_predict, index=x.index).to_timestamp()
plt.plot(trend, color="blue")
plt.plot(trend2, color="red")
[<matplotlib.lines.Line2D at 0x186e1360070>]
由于seasonal是周期性的，直接用shift即可：
index1 = pd.date_range(start="1949-01-01", end="1962-04-01", freq="MS")
seasonal = seasonal.reindex(index1)
seasonal = seasonal.shift(24)
plt.plot(seasonal)
[<matplotlib.lines.Line2D at 0x186e1375520>]
在预测残差的基础上，重新添加trend和seasonal
model_ARIMA = sm.tsa.arima.ARIMA(residual, order=(0, 1, 1)).fit()
predictions_ARIMA = model_ARIMA.predict(start="1950-01-01", end="1962-04-01")
predictions_ARIMA = predictions_ARIMA.add(trend2, fill_value=0).add(
seasonal, fill_value=0
)
predictions_ARIMA = np.exp(predictions_ARIMA)
plt.plot(data["Passengers"], color="blue")
plt.plot(predictions_ARIMA, color="red")
[<matplotlib.lines.Line2D at 0x186e141e400>]
ARIMA模型实战，基于网格搜索确定参数
确定ARIMA参数值
从上面的ACF和PACF可大致确定ARIMA的参数，p=偏自相关截尾的lag。如果偏自相关截尾后，自相关仍然拖尾，考虑增加p；q=自相关截尾的lag。如果自相关截尾后，偏自相关仍然拖尾，考虑增加q。以上仅是理论上的说法，具体最合适的p,q，还是要以AIC、BIC或预测结果残差最小来循环确定。ARIMA模型是通过优化AIC, BIC, HQIC等得出的, 一般来说AIC较小的模型拟合度较好,但拟合度较好不能说明预测能力好。
ts_week_log = ts_log
ts_week_log_diff = ts_log_diff
d = 1 # 一阶差分后平稳，所以d=1。p,q参数使用循环暴力调参确定最佳值
pdq_result = pd.DataFrame()
for p in range(0, 5):
for q in range(0, 5):
model = sm.tsa.arima.ARIMA(
ts_week_log, order=(p, d, q)
) # ARIMA的参数: order = p, d, q
results_ARIMA = model.fit()
pdq_result = pdq_result.append(
pd.DataFrame(
{
"p": [p],
"d": [d],
"q": [q],
"AIC": [results_ARIMA.aic],
"BIC": [results_ARIMA.bic],
"HQIC": [results_ARIMA.hqic],
}
)
)
print(results_ARIMA.summary()) # summary2()也可以使用
print(pdq_result)
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(0, 1, 0) Log Likelihood 117.220
Date: Thu, 05 Jan 2023 AIC -232.439
Time: 15:34:38 BIC -229.476
Sample: 01-01-1949 HQIC -231.235
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
sigma2 0.0114 0.002 5.778 0.000 0.008 0.015
===================================================================================
Ljung-Box (L1) (Q): 5.83 Jarque-Bera (JB): 6.78
Prob(Q): 0.02 Prob(JB): 0.03
Heteroskedasticity (H): 1.23 Skew: -0.09
Prob(H) (two-sided): 0.48 Kurtosis: 1.95
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(0, 1, 1) Log Likelihood 121.363
Date: Thu, 05 Jan 2023 AIC -238.725
Time: 15:34:38 BIC -232.800
Sample: 01-01-1949 HQIC -236.317
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ma.L1 0.2768 0.104 2.668 0.008 0.073 0.480
sigma2 0.0107 0.002 6.286 0.000 0.007 0.014
===================================================================================
Ljung-Box (L1) (Q): 0.15 Jarque-Bera (JB): 6.45
Prob(Q): 0.69 Prob(JB): 0.04
Heteroskedasticity (H): 1.16 Skew: 0.22
Prob(H) (two-sided): 0.61 Kurtosis: 2.05
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(0, 1, 2) Log Likelihood 123.423
Date: Thu, 05 Jan 2023 AIC -240.847
Time: 15:34:38 BIC -231.958
Sample: 01-01-1949 HQIC -237.235
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ma.L1 0.2343 0.085 2.757 0.006 0.068 0.401
ma.L2 -0.2831 0.078 -3.613 0.000 -0.437 -0.130
sigma2 0.0104 0.002 5.644 0.000 0.007 0.014
===================================================================================
Ljung-Box (L1) (Q): 0.01 Jarque-Bera (JB): 7.34
Prob(Q): 0.93 Prob(JB): 0.03
Heteroskedasticity (H): 1.17 Skew: -0.07
Prob(H) (two-sided): 0.60 Kurtosis: 1.90
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(0, 1, 3) Log Likelihood 124.765
Date: Thu, 05 Jan 2023 AIC -241.530
Time: 15:34:38 BIC -229.679
Sample: 01-01-1949 HQIC -236.714
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ma.L1 0.1132 0.091 1.238 0.216 -0.066 0.293
ma.L2 -0.4618 0.078 -5.886 0.000 -0.616 -0.308
ma.L3 -0.1759 0.088 -1.992 0.046 -0.349 -0.003
sigma2 0.0102 0.002 5.610 0.000 0.007 0.014
===================================================================================
Ljung-Box (L1) (Q): 0.08 Jarque-Bera (JB): 7.39
Prob(Q): 0.77 Prob(JB): 0.02
Heteroskedasticity (H): 1.24 Skew: -0.12
Prob(H) (two-sided): 0.46 Kurtosis: 1.91
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(0, 1, 4) Log Likelihood 135.145
Date: Thu, 05 Jan 2023 AIC -260.290
Time: 15:34:39 BIC -245.475
Sample: 01-01-1949 HQIC -254.270
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ma.L1 0.1898 0.096 1.986 0.047 0.002 0.377
ma.L2 -0.1798 0.104 -1.728 0.084 -0.384 0.024
ma.L3 -0.1395 0.089 -1.563 0.118 -0.314 0.035
ma.L4 -0.4043 0.109 -3.706 0.000 -0.618 -0.190
sigma2 0.0088 0.001 6.139 0.000 0.006 0.012
===================================================================================
Ljung-Box (L1) (Q): 0.45 Jarque-Bera (JB): 7.25
Prob(Q): 0.50 Prob(JB): 0.03
Heteroskedasticity (H): 1.03 Skew: 0.29
Prob(H) (two-sided): 0.91 Kurtosis: 2.06
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(1, 1, 0) Log Likelihood 120.299
Date: Thu, 05 Jan 2023 AIC -236.598
Time: 15:34:39 BIC -230.672
Sample: 01-01-1949 HQIC -234.190
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.2054 0.093 2.204 0.028 0.023 0.388
sigma2 0.0109 0.002 6.050 0.000 0.007 0.014
===================================================================================
Ljung-Box (L1) (Q): 0.13 Jarque-Bera (JB): 6.19
Prob(Q): 0.71 Prob(JB): 0.05
Heteroskedasticity (H): 1.16 Skew: 0.19
Prob(H) (two-sided): 0.60 Kurtosis: 2.05
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(1, 1, 1) Log Likelihood 124.313
Date: Thu, 05 Jan 2023 AIC -242.626
Time: 15:34:39 BIC -233.738
Sample: 01-01-1949 HQIC -239.014
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 -0.5773 0.164 -3.516 0.000 -0.899 -0.255
ma.L1 0.8478 0.098 8.684 0.000 0.656 1.039
sigma2 0.0103 0.002 5.992 0.000 0.007 0.014
===================================================================================
Ljung-Box (L1) (Q): 0.02 Jarque-Bera (JB): 5.94
Prob(Q): 0.90 Prob(JB): 0.05
Heteroskedasticity (H): 1.07 Skew: 0.04
Prob(H) (two-sided): 0.82 Kurtosis: 2.00
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(1, 1, 2) Log Likelihood 128.432
Date: Thu, 05 Jan 2023 AIC -248.865
Time: 15:34:39 BIC -237.013
Sample: 01-01-1949 HQIC -244.049
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.5432 0.133 4.069 0.000 0.282 0.805
ma.L1 -0.3822 0.148 -2.586 0.010 -0.672 -0.092
ma.L2 -0.4096 0.096 -4.266 0.000 -0.598 -0.221
sigma2 0.0097 0.002 5.646 0.000 0.006 0.013
===================================================================================
Ljung-Box (L1) (Q): 0.28 Jarque-Bera (JB): 8.01
Prob(Q): 0.60 Prob(JB): 0.02
Heteroskedasticity (H): 1.16 Skew: 0.03
Prob(H) (two-sided): 0.61 Kurtosis: 1.84
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(1, 1, 3) Log Likelihood 128.433
Date: Thu, 05 Jan 2023 AIC -246.866
Time: 15:34:39 BIC -232.052
Sample: 01-01-1949 HQIC -240.846
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.5402 0.325 1.664 0.096 -0.096 1.177
ma.L1 -0.3799 0.367 -1.034 0.301 -1.100 0.340
ma.L2 -0.4063 0.097 -4.190 0.000 -0.596 -0.216
ma.L3 -0.0044 0.187 -0.024 0.981 -0.370 0.362
sigma2 0.0097 0.002 5.647 0.000 0.006 0.013
===================================================================================
Ljung-Box (L1) (Q): 0.28 Jarque-Bera (JB): 8.01
Prob(Q): 0.59 Prob(JB): 0.02
Heteroskedasticity (H): 1.16 Skew: 0.04
Prob(H) (two-sided): 0.61 Kurtosis: 1.84
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(1, 1, 4) Log Likelihood 135.160
Date: Thu, 05 Jan 2023 AIC -258.320
Time: 15:34:40 BIC -240.543
Sample: 01-01-1949 HQIC -251.097
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 -0.0395 0.199 -0.199 0.842 -0.429 0.350
ma.L1 0.2245 0.185 1.212 0.226 -0.139 0.588
ma.L2 -0.1662 0.112 -1.479 0.139 -0.386 0.054
ma.L3 -0.1544 0.096 -1.605 0.108 -0.343 0.034
ma.L4 -0.4197 0.111 -3.770 0.000 -0.638 -0.201
sigma2 0.0088 0.001 6.189 0.000 0.006 0.012
===================================================================================
Ljung-Box (L1) (Q): 0.40 Jarque-Bera (JB): 7.10
Prob(Q): 0.53 Prob(JB): 0.03
Heteroskedasticity (H): 1.04 Skew: 0.30
Prob(H) (two-sided): 0.90 Kurtosis: 2.08
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(2, 1, 0) Log Likelihood 122.261
Date: Thu, 05 Jan 2023 AIC -238.522
Time: 15:34:40 BIC -229.633
Sample: 01-01-1949 HQIC -234.910
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.2414 0.101 2.381 0.017 0.043 0.440
ar.L2 -0.1667 0.100 -1.665 0.096 -0.363 0.030
sigma2 0.0106 0.002 6.020 0.000 0.007 0.014
===================================================================================
Ljung-Box (L1) (Q): 0.06 Jarque-Bera (JB): 6.58
Prob(Q): 0.81 Prob(JB): 0.04
Heteroskedasticity (H): 1.20 Skew: 0.14
Prob(H) (two-sided): 0.53 Kurtosis: 1.98
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(2, 1, 1) Log Likelihood 129.732
Date: Thu, 05 Jan 2023 AIC -251.464
Time: 15:34:40 BIC -239.612
Sample: 01-01-1949 HQIC -246.648
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.9792 0.100 9.839 0.000 0.784 1.174
ar.L2 -0.3740 0.118 -3.177 0.001 -0.605 -0.143
ma.L1 -0.8303 0.102 -8.118 0.000 -1.031 -0.630
sigma2 0.0095 0.002 6.134 0.000 0.006 0.013
===================================================================================
Ljung-Box (L1) (Q): 1.02 Jarque-Bera (JB): 8.63
Prob(Q): 0.31 Prob(JB): 0.01
Heteroskedasticity (H): 1.14 Skew: 0.22
Prob(H) (two-sided): 0.66 Kurtosis: 1.88
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(2, 1, 2) Log Likelihood 128.890
Date: Thu, 05 Jan 2023 AIC -247.780
Time: 15:34:41 BIC -232.965
Sample: 01-01-1949 HQIC -241.760
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.2432 0.166 1.469 0.142 -0.081 0.568
ar.L2 0.2677 0.199 1.344 0.179 -0.123 0.658
ma.L1 -0.0867 0.118 -0.732 0.464 -0.319 0.145
ma.L2 -0.6813 0.142 -4.786 0.000 -0.960 -0.402
sigma2 0.0096 0.002 5.730 0.000 0.006 0.013
===================================================================================
Ljung-Box (L1) (Q): 0.01 Jarque-Bera (JB): 7.04
Prob(Q): 0.94 Prob(JB): 0.03
Heteroskedasticity (H): 1.08 Skew: -0.06
Prob(H) (two-sided): 0.79 Kurtosis: 1.92
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(2, 1, 3) Log Likelihood 130.066
Date: Thu, 05 Jan 2023 AIC -248.132
Time: 15:34:41 BIC -230.355
Sample: 01-01-1949 HQIC -240.908
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 -0.1264 0.452 -0.279 0.780 -1.013 0.760
ar.L2 0.4325 0.158 2.742 0.006 0.123 0.742
ma.L1 0.2982 0.477 0.626 0.532 -0.636 1.233
ma.L2 -0.7567 0.118 -6.399 0.000 -0.989 -0.525
ma.L3 -0.2157 0.271 -0.797 0.426 -0.746 0.315
sigma2 0.0095 0.002 5.763 0.000 0.006 0.013
===================================================================================
Ljung-Box (L1) (Q): 0.06 Jarque-Bera (JB): 7.12
Prob(Q): 0.81 Prob(JB): 0.03
Heteroskedasticity (H): 1.07 Skew: -0.01
Prob(H) (two-sided): 0.82 Kurtosis: 1.91
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(2, 1, 4) Log Likelihood 141.432
Date: Thu, 05 Jan 2023 AIC -268.864
Time: 15:34:41 BIC -248.124
Sample: 01-01-1949 HQIC -260.436
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 -0.2815 0.120 -2.343 0.019 -0.517 -0.046
ar.L2 -0.3801 0.133 -2.865 0.004 -0.640 -0.120
ma.L1 0.4584 0.089 5.179 0.000 0.285 0.632
ma.L2 0.2640 0.110 2.401 0.016 0.048 0.480
ma.L3 -0.2113 0.089 -2.362 0.018 -0.387 -0.036
ma.L4 -0.6514 0.088 -7.415 0.000 -0.824 -0.479
sigma2 0.0078 0.001 7.653 0.000 0.006 0.010
===================================================================================
Ljung-Box (L1) (Q): 0.00 Jarque-Bera (JB): 1.98
Prob(Q): 0.97 Prob(JB): 0.37
Heteroskedasticity (H): 1.05 Skew: 0.19
Prob(H) (two-sided): 0.86 Kurtosis: 2.57
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(3, 1, 0) Log Likelihood 122.822
Date: Thu, 05 Jan 2023 AIC -237.644
Time: 15:34:42 BIC -225.793
Sample: 01-01-1949 HQIC -232.828
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.2269 0.103 2.200 0.028 0.025 0.429
ar.L2 -0.1444 0.109 -1.320 0.187 -0.359 0.070
ar.L3 -0.0892 0.084 -1.065 0.287 -0.254 0.075
sigma2 0.0105 0.002 6.007 0.000 0.007 0.014
===================================================================================
Ljung-Box (L1) (Q): 0.17 Jarque-Bera (JB): 6.87
Prob(Q): 0.68 Prob(JB): 0.03
Heteroskedasticity (H): 1.21 Skew: 0.24
Prob(H) (two-sided): 0.51 Kurtosis: 2.04
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(3, 1, 1) Log Likelihood 130.225
Date: Thu, 05 Jan 2023 AIC -250.450
Time: 15:34:42 BIC -235.636
Sample: 01-01-1949 HQIC -244.430
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.9318 0.130 7.181 0.000 0.677 1.186
ar.L2 -0.2911 0.157 -1.851 0.064 -0.599 0.017
ar.L3 -0.0888 0.115 -0.770 0.441 -0.315 0.137
ma.L1 -0.8124 0.124 -6.550 0.000 -1.056 -0.569
sigma2 0.0094 0.002 6.263 0.000 0.006 0.012
===================================================================================
Ljung-Box (L1) (Q): 0.78 Jarque-Bera (JB): 8.97
Prob(Q): 0.38 Prob(JB): 0.01
Heteroskedasticity (H): 1.12 Skew: 0.28
Prob(H) (two-sided): 0.70 Kurtosis: 1.90
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(3, 1, 2) Log Likelihood 131.900
Date: Thu, 05 Jan 2023 AIC -251.800
Time: 15:34:42 BIC -234.023
Sample: 01-01-1949 HQIC -244.576
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.1454 0.168 0.867 0.386 -0.183 0.474
ar.L2 0.4038 0.166 2.430 0.015 0.078 0.729
ar.L3 -0.2409 0.138 -1.751 0.080 -0.511 0.029
ma.L1 0.0634 0.152 0.418 0.676 -0.234 0.360
ma.L2 -0.7558 0.118 -6.391 0.000 -0.988 -0.524
sigma2 0.0092 0.002 5.920 0.000 0.006 0.012
===================================================================================
Ljung-Box (L1) (Q): 1.34 Jarque-Bera (JB): 7.80
Prob(Q): 0.25 Prob(JB): 0.02
Heteroskedasticity (H): 1.05 Skew: 0.11
Prob(H) (two-sided): 0.87 Kurtosis: 1.88
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(3, 1, 3) Log Likelihood 145.586
Date: Thu, 05 Jan 2023 AIC -277.172
Time: 15:34:43 BIC -256.432
Sample: 01-01-1949 HQIC -268.745
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.7297 0.077 9.424 0.000 0.578 0.881
ar.L2 0.6472 0.122 5.302 0.000 0.408 0.886
ar.L3 -0.8897 0.078 -11.469 0.000 -1.042 -0.738
ma.L1 -0.8285 0.300 -2.765 0.006 -1.416 -0.241
ma.L2 -0.8462 0.523 -1.617 0.106 -1.872 0.179
ma.L3 0.9782 0.272 3.593 0.000 0.445 1.512
sigma2 0.0073 0.002 3.463 0.001 0.003 0.011
===================================================================================
Ljung-Box (L1) (Q): 0.50 Jarque-Bera (JB): 7.28
Prob(Q): 0.48 Prob(JB): 0.03
Heteroskedasticity (H): 1.02 Skew: 0.27
Prob(H) (two-sided): 0.95 Kurtosis: 2.03
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(3, 1, 4) Log Likelihood 143.428
Date: Thu, 05 Jan 2023 AIC -270.855
Time: 15:34:43 BIC -247.153
Sample: 01-01-1949 HQIC -261.224
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.9615 0.177 5.447 0.000 0.616 1.308
ar.L2 -0.8704 0.178 -4.888 0.000 -1.219 -0.521
ar.L3 0.0482 0.172 0.280 0.780 -0.290 0.386
ma.L1 -0.8711 0.165 -5.270 0.000 -1.195 -0.547
ma.L2 0.6169 0.187 3.299 0.001 0.250 0.983
ma.L3 0.0966 0.189 0.512 0.609 -0.273 0.466
ma.L4 -0.4669 0.094 -4.962 0.000 -0.651 -0.282
sigma2 0.0074 0.001 5.803 0.000 0.005 0.010
===================================================================================
Ljung-Box (L1) (Q): 0.27 Jarque-Bera (JB): 6.36
Prob(Q): 0.60 Prob(JB): 0.04
Heteroskedasticity (H): 1.10 Skew: -0.02
Prob(H) (two-sided): 0.74 Kurtosis: 1.97
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(4, 1, 0) Log Likelihood 129.987
Date: Thu, 05 Jan 2023 AIC -249.974
Time: 15:34:44 BIC -235.159
Sample: 01-01-1949 HQIC -243.954
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.2037 0.092 2.219 0.026 0.024 0.384
ar.L2 -0.1935 0.099 -1.947 0.052 -0.388 0.001
ar.L3 -0.0181 0.079 -0.228 0.820 -0.174 0.137
ar.L4 -0.3136 0.078 -4.007 0.000 -0.467 -0.160
sigma2 0.0095 0.002 5.758 0.000 0.006 0.013
===================================================================================
Ljung-Box (L1) (Q): 0.01 Jarque-Bera (JB): 8.95
Prob(Q): 0.92 Prob(JB): 0.01
Heteroskedasticity (H): 1.18 Skew: 0.20
Prob(H) (two-sided): 0.57 Kurtosis: 1.84
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(4, 1, 1) Log Likelihood 130.063
Date: Thu, 05 Jan 2023 AIC -248.127
Time: 15:34:44 BIC -230.349
Sample: 01-01-1949 HQIC -240.903
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.0980 0.348 0.282 0.778 -0.584 0.780
ar.L2 -0.1691 0.147 -1.152 0.249 -0.457 0.119
ar.L3 -0.0337 0.091 -0.372 0.710 -0.211 0.144
ar.L4 -0.3227 0.085 -3.817 0.000 -0.488 -0.157
ma.L1 0.1172 0.383 0.306 0.760 -0.633 0.868
sigma2 0.0095 0.002 5.613 0.000 0.006 0.013
===================================================================================
Ljung-Box (L1) (Q): 0.04 Jarque-Bera (JB): 8.79
Prob(Q): 0.84 Prob(JB): 0.01
Heteroskedasticity (H): 1.18 Skew: 0.22
Prob(H) (two-sided): 0.57 Kurtosis: 1.86
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(4, 1, 2) Log Likelihood 135.218
Date: Thu, 05 Jan 2023 AIC -256.435
Time: 15:34:44 BIC -235.696
Sample: 01-01-1949 HQIC -248.008
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.1625 0.156 1.043 0.297 -0.143 0.468
ar.L2 0.3247 0.180 1.800 0.072 -0.029 0.678
ar.L3 -0.1475 0.107 -1.376 0.169 -0.358 0.063
ar.L4 -0.2461 0.109 -2.258 0.024 -0.460 -0.032
ma.L1 7.278e-06 0.140 5.2e-05 1.000 -0.274 0.274
ma.L2 -0.6176 0.182 -3.385 0.001 -0.975 -0.260
sigma2 0.0088 0.001 6.086 0.000 0.006 0.012
===================================================================================
Ljung-Box (L1) (Q): 0.62 Jarque-Bera (JB): 8.05
Prob(Q): 0.43 Prob(JB): 0.02
Heteroskedasticity (H): 0.96 Skew: 0.30
Prob(H) (two-sided): 0.88 Kurtosis: 2.01
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(4, 1, 3) Log Likelihood 149.371
Date: Thu, 05 Jan 2023 AIC -282.743
Time: 15:34:45 BIC -259.040
Sample: 01-01-1949 HQIC -273.111
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.2149 0.127 1.691 0.091 -0.034 0.464
ar.L2 -0.3598 0.105 -3.438 0.001 -0.565 -0.155
ar.L3 0.6298 0.100 6.307 0.000 0.434 0.826
ar.L4 -0.5320 0.122 -4.360 0.000 -0.771 -0.293
ma.L1 0.0333 0.109 0.306 0.760 -0.180 0.247
ma.L2 0.2540 0.107 2.366 0.018 0.044 0.464
ma.L3 -0.8206 0.150 -5.456 0.000 -1.115 -0.526
sigma2 0.0070 0.001 5.326 0.000 0.004 0.010
===================================================================================
Ljung-Box (L1) (Q): 0.79 Jarque-Bera (JB): 1.09
Prob(Q): 0.37 Prob(JB): 0.58
Heteroskedasticity (H): 1.10 Skew: -0.21
Prob(H) (two-sided): 0.74 Kurtosis: 2.93
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
SARIMAX Results
==============================================================================
Dep. Variable: Passengers No. Observations: 144
Model: ARIMA(4, 1, 4) Log Likelihood 155.529
Date: Thu, 05 Jan 2023 AIC -293.057
Time: 15:34:45 BIC -266.392
Sample: 01-01-1949 HQIC -282.222
- 12-01-1960
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.6332 0.044 14.291 0.000 0.546 0.720
ar.L2 -0.1149 0.059 -1.955 0.051 -0.230 0.000
ar.L3 0.5871 0.050 11.720 0.000 0.489 0.685
ar.L4 -0.8900 0.044 -20.382 0.000 -0.976 -0.804
ma.L1 -0.6381 0.228 -2.805 0.005 -1.084 -0.192
ma.L2 -0.2255 0.049 -4.618 0.000 -0.321 -0.130
ma.L3 -0.6134 0.124 -4.934 0.000 -0.857 -0.370
ma.L4 0.9806 0.249 3.940 0.000 0.493 1.468
sigma2 0.0062 0.002 3.708 0.000 0.003 0.009
===================================================================================
Ljung-Box (L1) (Q): 0.92 Jarque-Bera (JB): 5.73
Prob(Q): 0.34 Prob(JB): 0.06
Heteroskedasticity (H): 1.07 Skew: -0.07
Prob(H) (two-sided): 0.82 Kurtosis: 2.03
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
p d q AIC BIC HQIC
0 0 1 0 -232.439087 -229.476242 -231.235129
0 0 1 1 -238.725393 -232.799704 -236.317477
0 0 1 2 -240.846549 -231.958015 -237.234675
0 0 1 3 -241.530322 -229.678943 -236.714489
0 0 1 4 -260.289581 -245.475357 -254.269790
0 1 1 0 -236.597873 -230.672184 -234.189957
0 1 1 1 -242.626205 -233.737671 -239.014330
0 1 1 2 -248.864828 -237.013449 -244.048995
0 1 1 3 -246.866208 -232.051985 -240.846417
0 1 1 4 -258.320336 -240.543268 -251.096587
0 2 1 0 -238.521899 -229.633365 -234.910025
0 2 1 1 -251.463509 -239.612131 -246.647676
0 2 1 2 -247.779684 -232.965461 -241.759893
0 2 1 3 -248.131739 -230.354671 -240.907990
0 2 1 4 -268.863698 -248.123786 -260.435991
0 3 1 0 -237.644076 -225.792697 -232.828243
0 3 1 1 -250.449947 -235.635724 -244.430156
0 3 1 2 -251.799916 -234.022849 -244.576167
0 3 1 3 -277.172274 -256.432362 -268.744567
0 3 1 4 -270.855334 -247.152577 -261.223669
0 4 1 0 -249.973568 -235.159344 -243.953777
0 4 1 1 -248.126543 -230.349475 -240.902794
0 4 1 2 -256.435421 -235.695509 -248.007714
0 4 1 3 -282.742755 -259.039998 -273.111090
0 4 1 4 -293.057293 -266.391692 -282.221670
model = sm.tsa.arima.ARIMA(ts_week_log, order=(4, 1, 4)) # ARIMA的参数: order = p, d, q
results_ARIMA = model.fit()
plt.plot(ts_week_log_diff)
plt.plot(
results_ARIMA.fittedvalues, color="red"
) # fittedvalues返回的是d次差分后的序列, 因此和ts_week_log_diff比较
residuals = pd.DataFrame(results_ARIMA.resid) # 拟合的数据和原始数据ts_week_log的残差
plt.show()
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(residuals.values.squeeze(), lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(residuals, lags=40, ax=ax2)
因为ARIMA模型会使拟合和原始数据的残差变为白噪声，所以残差的核密度(概率密度)为正态分布：
print("Residuals Summary:\n", residuals.describe())
residuals.plot(kind="kde") # 概率密度图kde
Residuals Summary:
0
count 144.000000
mean 0.049256
std 0.399851
min -0.155273
25% -0.044578
50% 0.018517
75% 0.086474
max 4.718499
<AxesSubplot: ylabel='Density'>
将拟合的值逆向转化并可视化：
predictions_ARIMA_diff = pd.Series(
results_ARIMA.fittedvalues, copy=True
) # 复制Log和差分后的预测值
predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum() # 一阶差分的逆向转换为累计和
predictions_ARIMA_log = pd.Series(
ts_week_log.iloc[0], index=ts_week_log.index
) # 将一个新序列的每个元素都填充为ts_week_log的第一个元素
predictions_ARIMA_log = predictions_ARIMA_log.add(
predictions_ARIMA_diff_cumsum, fill_value=0
) # 将累计和追加到第一个元素上，还原log序列。其中的空值填充为0
predictions_ARIMA = np.exp(predictions_ARIMA_log) # 进行指数运算，还原为原始值
plt.plot(ts_week)
plt.plot(predictions_ARIMA)
plt.title("RMSE:%.4f" % np.sqrt(sum((predictions_ARIMA - ts_week) ** 2) / len(ts_week)))
plt.show()
from sklearn.metrics import mean_squared_error
pdq = (4, 1, 4)
test_point = 15 # 测试点数
size = len(ts_week_log) - 15
train, test = ts_week_log[0:size], ts_week_log[size : len(ts_week_log)] # 切分数据集
model = sm.tsa.arima.ARIMA(train, order=pdq)
model_fit = model.fit()
predictions = model_fit.predict(
start=size, end=size + test_point - 1, dynamic=True
).tolist()
error = mean_squared_error(test, predictions)
print("Printing Mean Squared Error of Predictions...")
predictions_series = pd.Series(
predictions, index=test.index
) # predictions没有Index，无法在图上显示出来，必须先追加Index
fig, ax = plt.subplots()
ax.set(title="Spot Exchange Rate, Euro into USD", xlabel="Date", ylabel="Euro into USD")
ax.plot(ts_week[-50:], "o", label="observed")
legend = ax.legend(loc="upper left")
legend.get_frame().set_facecolor("w")
Printing Mean Squared Error of Predictions...
SARIMAX
SARIMAX模型简介
ARIMA是用于单变量时间序列数据预测的最广泛使用的预测方法之一。尽管该方法可以处理具有趋势的数据，但它不支持具有季节性组件的时间序列。支持对该系列季节性组成部分进行直接建模的ARIMA扩展称为(SARIMA(Seasonal Auto Regressive Integrated Moving Averages with eXogenous regressors.)。它增加了三个新的参数来指定系列季节性成分的自回归（AR），差分（I）和移动平均（MA），以及季节性周期的附加参数。
ARIMA -> SARIMA -> SARIMAX：
S是Seasonal，就是季节性、周期性的意思
X是eXogenous，外部信息的意思
所以SARIMAX和ARIMA比就是增加了周期性，并且使用了外部信息来增强模型的预测能力
SARIMAX模型的使用
SARIMAX模型与ARIMA最大的不同是，会多出季节性参数。需要提供的参数为：
趋势参数：（与ARIMA模型相同）
p：趋势自回归阶数
d：趋势差分阶数
q：趋势移动平均阶数
季节性参数：
P：季节性自回归阶数
D：季节性差分阶数
Q：季节性移动平均阶数
m：单个季节期间的时间步数
季节性ARIMA模型使用等于季节数的滞后差异来消除加性季节效应。与滞后1差分去除趋势一样，滞后差分引入移动平均项。季节性ARIMA模型包括滞后S处的自回归和移动平均项。SARIMA模型的符号指定为：SARIMA(p,d,q)(P,D,Q)m，其中 m m m参数会影响 P P P， D D D和 Q Q Q参数。例如，月度数据的 m m m为12表示每年的季节周期。
如何确定季节性部分：
假如序列有显著是季节性模式，则需要用季节性差分，但是不要使用两次季节性差分或总过超过两次差分（季节性和非季节性）
假如差分序列在滞后s处有正的ACF，s是季节的长度，则需要加入SAR项；假如差分序列在滞后s处有负的ACF，则需要加入SMA项，如果使用了季节性差异，则后一种情况很可能发生，如果数据具有稳定和合乎逻辑的季节性模式。如果没有使用季节性差异，前者很可能会发生，只有当季节性模式不稳定时才适用。应该尽量避免在同一模型中使用多于一个或两个季节性参数（SAR + SMA），因为这可能导致过度拟合数据和/或估算中的问题。
SARIMAX模型实战
实验数据：美国夏威夷莫纳罗亚火山天文台空气样本中大气二氧化碳数据（该数据收集了1958年3月至2001年12月的二氧化碳样本数据。重要的是statsmodels中自带此部分数据）
加载数据与序列可视化
import itertools
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
plt.style.use("fivethirtyeight")
data = sm.datasets.co2.load_pandas().data
data = data["co2"].resample("MS").mean()
data = data.fillna(data.bfill())
print(data.head())
data.plot(figsize=(12, 6))
plt.show()
1958-03-01 316.100000
1958-04-01 317.200000
1958-05-01 317.433333
1958-06-01 315.625000
1958-07-01 315.625000
Freq: MS, Name: co2, dtype: float64
从绘制的图形看，整体的数据呈上升趋势，并且存在一定的季节性波动。
为模型选择最优参数
在上面的文章中详细介绍了ARIMA模型及，对应的(p,d,q）参数，方案是采用ACF和PACF图来确定具体的参数值。这里采用了另外一种方法“网格搜索”，通过迭代不同参数的组合来获取最佳参数。
# Define the p, d and q parameters to take any value between 0 and 2
p = d = q = range(0, 2)
# Generate all different combinations of p, q and q triplets
pdq = list(itertools.product(p, d, q))
print(pdq)
# Generate all different combinations of seasonal p, q and q triplets
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]
print(seasonal_pdq)
[(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1), (1, 1, 0), (1, 1, 1)]
[(0, 0, 0, 12), (0, 0, 1, 12), (0, 1, 0, 12), (0, 1, 1, 12), (1, 0, 0, 12), (1, 0, 1, 12), (1, 1, 0, 12), (1, 1, 1, 12)]
接下来我们就可以通过不同的参数组合来对SARIMAX进行评估。在评估和比较不同参数的模型时，我们可以使用AIC值来评估。该值将在模型中直接返回。在这里AIC的值越小，则模型越优。
import warnings
warnings.filterwarnings("ignore") # specify to ignore warning messages
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
model = sm.tsa.statespace.SARIMAX(
data,
order=param,
seasonal_order=param_seasonal,
enforce_stationarity=False,
enforce_invertibility=False,
)
results = model.fit()
print("ARIMA{}x{} - AIC:{}".format(param, param_seasonal, results.aic))
except:
continue
ARIMA(0, 0, 0)x(0, 0, 0, 12) - AIC:7612.583429881011
ARIMA(0, 0, 0)x(0, 0, 1, 12) - AIC:6787.343624032837
ARIMA(0, 0, 0)x(0, 1, 0, 12) - AIC:1854.8282341411787
ARIMA(0, 0, 0)x(0, 1, 1, 12) - AIC:1596.7111727642252
ARIMA(0, 0, 0)x(1, 0, 0, 12) - AIC:1058.9388921320046
ARIMA(0, 0, 0)x(1, 0, 1, 12) - AIC:1056.2878392776897
ARIMA(0, 0, 0)x(1, 1, 0, 12) - AIC:1361.657897807711
ARIMA(0, 0, 0)x(1, 1, 1, 12) - AIC:1044.764791294482
ARIMA(0, 0, 1)x(0, 0, 0, 12) - AIC:6881.048755640422
ARIMA(0, 0, 1)x(0, 0, 1, 12) - AIC:6072.662327853603
ARIMA(0, 0, 1)x(0, 1, 0, 12) - AIC:1379.1941067380717
ARIMA(0, 0, 1)x(0, 1, 1, 12) - AIC:1241.4174716779276
ARIMA(0, 0, 1)x(1, 0, 0, 12) - AIC:1089.645943257476
ARIMA(0, 0, 1)x(1, 0, 1, 12) - AIC:780.4318116365525
ARIMA(0, 0, 1)x(1, 1, 0, 12) - AIC:1119.5957893616242
ARIMA(0, 0, 1)x(1, 1, 1, 12) - AIC:807.0912988010436
ARIMA(0, 1, 0)x(0, 0, 0, 12) - AIC:1675.8086923024289
ARIMA(0, 1, 0)x(0, 0, 1, 12) - AIC:1240.221119919406
ARIMA(0, 1, 0)x(0, 1, 0, 12) - AIC:633.4425587178921
ARIMA(0, 1, 0)x(0, 1, 1, 12) - AIC:337.79385466014884
ARIMA(0, 1, 0)x(1, 0, 0, 12) - AIC:619.9501759055873
ARIMA(0, 1, 0)x(1, 0, 1, 12) - AIC:376.92837604498965
ARIMA(0, 1, 0)x(1, 1, 0, 12) - AIC:478.32969079023997
ARIMA(0, 1, 0)x(1, 1, 1, 12) - AIC:323.3246657433406
ARIMA(0, 1, 1)x(0, 0, 0, 12) - AIC:1371.1872602335284
ARIMA(0, 1, 1)x(0, 0, 1, 12) - AIC:1101.8410734304073
ARIMA(0, 1, 1)x(0, 1, 0, 12) - AIC:587.9479709956827
ARIMA(0, 1, 1)x(0, 1, 1, 12) - AIC:302.49490000862465
ARIMA(0, 1, 1)x(1, 0, 0, 12) - AIC:584.4333533219732
ARIMA(0, 1, 1)x(1, 0, 1, 12) - AIC:337.19990508929817
ARIMA(0, 1, 1)x(1, 1, 0, 12) - AIC:433.0863608054848
ARIMA(0, 1, 1)x(1, 1, 1, 12) - AIC:281.51900551307494
ARIMA(1, 0, 0)x(0, 0, 0, 12) - AIC:1676.8881767362063
ARIMA(1, 0, 0)x(0, 0, 1, 12) - AIC:1241.9354423854566
ARIMA(1, 0, 0)x(0, 1, 0, 12) - AIC:624.2602351237664
ARIMA(1, 0, 0)x(0, 1, 1, 12) - AIC:341.2896611169605
ARIMA(1, 0, 0)x(1, 0, 0, 12) - AIC:579.3896670400418
ARIMA(1, 0, 0)x(1, 0, 1, 12) - AIC:370.5914826521988
ARIMA(1, 0, 0)x(1, 1, 0, 12) - AIC:476.0500427293109
ARIMA(1, 0, 0)x(1, 1, 1, 12) - AIC:327.58085031215364
ARIMA(1, 0, 1)x(0, 0, 0, 12) - AIC:1372.6085881647123
ARIMA(1, 0, 1)x(0, 0, 1, 12) - AIC:1199.4888194376067
ARIMA(1, 0, 1)x(0, 1, 0, 12) - AIC:586.4485732576351
ARIMA(1, 0, 1)x(0, 1, 1, 12) - AIC:305.6273829633814
ARIMA(1, 0, 1)x(1, 0, 0, 12) - AIC:587.5093074043984
ARIMA(1, 0, 1)x(1, 0, 1, 12) - AIC:398.973550780749
ARIMA(1, 0, 1)x(1, 1, 0, 12) - AIC:433.5469464656961
ARIMA(1, 0, 1)x(1, 1, 1, 12) - AIC:284.35969839629036
ARIMA(1, 1, 0)x(0, 0, 0, 12) - AIC:1324.3111127324544
ARIMA(1, 1, 0)x(0, 0, 1, 12) - AIC:1060.9351914437214
ARIMA(1, 1, 0)x(0, 1, 0, 12) - AIC:600.7412682594972
ARIMA(1, 1, 0)x(0, 1, 1, 12) - AIC:312.1329633875623
ARIMA(1, 1, 0)x(1, 0, 0, 12) - AIC:593.6637754760593
ARIMA(1, 1, 0)x(1, 0, 1, 12) - AIC:349.2091457479196
ARIMA(1, 1, 0)x(1, 1, 0, 12) - AIC:440.1375883655929
ARIMA(1, 1, 0)x(1, 1, 1, 12) - AIC:293.7426223158166
ARIMA(1, 1, 1)x(0, 0, 0, 12) - AIC:1262.6545542464662
ARIMA(1, 1, 1)x(0, 0, 1, 12) - AIC:1052.0636724057658
ARIMA(1, 1, 1)x(0, 1, 0, 12) - AIC:581.3099934822465
ARIMA(1, 1, 1)x(0, 1, 1, 12) - AIC:295.9374057357466
ARIMA(1, 1, 1)x(1, 0, 0, 12) - AIC:576.8647111794348
ARIMA(1, 1, 1)x(1, 0, 1, 12) - AIC:327.9049121396483
ARIMA(1, 1, 1)x(1, 1, 0, 12) - AIC:444.1243686500529
ARIMA(1, 1, 1)x(1, 1, 1, 12) - AIC:277.78022026039173
从上面返回的结果中，我们发现ARIMA(1, 1, 1)x(1, 1, 1, 12)返回的AIC为277.78，值最小。因此，我们应该认为这是我们所考虑的所有模型中的最佳选择。
拟合ARIMA模型
通过“网格搜索”我们找到了最佳拟合模型的参数，接下里我们将最佳参数值输入到一个新的 SARIMAX 模型：
model = sm.tsa.statespace.SARIMAX(
data,
order=(1, 1, 1),
seasonal_order=(1, 1, 1, 12),
enforce_stationarity=False,
enforce_invertibility=False,
)
results = model.fit()
print(results.summary().tables[1])
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.3182 0.092 3.442 0.001 0.137 0.499
ma.L1 -0.6254 0.077 -8.163 0.000 -0.776 -0.475
ar.S.L12 0.0010 0.001 1.732 0.083 -0.000 0.002
ma.S.L12 -0.8769 0.026 -33.811 0.000 -0.928 -0.826
sigma2 0.0972 0.004 22.632 0.000 0.089 0.106
==============================================================================
从 SARIMAX 输出结果得到的 summary 属性返回大量信息，但我们将把注意力集中在coef列上。在这里每列的P值都比0.05小或接近0.05，因此我们可以保留模型中的所有权重。在拟合季节 ARIMA 模型(以及其他相关的模型)时，重要的是要运行模型诊断，以确保模型所做的假设是正确的。 plot_diagnostics()允许我们快速生成模型诊断并调查任何异常行为。
results.plot_diagnostics(figsize=(12, 12))
plt.show()
这里我们主要要确保残差稳定，且平均分布为0。上面的模型诊断表明：
右上角的红色KDE线和黄色N(0,1)线非常的接近（N(0,1)表示，平均值为0，标准差为1），表明残差是正态分布。
左下角对的QQ分位图表明，从N(0,1)的标准正态分布抽取的样本，残差(蓝点)的有序分布服从线性趋势。 同样，这是一个强烈的迹象表明，残差是正态分布。
随着时间的推移（左上图）残差不会显示任何明显的季节性，似乎是白噪声。这通过右下角的自相关（即相关图）证实了这一点，它表明时间序列的残差与其自身的滞后版本有很低的相关性。
这些观察结果使我们得出结论，我们的模型产生了令人满意的结果，可以帮助我们预测未来的数据。
备注：虽然我们得到了令人满意的拟合参数，但是ARIMA模型的一些参数还是可以改变，以改进我们的模型拟合。 例如，我们的网格搜索只考虑受限制的参数组合集，所以如果我们扩大网格搜索，我们可能会找到更好的模型。
验证预测
我们已经获得了我们时间序列的模型，现在可以用来产生预测。 我们首先将预测值与时间序列的实际值进行比较，这将有助于我们了解我们的预测的准确性。 get_prediction()和conf_int()属性允许我们获得时间序列预测的值和相关的置信区间。
pred = results.get_prediction(start=pd.to_datetime("1998-01-01"), dynamic=False)
pred_ci = pred.conf_int()
代码中设置从1998年1月开始预测，通过dynamic=False设置后续的每一个都是使用1998年1月之前的数据预测生产的。
我们可以绘制二氧化碳时间序列的实际值和预测值，以评估我们做得如何。
ax = data["1990":].plot(label="Observed", figsize=(12, 6))
pred.predicted_mean.plot(ax=ax, label="One-step ahead Forecast", alpha=0.7)
ax.fill_between(
pred_ci.index, pred_ci.iloc[:, 0], pred_ci.iloc[:, 1], color="k", alpha=0.2
)
ax.set_xlabel("Date")
ax.set_ylabel("CO2 Levels")
plt.legend()
plt.show()
总体而言，我们的预测值与真实值保持一致，呈现总体增长趋势。量化我们的预测的准确性也是有用的。 我们将使用MSE（均方误差）来评估预测的准确性。
data_forecasted = pred.predicted_mean
data_truth = data["1998-01-01":]
# Compute the mean square error
mse = ((data_forecasted - data_truth) ** 2).mean()
print("The Mean Squared Error of our forecasts is {}".format(round(mse, 2)))
The Mean Squared Error of our forecasts is 0.07
这里预测的MSE值为0.07，是接近0的非常低的值。0的MSE是估计器将以完美的精度预测参数的观测值，这将是一个理想的场景但通常不可能。在这种情况下，我们只使用时间序列中的信息到某一点，之后，使用先前预测时间点的值生成预测。然而，使用动态预测可以获得更好地表达我们的真实预测能力。 在下面的代码块中，我们指定从1998年1月起开始计算动态预测和置信区间。
pred_dynamic = results.get_prediction(
start=pd.to_datetime("1998-01-01"), dynamic=True, full_results=True
)
pred_dynamic_ci = pred_dynamic.conf_int()
ax = data["1990":].plot(label="Observed", figsize=(12, 6))
pred_dynamic.predicted_mean.plot(label="Dynamic Forecast", ax=ax)
ax.fill_between(
pred_dynamic_ci.index,
pred_dynamic_ci.iloc[:, 0],
pred_dynamic_ci.iloc[:, 1],
color="k",
alpha=0.25,
)
ax.fill_betweenx(
ax.get_ylim(), pd.to_datetime("1998-01-01"), data.index[-1], alpha=0.1, zorder=-1
)
ax.set_xlabel("Date")
ax.set_ylabel("CO2 Levels")
plt.legend()
plt.show()
我们看到即使使用动态预测，总体预测也是准确的。 所有预测值（红线）与实际值（蓝线）相当吻合，并且在我们预测的置信区间内。我们再次通过计算MSE量化我们预测的准确性：
# Extract the predicted and true values of our time series
data_forecasted = pred_dynamic.predicted_mean
data_truth = data["1998-01-01":]
# Compute the mean square error
mse = ((data_forecasted - data_truth) ** 2).mean()
print("The Mean Squared Error of our forecasts is {}".format(round(mse, 2)))
The Mean Squared Error of our forecasts is 1.01
从动态预测获得的预测值产生的MSE 为1.01。要比非动态要高，这个也在预料之中，原因是可供参考的历史数据较少。
采用动态或非动态的预测都证实了时间序列模型的有效。然而，关于时间序列预测的大部分兴趣是能够及时预测未来价值观。
预测及可视化
在本教程的最后一步，我们将介绍如何利用季节性ARIMA时间序列模型来预测未来的价值。 我们的时间序列对象的get_forecast()属性可以计算预先指定数量的“步数”的预测值。
# Get forecast 500 steps ahead in future
pred_uc = results.get_forecast(steps=500)
# Get confidence intervals of forecasts
pred_ci = pred_uc.conf_int()
# 我们可以使用此代码的输出绘制其未来值的时间序列和预测。
ax = data.plot(label="Observed", figsize=(12, 6))
pred_uc.predicted_mean.plot(ax=ax, label="Forecast")
ax.fill_between(
pred_ci.index, pred_ci.iloc[:, 0], pred_ci.iloc[:, 1], color="k", alpha=0.25
)
ax.set_xlabel("Date")
ax.set_ylabel("CO2 Levels")
plt.legend()
plt.show()
我们所作的预测和相关的置信区间现在都可以用来进一步理解时间序列，并预测将会发生什么。 我们的预测显示，预计时间序列将以稳定的速度持续增长。
4036
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
