​
目录

​
引言

金融时间序列分析在最近几十年一直是热门研究课题之一。在本指南中，我们将通过真实的金融数据来说明金融时间序列的典型事实。为了刻画这些事实，需要不同于 Box-Jenkins 模型的新模型。因此，ARCH 模型由 R. F. Engle 于 1982 年首次提出，此后被大量学者所拓展。我们还将演示如何使用 Python 及其库来实现 ARCH。

众所周知，许多时间序列都具有 ARCH 效应，即虽然（建模残差）序列是白噪声，但其平方序列可能存在自相关。更重要的是，在实践中发现大量金融时间序列具有这一特性，因此 ARCH 效应已成为金融时间序列的典型事实之一。

​
金融时间序列的典型事实

现在我们简要列出并描述金融收益序列的几个重要典型事实（特征）

  • 厚尾（重尾）：收益的分布密度函数通常比相应的正态分布密度函数的尾部更厚（更重）。

  • ARCH 效应：虽然收益序列通常可以看作是白噪声，但其平方（和绝对值）序列通常可能存在自相关，并且这些自相关几乎不是负的。

  • 波动性聚集：大的收益变化倾向于在时间上聚集，小的变化倾向于跟随小的变化。

  • 不对称性：众所周知，资产收益的分布略微负偏。一种可能的解释是交易者对不利信息的反应比对有利信息的反应更强烈。

​
ARCH 模型的定义

具体而言，我们给出 ARCH 模型的定义如下。

定义 1. 阶数为 p≥1p≥1 的 ARCH(p)\text{ARCH(p)} 模型形式如下

