ARIMA模型

简介

自回归积分滑动平均（ARIMA）模型由Box和Jenkins（1970年）推广普及。ARIMA模型将单变量时间序列描述为自回归（AR）和滑动平均（MA）滞后的组合，这些滞后捕捉了时间序列内的自相关性。积分阶数表示序列经过多少次差分以获得平稳序列。

我们为某个时间序列数据\(y_{t}\)编写一个\(ARIMA(p,d,q)\)模型，其中\(p\)是自回归滞后项的数量，\(d\)是差分阶数，\(q\)是移动平均滞后项的数量，表达式如下：

\[\Delta^{D}y_{t} = \sum^{p}_{i=1}\phi_{i}\Delta^{D}y_{t-i} + \sum^{q}_{j=1}\theta_{j}\epsilon_{t-j} + \epsilon_{t}\]

\[\epsilon_{t} \sim N\left(0,\sigma^{2}\right)\]

ARIMA模型与时间序列的Box-Jenkins方法相关联。根据这种方法，您应该对序列进行差分直到其平稳，然后使用信息准则和自相关图来选择\(ARIMA\)过程的适当滞后阶数。接着，您应用推断来获得潜在变量的估计值，并检查模型以查看模型是否捕捉到了时间序列中的自相关性。例如，您可以绘制模型残差的自相关图。一旦您满意，就可以使用该模型进行回顾和预测。

示例

我们将为年度太阳黑子数据运行一个ARIMA模型。首先加载数据：

import numpy as np
import pandas as pd
import pyflux as pf
from datetime import datetime
import matplotlib.pyplot as plt
%matplotlib inline

data = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/datasets/sunspot.year.csv')
data.index = data['time'].values

plt.figure(figsize=(15,5))
plt.plot(data.index,data['sunspot.year'])
plt.ylabel('Sunspots')
plt.title('Yearly Sunspot Data');

我们可以按照以下方式构建一个ARIMA模型，指定所需的模型阶数，以及包含数据的pandas DataFrame或numpy数组。这里我们指定了一个任意的\(ARIMA(4,0,4)\)模型：

model = pf.ARIMA(data=data, ar=4, ma=4, target='sunspot.year', family=pf.Normal())

接下来我们估计潜在变量。对于这个例子，我们将使用最大似然点质量估计 \(z^{MLE}\)：

x = model.fit("MLE")
x.summary()

ARIMA(4,0,4)
======================================== ============================================
Dependent Variable: sunspot.year         Method: MLE
Start Date: 1704                         Log Likelihood: -1189.488
End Date: 1988                           AIC: 2398.9759
Number of observations: 285              BIC: 2435.5008
=====================================================================================
Latent Variable      Estimate   Std Error  z        P>|z|    95% C.I.
==================== ========== ========== ======== ======== ========================
Constant             8.0092     3.2275     2.4816   0.0131   (1.6834 | 14.3351)
AR(1)                1.6255     0.0367     44.2529  0.0      (1.5535 | 1.6975)
AR(2)                -0.4345    0.2455     -1.7701  0.0767   (-0.9157 | 0.0466)
AR(3)                -0.8819    0.2295     -3.8432  0.0001   (-1.3317 | -0.4322)
AR(4)                0.5261     0.0429     12.2515  0.0      (0.4419 | 0.6103)
MA(1)                -0.5061    0.0383     -13.2153 0.0      (-0.5812 | -0.4311)
MA(2)                -0.481     0.1361     -3.533   0.0004   (-0.7478 | -0.2142)
MA(3)                0.2511     0.1093     2.2979   0.0216   (0.0369 | 0.4653)
MA(4)                0.2846     0.0602     4.7242   0.0      (0.1665 | 0.4027)
Sigma                15.7944
=====================================================================================

我们可以绘制潜在变量 \(z^{MLE}\)：使用 plot_z() 方法：

model.plot_z(figsize=(15,5))

我们可以使用plot_fit()来绘制样本内拟合图：

model.plot_fit(figsize=(15,10))

我们可以通过使用plot_predict_is()方法进行滚动样本内预测来了解模型的性能：

model.plot_predict_is(h=50, figsize=(15,5))

如果我们想要绘制预测结果，可以使用plot_predict()方法：

model.plot_predict(h=20,past_values=20,figsize=(15,5))

如果我们希望以DataFrame格式获取预测结果，可以直接使用predict()方法。

类描述

class ARIMA(data, ar, ma, integ, target, family)

自回归积分滑动平均模型（ARIMA）。

参数	类型	描述
data	pd.DataFrame or np.ndarray	包含单变量时间序列
ar	int	自回归滞后阶数
ma	int	移动平均滞后项的数量
integ	int	数据差分次数 (默认: 0)
target	string or int	指定使用DataFrame/array中的哪一列。
family	pf.Family instance	The distribution for the time series, e.g pf.Normal()

属性

latent_variables

一个包含模型潜在变量信息的pf.LatentVariables()对象，包括先验设置、任何拟合值、初始值和其他潜在变量信息。当模型被拟合时，这里就是潜在变量被更新/存储的地方。有关此对象内属性的信息以及访问潜在变量信息的方法，请参阅潜在变量文档。

方法

adjust_prior(index, prior)

调整模型潜在变量的先验分布。潜在变量及其索引可以通过打印附加到模型实例的latent_variables属性来查看。

参数	类型	描述
index	int	要更改的潜变量索引
prior	pf.Family instance	Prior distribution, e.g. pf.Normal()

返回: void - 修改模型的 latent_variables 属性

fit(method, **kwargs)

估计模型的潜在变量。用户选择一个推断选项，该方法会返回一个结果对象，同时更新模型的latent_variables属性。

参数	类型	描述
method	str	推断选项：例如 'M-H' 或 'MLE'

请参阅文档中的贝叶斯推断和经典推断部分，了解完整的推断选项列表。可以输入与所选特定推断模式相关的可选参数。

返回: 包含估计潜在变量信息的pf.Results实例

plot_fit(**kwargs)

绘制模型对数据的拟合情况。可选参数包括figsize，即绘图图形的尺寸。

返回 : void - 显示一个matplotlib绘图

plot_ppc(T, nsims)

绘制后验预测检查的直方图，使用用户选择的差异度量。此方法仅在使用贝叶斯推断进行拟合时有效。

参数	类型	描述
T	function	Discrepancy, e.g. np.mean or np.max
nsims	int	PPC需要进行多少次模拟

返回值: void - 显示一个matplotlib绘图

plot_predict(h, past_values, intervals, **kwargs)

绘制模型的预测结果，并附带置信区间。

参数	类型	描述
h	int	预测向前多少步
past_values	int	要绘制的历史数据点数量
intervals	boolean	是否绘制区间

可选参数包括figsize - 图表绘制的尺寸。请注意 如果您使用最大似然估计或变分推断，显示的区间将不会 反映潜在变量的不确定性。只有Metropolis-Hastings方法能提供完全贝叶斯 预测区间。由于平均场推断的局限性（无法考虑后验相关性）， 变分推断的贝叶斯区间不予显示。

返回 : void - 显示一个matplotlib绘图

plot_predict_is(h, fit_once, fit_method, **kwargs)

绘制模型在样本内的滚动预测。这意味着用户假装数据的最后一部分是样本外的，并在每个时间段后进行预测并评估其表现。用户可以选择是在开始时一次性拟合参数，还是在每个时间步都进行拟合。

参数	类型	描述
h	int	使用多少个先前的时间步
fit_once	boolean	是否只拟合一次，还是每个时间步都拟合
fit_method	str	选择哪种推断方法，例如'MLE'

可选参数包括figsize - 要绘制的图形尺寸。h是一个整数，表示要模拟性能的前几步。

返回 : void - 显示一个matplotlib绘图

plot_sample(nsims, plot_data=True)

绘制模型后验预测密度的样本。此方法仅适用于通过贝叶斯推断拟合的模型。

参数	类型	描述
nsims	int	需要抽取的样本数量
plot_data	boolean	是否同时绘制真实数据

返回 : void - 显示一个matplotlib绘图

plot_z(indices, figsize)

返回潜在变量及其相关不确定性的绘图。

参数	类型	描述
indices	int or list	要绘制的潜变量索引
figsize	tuple	matplotlib图形的大小

返回 : void - 显示一个matplotlib绘图

ppc(T, nsims)

返回后验预测检验的p值。此方法仅在您使用贝叶斯推断进行拟合时才有效。

参数	类型	描述
T	function	Discrepancy, e.g. np.mean or np.max
nsims	int	PPC需要进行多少次模拟

返回值: int - 差异检验的p值

predict(h, intervals=False)

返回模型预测结果的DataFrame。

参数	类型	描述
h	int	预测向前多少步
intervals	boolean	是否返回预测区间

请注意，如果您使用最大似然估计或变分推断，显示的区间将不会反映潜在变量的不确定性。只有Metropolis-Hastings方法能提供完全贝叶斯预测区间。由于平均场推断的局限性（未考虑后验相关性），变分推断的贝叶斯区间未予显示。

返回 : pd.DataFrame - 模型预测结果

predict_is(h, fit_once, fit_method)

返回模型样本内滚动预测的DataFrame。

参数	类型	描述
h	int	使用多少个先前的时间步
fit_once	boolean	是否只拟合一次，还是每个时间步都拟合
fit_method	str	选择哪种推断方法，例如'MLE'

返回 : pd.DataFrame - 模型预测结果

sample(nsims)

返回从后验预测密度中抽取的数据的np.ndarray数组。此方法仅在您使用贝叶斯推断拟合模型时才有效。

参数	类型	描述
nsims	int	进行多少次后验抽样

返回值 : np.ndarray - 来自后验预测密度的样本。

参考文献

Box, G; Jenkins, G. (1970). 时间序列分析：预测与控制. 旧金山: Holden-Day.

汉密尔顿, J.D. (1994). 时间序列分析. 泰勒 & 弗朗西斯美国.