{Xt=σtεtσt2=ω+α1Xt−12+α2Xt−22+⋯+αpXt−p2 \begin{equation} \left\{ \begin{array}{ll} X_t =\sigma_t \varepsilon_t \\ \sigma_{t}^2 =\omega+ \alpha_1 X_{t-1}^2 + \alpha_2 X_{t-2}^2 + \cdots+ \alpha_p X_{t-p}^2 \\ \end{array} \right. \end{equation}

其中 ω≥0,αi≥0\omega ≥ 0, \alpha_i ≥ 0 和 αp>0\alpha_p > 0 是常数，εt∼iid(0,1)\varepsilon_t \sim iid(0, 1)，且 εt\varepsilon_t 独立于 {Xk;k≤t−1}\{X_k;k ≤ t − 1 \}。如果随机过程 XtX_t 满足方程 (1)，则称其为 ARCH(p)ARCH(p) 过程。

根据定义 1，σt2\sigma_{t}^2（以及 σt\sigma_t）独立于 εt\varepsilon_t。此外，通常进一步假设 εt∼N(0,1)\varepsilon_t \sim N(0, 1)。然而，有时为了捕捉金融时间序列的更多特征，我们需要进一步假设 εt\varepsilon_t 服从标准化（偏）学生 T 分布或广义误差分布。

设 Fs\mathscr{F}_s 表示由 {Xk;k≤s}\{X_k;k ≤ s \} 生成的信息集，即 σ(Xk;k≤s)\sigma(X_k;k ≤ s) 域。容易看出，对于任意 s<ts <t，Fs\mathscr{F}_s 独立于 εt\varepsilon_t。根据定义 1 和条件数学期望的性质，我们有

E(Xt∣Ft−1)=E(σtεt∣Ft−1)=σtE(εt∣Ft−1)=σtE(εt)=0 \begin{equation} E(X_t|\mathscr{F}_{t−1}) = E(\sigma_t \varepsilon_t|\mathscr{F}_{t−1}) = \sigma_t E( \varepsilon_t|\mathscr{F}_{t−1}) = \sigma_t E(\varepsilon_t) = 0 \tag 2 \end{equation}

和

Var(Xt2∣Ft−1)=E(Xt2∣Ft−1)=E(σt2εt2∣Ft−1)=σt2E(εt2∣Ft−1)=σt2E(εt2)=σt2. \text{Var}(X_{t}^2| \mathscr{F}_{t−1}) = E(X_{t}^2|\mathscr{F}_{t−1}) = E(\sigma_{t}^2 \varepsilon_{t}^2|\mathscr{F}_{t−1}) = \sigma_{t}^2 E(\varepsilon_{t}^2|\mathscr{F}_{t−1}) = \sigma_{t}^2 E(\varepsilon_{t}^2) = \sigma_{t}^2.

这意味着 σt2\sigma_{t}^2 是 XtX_t 的条件方差，并且它根据 {Xk2;t−p≤k≤t−1}\{X_{k}^2; t −p ≤ k ≤ t −1\} 的先前值演变，类似于 AR(p)\text{AR}(p) 模型。因此，模型 (1) 被命名为 ARCH(p)\text{ARCH}(p) 模型。

作为 ARCH(p) 模型的一个示例，我们考虑 ARCH(1) 模型

{Xt=σtεtσt2=ω+α1Xt−12 \begin{equation} \left\{ \begin{array}{ll} \tag 3 X_t =\sigma_t \varepsilon_t \\ \sigma_{t}^2 =\omega+ \alpha_1 X_{t-1}^2 \\ \end{array} \right. \end{equation}

明确地说，无条件均值 E(Xt)=E(σtεt)=E(σt)E(εt)=0.E(X_t) = E(\sigma_t \varepsilon_t) = E(\sigma_t) E(\varepsilon_t) = 0.。

此外，ARCH(1) 模型可以表示为

Xt2=σt2+Xt2−σt2=ω+α1Xt−12+σt2εt2−σt2=ω+α1Xt2+ηtX_{t}^2 =\sigma_{t}^2 +X_{t}^2 − \sigma_{t}^2 =\omega +\alpha_1 X_{t-1}^2 +\sigma_{t}^2 \varepsilon_{t}^2 −\sigma_{t}^2 =\omega +\alpha_1 X_{t}^2 +\eta_t

即，

Xt2=ω+α1Xt2+ηt \begin{equation} X_{t}^2 =\omega +\alpha_1 X_{t}^2 +\eta_t \tag 4 \end{equation}

其中 ηt=σt2(εt2−1)\eta_t = \sigma_{t}^2(\varepsilon_{t}^2 − 1)。可以证明 $\eta_t$ 是一个新的白噪声，这留给读者作为练习。因此，如果 $0 < \alpha_1 < 1$，方程 (4) 是序列 Xt2 的平稳 AR(1) 模型。因此，无条件方差

Var(Xt)=E(Xt2)=E(ω+α1Xt−12+ηt)=ω+α1E(Xt2),Var ( X_t ) = E( X_{t}^2 ) = E(\omega+ \alpha_1 X_{t-1}^2 + \eta_t ) = \omega+ \alpha_1 E( X_{t}^2 ) ,

即 Var(Xt)=E(Xt2)=ω1−α1Var (X_t) = E (X_{t}^2 ) =\frac{\omega}{1-\alpha_1}

此外，对于 $h > 0$，根据条件数学期望的性质和 (2) 式，我们有

E(Xt+hXt)=E(E(Xt+hXt∣Ft+h−1))=E(XtE(Xt+h∣Ft+h−1))=0.E(X_{t+h} X_t) = E(E(X_{t+h} X_t|\mathscr{F}_{t+h-1})) = E(X_t E(X_{t+h}|\mathscr{F}_{t+h-1})) = 0.

综上所述，如果 $0 < \alpha_1 < 1$，我们有

  • 任何由公式 (3) 定义的 ARCH(1) 过程 $\{X_t \}$ 服从白噪声 $WN(0, \omega/(1 − \alpha_1))$。

  • 由于 $X_{t}^2$ 是由 (4) 定义的 AR(1) 过程，因此 $\operatorname{Corr}(X_{t}^2,X_{t+h}^2) = \alpha_{1}^{|h|} > 0$，这揭示了 ARCH 效应。

  • 显然，$E(\eta_t|\mathscr{F}_s)=0$ 对于任何 $t>s$ 成立，且结合方程 (4)，对于任何 $k>1$： $\operatorname{Var}(X_{t+k}|\mathscr{F}_t)=E(X_{t+K}^2|\mathscr{F}_t)$ $= E(\omega + \alpha_1 X_{t+k-1}+ \eta_{t+k}|\mathscr{F}_t )$ $= \omega + \alpha_1 \operatorname{Var}(X_{t+k−1}|\mathscr{F}_t),$

这反映了波动率聚类现象，即大（小）波动之后是大（小）波动。

此外，我们可以证明由方程 (3) 定义的 Xt 具有比相应正态分布更重的尾部。最后，请注意 ARCH(1) 模型的这些性质可以推广到 ARCH(p) 模型。

​
自回归条件异方差 (ARCH) 模型的优点和缺点

优点缺点
- ARCH 模型对于金融时间序列中的波动率建模很有用，这对于投资决策和风险管理非常重要。- ARCH 模型假设预测误差是独立同分布的，这在某些情况下可能不切实际。
- ARCH 模型考虑了异方差性，这意味着它可以对随时间变化的方差的时间序列进行建模。- ARCH 模型可能难以拟合具有许多参数的数据，这可能需要大量数据或高级估计技术。
- ARCH 模型相对易于使用，可以使用标准计量经济学软件实现。- ARCH 模型没有考虑时间序列的均值和方差之间可能存在的关系，这在某些情况下可能很重要。

注意

ARCH 模型是金融时间序列中波动率建模的有用工具，但与任何计量经济学模型一样，它也有局限性，应根据所建模数据的具体特征谨慎使用。

​
自回归条件异方差 (ARCH) 应用

  • 金融 - ARCH 模型广泛应用于金融领域，对股票价格、汇率、利率等金融时间序列的波动率进行建模。

  • 经济学 - ARCH 模型可用于对 GDP、通货膨胀、失业率等经济数据的波动率进行建模。

  • 工程 - ARCH 模型可用于工程领域，对能源、气候、污染、工业生产等相关数据的波动率进行建模。

  • 社会科学 - ARCH 模型可用于社会科学领域，对人口学、健康、教育等相关数据的波动率进行建模。

  • 生物学 - ARCH 模型可用于生物学领域，对进化、遗传学、流行病学等相关数据的波动率进行建模。

​
加载库和数据

提示

需要 Statsforecast。要安装，请参阅 说明。

接下来，我们导入绘图库并配置绘图样式。

import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
plt.style.use('fivethirtyeight')
plt.rcParams['lines.linewidth'] = 1.5
dark_style = {
    'figure.facecolor': '#212946',
    'axes.facecolor': '#212946',
    'savefig.facecolor':'#212946',
    'axes.grid': True,
    'axes.grid.which': 'both',
    'axes.spines.left': False,
    'axes.spines.right': False,
    'axes.spines.top': False,
    'axes.spines.bottom': False,
    'grid.color': '#2A3459',
    'grid.linewidth': '1',
    'text.color': '0.9',
    'axes.labelcolor': '0.9',
    'xtick.color': '0.9',
    'ytick.color': '0.9',
    'font.size': 12 }
plt.rcParams.update(dark_style)

from pylab import rcParams
rcParams['figure.figsize'] = (18,7)

​
读取数据

让我们从 Yahoo Finance 网站获取 S&P500 股票数据。

import datetime

import pandas as pd
import time
import yfinance as yf

ticker = '^GSPC'
period1 = datetime.datetime(2015, 1, 1)
period2 = datetime.datetime(2023, 9, 22)
interval = '1d' # 1d, 1m

SP_500 = yf.download(ticker, start=period1, end=period2, interval=interval, progress=False)
SP_500 = SP_500.reset_index()

SP_500.head()
价格日期调整后收盘价收盘价最高价最低价开盘价成交量
股票代码^GSPC^GSPC^GSPC^GSPC^GSPC^GSPC
02015-01-02 00:00:00+00:002058.1999512058.1999512072.3601072046.0400392058.8999022708700000
12015-01-05 00:00:00+00:002020.5799562020.5799562054.4399412017.3399662054.4399413799120000
22015-01-06 00:00:00+00:002002.6099852002.6099852030.2500001992.4399412022.1500244460110000
32015-01-07 00:00:00+00:002025.9000242025.9000242029.6099852005.5500492005.5500493805480000
42015-01-08 00:00:00+00:002062.1398932062.1398932064.0800782030.6099852030.6099853934010000
df=SP_500[["Date","Close"]].copy()

StatsForecast 的输入始终是一个具有三列的长格式数据框：unique_id、ds 和 y

  • unique_id（字符串、整数或类别）表示序列的标识符。

  • ds（日期戳）列应为 Pandas 预期的格式，日期最好是 YYYY-MM-DD，时间戳最好是 YYYY-MM-DD HH:MM:SS。

  • y（数值）表示我们希望预测的测量值。

df["unique_id"]="1"
df.columns=["ds", "y", "unique_id"]
df.head()
dsyunique_id
02015-01-02 00:00:00+00:002058.1999511
12015-01-05 00:00:00+00:002020.5799561
22015-01-06 00:00:00+00:002002.6099851
32015-01-07 00:00:00+00:002025.9000241
42015-01-08 00:00:00+00:002062.1398931
print(df.dtypes)

ds           datetime64[ns]
y                   float64
unique_id            object
dtype: object

​
使用 plot 方法探索数据

使用 StatsForecast 类的 plot 方法绘制序列。此方法会从数据集中随机打印一个序列，对于基本探索性数据分析 (EDA) 非常有用。

from statsforecast import StatsForecast

StatsForecast.plot(df)

​
增广迪基-福勒检验

增广迪基-福勒 (ADF) 检验是一种统计检验，用于确定时间序列数据中是否存在单位根。单位根可能导致时间序列分析中出现不可预测的结果。在单位根检验中，会建立一个零假设，以确定时间序列数据受趋势影响的强度。通过接受零假设，我们接受时间序列数据不平稳的证据。通过拒绝零假设或接受备择假设，我们接受时间序列数据是由平稳过程生成的证据。此过程也称为平稳趋势。ADF 检验统计量的值为负。ADF 值越低，表示对零假设的拒绝越强。

增广迪基-福勒检验是一种常用的统计检验，用于检验给定时间序列是否平稳。我们可以通过定义零假设和备择假设来实现这一点。

零假设：时间序列不平稳。它给出了一个时间相关的趋势。备择假设：时间序列平稳。换句话说，该序列不依赖于时间。

ADF 或 t 统计量 < 临界值：拒绝零假设，时间序列是平稳的。ADF 或 t 统计量 > 临界值：未能拒绝零假设，时间序列是不平稳的。

让我们检查一下我们正在分析的序列是否是平稳序列。让我们创建一个函数，使用 Dickey Fuller 检验进行检查。

from statsmodels.tsa.stattools import adfuller

def Augmented_Dickey_Fuller_Test_func(series , column_name):
    print (f'Dickey-Fuller test results for columns: {column_name}')
    dftest = adfuller(series, autolag='AIC')
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','No Lags Used','Number of observations used'])
    for key,value in dftest[4].items():
       dfoutput['Critical Value (%s)'%key] = value
    print (dfoutput)
    if dftest[1] <= 0.05:
        print("Conclusion:====>")
        print("Reject the null hypothesis")
        print("The data is stationary")
    else:
        print("Conclusion:====>")
        print("The null hypothesis cannot be rejected")
        print("The data is not stationary")

Augmented_Dickey_Fuller_Test_func(df["y"],'S&P500')

Dickey-Fuller test results for columns: S&P500
Test Statistic          -0.814971
p-value                  0.814685
No Lags Used            10.000000
                          ...    
Critical Value (1%)     -3.433341
Critical Value (5%)     -2.862861
Critical Value (10%)    -2.567473
Length: 7, dtype: float64
Conclusion:====>
The null hypothesis cannot be rejected
The data is not stationary

在前面的结果中，我们可以看到 Augmented_Dickey_Fuller 检验给出了 p-value 0.864700，这告诉我们不能拒绝零假设，另一方面，我们序列的数据是不平稳的。

我们需要对时间序列进行差分，以便将数据转换为平稳数据。

​
收益率序列

自 20 世纪 70 年代以来，随着计算机和互联网技术的发展，金融业非常繁荣。金融产品（包括各种衍生品）的交易产生了大量数据，形成了金融时间序列。对于金融而言，金融产品的收益率是最令人感兴趣的，因此我们关注的重点是收益率序列。如果 $P_t$ 是某个金融产品在时间 t 的收盘价，则该产品的收益率是

Xt=(Pt−Pt−1)Pt−1≈log(Pt)−log(Pt−1).X_t = \frac{(P_t − P_{t−1})}{P_{t−1}} ≈ log(P_t ) − log(P_{t−1}).log(Pt​)−log(Pt−1​).

对收益率序列 $\{X_t \}$ 进行了大量独立研究。并总结了许多金融工具、市场和时间段共同的重要特征。请注意，如果您购买金融产品，则该产品成为您的资产，其收益率成为您的资产收益率。现在让我们看以下示例。

我们可以使用 pandas 的 DataFrame.pct_change() 函数来估算收益率序列。pct_change() 函数有一个 periods 参数，其默认值为 1。如果要计算 30 天收益率，则必须将该值更改为 30。

df['return'] = 100 * df["y"].pct_change()
df.dropna(inplace=True, how='any')
df.head()
dsyunique_id收益率
12015-01-05 00:00:00+00:002020.5799561-1.827811
22015-01-06 00:00:00+00:002002.6099851-0.889347
32015-01-07 00:00:00+00:002025.90002411.162984
42015-01-08 00:00:00+00:002062.13989311.788828
52015-01-09 00:00:00+00:002044.8100591-0.840381
import plotly.express as px
fig = px.line(df, x=df["ds"], y="return",title="SP500 Return Chart",template = "plotly_dark")
fig.show()

​
创建平方收益率

df['sq_return'] = df["return"].mul(df["return"])
df.head()
dsyunique_id收益率平方收益率
12015-01-05 00:00:00+00:002020.5799561-1.8278113.340891
22015-01-06 00:00:00+00:002002.6099851-0.8893470.790938
32015-01-07 00:00:00+00:002025.90002411.1629841.352532
42015-01-08 00:00:00+00:002062.13989311.7888283.199906
52015-01-09 00:00:00+00:002044.8100591-0.8403810.706240

​
收益率 vs 平方收益率

from plotly.subplots import make_subplots
import plotly.graph_objects as go

fig = make_subplots(rows=1, cols=2)

fig.add_trace(go.Scatter(x=df["ds"], y=df["return"],
                         mode='lines',
                         name='return'),
row=1, col=1
)


fig.add_trace(go.Scatter(x=df["ds"], y=df["sq_return"],
                         mode='lines',
                         name='sq_return'), 
    row=1, col=2
)

fig.update_layout(height=600, width=800, title_text="Returns vs Squared Returns", template = "plotly_dark")
fig.show()


from scipy.stats import probplot, moment
from statsmodels.tsa.stattools import adfuller, q_stat, acf
import numpy as np
import seaborn as sns

def plot_correlogram(x, lags=None, title=None):    
    lags = min(10, int(len(x)/5)) if lags is None else lags
    fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(14, 8))
    x.plot(ax=axes[0][0], title='Return')
    x.rolling(21).mean().plot(ax=axes[0][0], c='k', lw=1)
    q_p = np.max(q_stat(acf(x, nlags=lags), len(x))[1])
    stats = f'Q-Stat: {np.max(q_p):>8.2f}\nADF: {adfuller(x)[1]:>11.2f}'
    axes[0][0].text(x=.02, y=.85, s=stats, transform=axes[0][0].transAxes)
    probplot(x, plot=axes[0][1])
    mean, var, skew, kurtosis = moment(x, moment=[1, 2, 3, 4])
    s = f'Mean: {mean:>12.2f}\nSD: {np.sqrt(var):>16.2f}\nSkew: {skew:12.2f}\nKurtosis:{kurtosis:9.2f}'
    axes[0][1].text(x=.02, y=.75, s=s, transform=axes[0][1].transAxes)
    plot_acf(x=x, lags=lags, zero=False, ax=axes[1][0])
    plot_pacf(x, lags=lags, zero=False, ax=axes[1][1])
    axes[1][0].set_xlabel('Lag')
    axes[1][1].set_xlabel('Lag')
    fig.suptitle(title+ f'Dickey-Fuller: {adfuller(x)[1]:>11.2f}', fontsize=14)
    sns.despine()
    fig.tight_layout()
    fig.subplots_adjust(top=.9)

plot_correlogram(df["return"], lags=30, title="Time Series Analysis plot \n")

​
Ljung-Box 检验

Ljung-Box 是一种自相关性检验，我们可以结合 ACF 和 PACF 图使用它。Ljung-Box 检验接收我们的数据，可选地包括要检验的滞后值或要考虑的最大滞后值，以及是否计算 Box-Pierce 统计量。Ljung-Box 和 Box-Pierce 是两种相似的检验统计量 Q，它们与卡方分布进行比较，以确定序列是否为白噪声。我们可以在模型的残差上使用 Ljung-Box 检验来查找自相关性，理想情况下，我们的残差应该是白噪声。

  • 零假设：数据独立分布，无自相关。
  • 备择假设：数据不独立分布；它们表现出序列相关。

带有 Box-Pierce 选项的 Ljung-Box 检验将为每个滞后返回 Ljung-Box 检验统计量、Ljung-Box p 值、Box-Pierce 检验统计量和 Box-Pierce p 值。

如果 $p<\alpha (0.05)$，我们拒绝零假设。

from statsmodels.stats.diagnostic import acorr_ljungbox

ljung_res = acorr_ljungbox(df["return"], lags= 40, boxpierce=True)

ljung_res.head()
lb_statlb_pvaluebp_statbp_pvalue
149.2222732.285409e-1249.1551832.364927e-12
262.9913482.097020e-1462.8992342.195861e-14
363.9449448.433622e-1463.8506638.834380e-14
474.3436522.742989e-1574.2210242.911751e-15
580.2348627.494100e-1680.0934988.022242e-16

​
将数据分割成训练集和测试集

让我们将数据分成几个集合

  1. 用于训练 ARCH 模型的数据
  2. 用于测试模型的数据

对于测试数据，我们将使用最后 30 天的数据来测试和评估模型的性能。

df=df[["ds","unique_id","return"]]
df.columns=["ds", "unique_id", "y"]

train = df[df.ds<='2023-05-24'] # Let's forecast the last 30 days
test = df[df.ds>'2023-05-24']

train.shape, test.shape

((2112, 3), (82, 3))

现在让我们绘制训练数据和测试数据。

sns.lineplot(train,x="ds", y="y", label="Train")
sns.lineplot(test, x="ds", y="y", label="Test")
plt.show()

​
使用 StatsForecast 实现 ARCH

要了解更多关于 ARCH Model 函数的参数，请参见下文列表。有关更多信息，请访问文档。

p : int
    Number of lagged versions of the series.
alias : str
    Custom name of the model.
prediction_intervals : Optional[ConformalIntervals]
    Information to compute conformal prediction intervals.
    By default, the model will compute the native prediction
    intervals.

​
加载库

from statsforecast import StatsForecast 
from statsforecast.models import ARCH

​
构建模型

导入并实例化模型。设置参数有时比较棘手。大师 Rob Hyndmann 的这篇关于季节性周期的文章可能会有所帮助。season_length。

season_length = 7 # Daily data
horizon = len(test) # number of predictions biasadj=True, include_drift=True,

models = [ARCH(p=2)]

我们通过实例化一个带有以下参数的新 StatsForecast 对象来拟合模型

models：模型列表。从 models 中选择您需要的模型并导入。

  • freq: 一个字符串，表示数据的频率。（参见 pandas 的可用频率。）

  • n_jobs: n_jobs：int，用于并行处理的作业数，使用 -1 表示所有核心。

  • fallback_model: 如果模型失败时使用的模型。

任何设置都传递给构造函数。然后调用其 fit 方法并传入历史数据框。

sf = StatsForecast(models=models,
                   freq='C', # custom business day frequency
                   )

​
拟合模型

sf.fit(df=train)

StatsForecast(models=[ARCH(2)])

让我们看看 ARCH 模型的结果。我们可以使用以下指令观察它：

result=sf.fitted_[0,0].model_
result

{'p': 2,
 'q': 0,
 'coeff': array([0.44321058, 0.34706751, 0.35172097]),
 'message': 'Optimization terminated successfully',
 'y_vals': array([-1.12220267, -0.73186003]),
 'sigma2_vals': array([1.38768694,        nan, 1.89278112, ..., 0.76423271, 0.45064684,
        0.88037072]),
 'fitted': array([        nan,         nan,  2.23474807, ..., -1.48033228,
         1.10018999, -0.98050166]),
 'actual_residuals': array([        nan,         nan, -1.07176381, ...,  1.49583575,
        -2.22239266,  0.24864162])}

现在让我们可视化模型的残差。

如我们所见，上面获得的结果是字典形式的输出，要从字典中提取每个元素，我们将使用 .get() 函数提取元素，然后将其保存在 pd.DataFrame() 中。

residual=pd.DataFrame(result.get("actual_residuals"), columns=["residual Model"])
residual
残差模型
0NaN
1NaN
2-1.071764
……
21091.495836
2110-2.222393
21110.248642
import scipy.stats as stats

fig, axs = plt.subplots(nrows=2, ncols=2)

# plot[1,1]
residual.plot(ax=axs[0,0])
axs[0,0].set_title("Residuals");

# plot
sns.distplot(residual, ax=axs[0,1]);
axs[0,1].set_title("Density plot - Residual");

# plot
stats.probplot(residual["residual Model"], dist="norm", plot=axs[1,0])
axs[1,0].set_title('Plot Q-Q')

# plot
plot_acf(residual,  lags=35, ax=axs[1,1],color="fuchsia")
axs[1,1].set_title("Autocorrelation");

plt.show();

​
预测方法

如果您希望在具有多个序列或模型的产品设置中提高速度，我们建议使用 StatsForecast.forecast 方法，而不是 .fit 和 .predict。

主要区别在于 .forecast 不存储拟合值，并且在分布式环境中具有高度可伸缩性。

forecast 方法接受两个参数：预测未来 h（horizon）步和 level。

  • h (int): 表示未来 h 步的预测。在此示例中，为未来 12 个月。

  • level (list of floats): 此可选参数用于概率预测。设置预测区间的水平（或置信百分位数）。例如，level=[90] 表示模型期望真实值在 90% 的时间内位于该区间内。

此处的 forecast 对象是一个新的数据框，其中包含一列模型名称和 y hat 值，以及不确定性区间的列。根据您的计算机，此步骤大约需要 1 分钟。

Y_hat = sf.forecast(df=train, h=horizon, fitted=True)
Y_hat
unique_iddsARCH(2)
012023-05-25 00:00:00+00:001.681839
112023-05-26 00:00:00+00:00-0.777029
212023-05-29 00:00:00+00:00-0.677962
…………
7912023-09-13 00:00:00+00:000.695591
8012023-09-14 00:00:00+00:00-0.176075
8112023-09-15 00:00:00+00:00-0.158605
values=sf.forecast_fitted_values()
values.head()
unique_iddsyARCH(2)
012015-01-05 00:00:00+00:00-1.827811NaN
112015-01-06 00:00:00+00:00-0.889347NaN
212015-01-07 00:00:00+00:001.1629842.234748
312015-01-08 00:00:00+00:001.788828-0.667577
412015-01-09 00:00:00+00:00-0.840381-0.752438

使用 forecast 方法添加 95% 置信区间

sf.forecast(df=train, h=horizon, level=[95])
unique_iddsARCH(2)ARCH(2)-lo-95ARCH(2)-hi-95
012023-05-25 00:00:00+00:001.681839-0.4193263.783003
112023-05-26 00:00:00+00:00-0.777029-3.9390542.384996
212023-05-29 00:00:00+00:00-0.677962-3.9072622.551338
………………
7912023-09-13 00:00:00+00:000.695591-0.9375852.328766
8012023-09-14 00:00:00+00:00-0.176075-1.4053591.053210
8112023-09-15 00:00:00+00:00-0.158605-1.3819151.064705
# Merge the forecasts with the true values
Y_hat1 = test.merge(Y_hat, how='left', on=['unique_id', 'ds'])
Y_hat1
dsunique_idyARCH(2)
02023-05-25 00:00:00+00:0010.8757581.681839
12023-05-26 00:00:00+00:0011.304909-0.777029
22023-05-30 00:00:00+00:0010.001660-0.968703
……………
792023-09-19 00:00:00+00:001-0.215101NaN
802023-09-20 00:00:00+00:001-0.939479NaN
812023-09-21 00:00:00+00:001-1.640093NaN
# Merge the forecasts with the true values

fig, ax = plt.subplots(1, 1)
plot_df = pd.concat([train, Y_hat1]).set_index('ds')
plot_df[['y', "ARCH(2)"]].plot(ax=ax, linewidth=2)
ax.set_title(' Forecast', fontsize=22)
ax.set_ylabel('Year ', fontsize=20)
ax.set_xlabel('Timestamp [t]', fontsize=20)
ax.legend(prop={'size': 15})
ax.grid(True)
plt.show()

​
带有置信区间的预测方法

使用 predict 方法生成预测。

predict 方法接受两个参数：预测未来的 h（horizon）和 level。

  • h (int): 表示未来 h 步的预测。在此示例中，为未来 12 个月。

  • level (list of floats): 此可选参数用于概率预测。设置预测区间的水平（或置信百分位数）。例如，level=[95] 表示模型期望真实值在 95% 的时间内位于该区间内。

此处的 forecast 对象是一个新的数据框，其中包含一列模型名称和 y hat 值，以及不确定性区间的列。

此步骤应不到 1 秒。

sf.predict(h=horizon)
unique_iddsARCH(2)
012023-05-25 00:00:00+00:001.681839
112023-05-26 00:00:00+00:00-0.777029
212023-05-29 00:00:00+00:00-0.677962
…………
7912023-09-13 00:00:00+00:000.695591
8012023-09-14 00:00:00+00:00-0.176075
8112023-09-15 00:00:00+00:00-0.158605
forecast_df = sf.predict(h=horizon, level=[80,95]) 
forecast_df
unique_iddsARCH(2)ARCH(2)-lo-95ARCH(2)-lo-80ARCH(2)-hi-80ARCH(2)-hi-95
012023-05-25 00:00:00+00:001.681839-0.4193260.3079613.0557163.783003
112023-05-26 00:00:00+00:00-0.777029-3.939054-2.8445661.2905082.384996
212023-05-29 00:00:00+00:00-0.677962-3.907262-2.7894881.4335642.551338
……………………
7912023-09-13 00:00:00+00:000.695591-0.937585-0.3722851.7634672.328766
8012023-09-14 00:00:00+00:00-0.176075-1.405359-0.9798600.6277111.053210
8112023-09-15 00:00:00+00:00-0.158605-1.381915-0.9584850.6412741.064705

我们可以使用 pandas 函数 pd.concat() 将预测结果与历史数据连接起来，然后使用该结果进行绘图。

df_plot=pd.concat([df, forecast_df]).set_index('ds').tail(220)
df_plot
unique_idyARCH(2)ARCH(2)-lo-95ARCH(2)-lo-80ARCH(2)-hi-80ARCH(2)-hi-95
ds
2023-03-07 00:00:00+00:001-1.532692NaNNaNNaNNaNNaN
2023-03-08 00:00:00+00:0010.141479NaNNaNNaNNaNNaN
2023-03-09 00:00:00+00:001-1.845936NaNNaNNaNNaNNaN
……………………
2023-09-13 00:00:00+00:001NaN0.695591-0.937585-0.3722851.7634672.328766
2023-09-14 00:00:00+00:001NaN-0.176075-1.405359-0.9798600.6277111.053210
2023-09-15 00:00:00+00:001NaN-0.158605-1.381915-0.9584850.6412741.064705
sf.plot(train, test.merge(forecast_df), level=[80, 95], max_insample_length=120)

让我们使用 Statsforecast 中自带的 plot 函数绘制相同的图，如下所示。

​
交叉验证

在前面的步骤中，我们使用历史数据来预测未来。但是，为了评估其准确性，我们还想知道模型在过去会表现如何。要评估模型在数据上的准确性和稳健性，请执行交叉验证。

对于时间序列数据，交叉验证是通过在历史数据上定义一个滑动窗口并预测其后续周期来完成的。这种形式的交叉验证使我们能够在更广泛的时间实例上更好地估计模型的预测能力，同时保持训练集中的数据是连续的，这是模型所要求的。

下图描绘了这种交叉验证策略

​
执行时间序列交叉验证

时间序列模型的交叉验证被认为是最佳实践，但大多数实现速度非常慢。Statsforecast 库将交叉验证实现为分布式操作，从而减少了执行过程的时间。如果数据集很大，您还可以使用 Ray、Dask 或 Spark 在分布式集群中执行交叉验证。

在此示例中，我们希望评估模型在过去 5 个月（n_windows=5）的性能，每隔 12 个月（step_size=12）进行预测。根据您的计算机，此步骤大约需要 1 分钟。

StatsForecast 类的 cross_validation 方法接受以下参数。

  • df: 训练数据框

  • h (int): 表示预测未来 h 步。在此示例中，为未来 12 个月。

  • step_size (int): 每个窗口之间的步长。换句话说：您希望多久运行一次预测过程。

  • n_windows(int): 用于交叉验证的窗口数。换句话说：您希望评估过去多少次预测过程。

crossvalidation_df = sf.cross_validation(df=train,
                                         h=horizon,
                                         step_size=6,
                                         n_windows=5)

crossvaldation_df 对象是一个新的数据框，包含以下列：

  • unique_id: 序列标识符
  • ds: 日期戳或时间索引
  • cutoff: n_windows 的最后一个日期戳或时间索引。
  • y: 真实值
  • "model": 包含模型名称和拟合值的列。
crossvalidation_df
unique_iddscutoffyARCH(2)
012022-12-21 00:00:00+00:002022-12-20 00:00:00+00:001.4867991.382105
112022-12-22 00:00:00+00:002022-12-20 00:00:00+00:00-1.445170-0.651618
212022-12-23 00:00:00+00:002022-12-20 00:00:00+00:000.586810-0.595213
………………
40712023-05-22 00:00:00+00:002023-01-26 00:00:00+00:000.0155030.693070
40812023-05-23 00:00:00+00:002023-01-26 00:00:00+00:00-1.122203-0.176181
40912023-05-24 00:00:00+00:002023-01-26 00:00:00+00:00-0.731860-0.157522

​
模型评估

现在我们将根据预测结果评估模型，我们将使用不同类型的指标 MAE、MAPE、MASE、RMSE、SMAPE 来评估准确性。

from functools import partial

from utilsforecast.evaluation import evaluate
from utilsforecast.losses import mae, mape, mase, rmse, smape

evaluate(
    test.merge(Y_hat),
    train_df=train,
    metrics=[mae, mape, partial(mase, seasonality=5), rmse, smape],
    agg_fn='mean',
)
metricARCH(2)
0mae0.949721
1mape11.789856
2mase0.875298
3rmse1.164914
4smape0.725702

​
参考资料

  1. Changquan Huang • Alla Petukhina. Springer series (2022). Applied Time Series Analysis and Forecasting with Python.
  2. Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica: Journal of the econometric society, 987-1007..
  3. James D. Hamilton. Time Series Analysis Princeton University Press, Princeton, New Jersey, 1st Edition, 1994.
  4. Nixtla Parameters.
  5. Pandas 可用频率.
  6. Rob J. Hyndman 和 George Athanasopoulos (2018). “Forecasting Principles and Practice (3rd ed)”.
  7. 季节性周期 - Rob J Hyndman.