GARCH模型


# 此单元格不会被渲染，但用于隐藏警告。

import warnings
warnings.filterwarnings("ignore")

import logging
logging.getLogger('statsforecast').setLevel(logging.ERROR)

import pandas as pd
pd.set_option('display.max_rows', 6)

使用Statsforecast的GARCH模型的逐步指南。

在本教程中，我们将熟悉主要的StatsForecast类及其一些相关方法，如StatsForecast.plot、StatsForecast.forecast和StatsForecast.cross_validation。

目录

介绍

广义自回归条件异方差（GARCH）模型是一种用于建模和预测金融与经济时间序列波动性的统计技术。该模型由罗伯特·恩格尔于1982年开发，是对1988年安德鲁·洛和克雷格·麦金莱提出的自回归条件异方差（ARCH）模型的扩展。

GARCH模型允许捕捉时间序列数据中的条件异方差性，即时间序列方差随时间变化而波动的特征。此特性在金融数据分析中尤为重要，因为波动性可以是衡量风险的一个重要指标。

GARCH模型已成为金融时间序列分析的基本工具，并广泛应用于从风险管理到预测股票及其他金融资产价格等各种应用中。

GARCH模型的定义

定义 1. 一个\(\text{GARCH}(p,q)\)模型，其阶数为\((p≥1,q≥0)\)，格式为

\[ \begin{equation} \begin{cases} X_t = \sigma_t \varepsilon_t\\ \sigma_{t}^2 = \omega + \sum_{i=1}^{p} \alpha_i X_{t-i}^2 + \sum_{j=1}^{q} \beta_j \sigma_{t-j}^2 \end{cases} \end{equation} \]

其中，\(\omega ≥0,\alpha_i ≥0,\beta_j ≥0,\alpha_p >0\)，且\(\beta_q >0\)是常数，\(\varepsilon_t \sim iid(0,1)\)，并且\(\varepsilon_t\)与\(\{X_k;k ≤ t − 1 \}\)独立。一个随机过程\(X_t\)如果满足方程(1)，则称其为\(\text{GARCH}(p, q )\)过程。

在实际应用中，发现对于某些时间序列，仅当阶数\(p\)较大时，(1)定义的\(\text{ARCH}(p)\)模型才能提供足够的拟合。通过允许过去的波动性影响当前的波动性，(1)可以得出一个更加简约的模型。这就是我们需要GARCH模型的原因。此外，请注意阶数\(p ≥ 1\)的条件。定义1中的GARCH模型具有如下性质。

命题 1. 如果 \(X_t\) 是在 (1) 中定义的 \(\text{GARCH}(p, q)\) 过程，并且 \(\sum_{i=1}^{p} \alpha_{i} + \sum_{j=1}^{q} \beta_j <1\)，那么以下命题成立。

  • \(X_{t}^2\) 服从 \(\text{ARMA}(m, q)\) 模型

\[X_{t}^2=\omega +\sum_{i=1}^{m} (\alpha_i + \beta_i )X_{t-i}^2 + \eta_t − \sum_{j=1}^q \beta_j \eta_{t-j}\]

其中，对于 \(i >p\)，\(\alpha_i =0\)；对于 \(j >q\)，\(\beta_j =0\)；\(m=\max(p,q)\)，且 \(\eta_t =\sigma_{t}^2 (\varepsilon_{t}^2 −1)\)。

  • \(X_t\) 是白噪声，其具有

\[E(X)=0, E(X_{t+h} X_t )=0 \ \ \text{对于} \ 任何 \ \ h \neq 0, Var(X_t)= \frac{\omega}{1-\sum_{i=1}^{m} (\alpha_i + \beta_i )}\]

  • \(\sigma_{t}^2\) 是 \(X_t\) 的条件方差，也就是说，我们有

\[E(X_t|\mathscr{F}_{t−1}) = 0, \sigma_{t}^2 = Var(X_{t}^2|\mathscr{F}_{t−1}).\]

  • 模型 (1) 反映了厚尾和波动聚集的特征。

虽然

广义自回归条件异方差（GARCH）模型的优缺点

优点 缺点
1. 灵活的模型：GARCH模型灵活，可以匹配不同类型的时间序列数据及其不同的波动模式。 1. 需要大量数据：GARCH模型需要大量数据才能准确估计模型参数。
2. 模拟波动的能力：GARCH模型能够模拟时间序列的波动性和异方差性，从而提高预测的准确性。 2. 对模型规范敏感：GARCH模型对模型规范非常敏感，如果规范不当，估计可能会很困难。
3. 包含历史信息：GARCH模型包含时间序列波动性的历史信息，这使得它在预测未来波动时非常有用。 3. 计算开销大：GARCH模型可能计算成本高，尤其是在使用更复杂模型时。
4. 允许包含外生变量：GARCH模型可以扩展以包含外生变量，这可以提高预测的准确性。 4. 不考虑极端事件：GARCH模型不考虑时间序列中的极端或意外事件，这可能会影响高波动情况下预测的准确性。
5. GARCH模型可以模拟条件异方差，即时间序列方差的变化是时间和时间序列自身以前值的函数。 5. GARCH模型假设时间序列误差服从正态分布，这在实践中可能并不成立。如果误差不服从正态分布，模型可能会产生不准确的波动性估计。
6. GARCH模型可用于估计投资组合的风险价值（VaR）和条件风险价值（CVaR）。

广义自回归条件异方差（GARCH）模型可以应用于多个领域

广义自回归条件异方差（GARCH）模型可以应用于多个需要对时间序列波动性进行建模和预测的领域。GARCH模型可以应用的一些领域包括：

  1. 金融市场： GARCH模型被广泛用于建模金融资产（如股票、债券、货币等）的收益波动（风险）。它允许捕捉波动性变化的特性。

  2. 商品价格： 原材料的价格，例如石油、黄金、谷物等，表现出可以用GARCH建模的条件波动性。

  3. 信用风险： 贷款和债券的违约风险在时间上也表现出与GARCH模型相适应的波动性。

  4. 经济时间序列： 宏观经济指标（如通货膨胀、GDP、失业率等）具有可用GARCH建模的条件波动性。

  5. 隐含波动性： GARCH模型可以估计金融期权中的隐含波动性。

  6. 预测： GARCH允许对任何时间序列进行条件波动性的预测。

  7. 风险分析： GARCH对于测量和管理投资组合和资产的风险非常有用。

  8. 金融： GARCH模型在金融领域被广泛用于建模金融资产（如股票、债券和货币）的价格波动性。

  9. 经济学： GARCH模型在经济学中用于建模商品和服务价格、通货膨胀和其他经济指标的波动性。

  10. 环境科学： GARCH模型在环境科学中用于建模温度、降水量和空气质量等变量的波动性。

  11. 社会科学： GARCH模型在社会科学中用于建模犯罪、迁移和就业等变量的波动性。

  12. 工程： GARCH模型在工程中用于建模电能需求、工业生产和车辆交通等变量的波动性。

  13. 健康科学： GARCH模型在健康科学中用于建模传染病病例数量和药品价格等变量的波动性。

GARCH模型适用于任何需要对时间序列中的异质条件波动性进行建模和预测的环境，尤其是在金融和经济学领域。

加载库和数据

Tip

将需要使用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)

读取数据

让我们从雅虎财经网站获取S&P500股票数据。

import pandas as pd
import time
from datetime import datetime
import yfinance as yf

ticker = '^GSPC'
period1 = int(time.mktime(datetime(2015, 1, 1, 23, 59).timetuple()))
period2 = int(time.mktime(datetime.now().timetuple()))
interval = '1d' # 1天，1个月

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

SP_500.head()
[*********************100%***********************]  1 of 1 completed
Date Open High Low Close Adj Close Volume
0 2015-01-02 2058.899902 2072.360107 2046.040039 2058.199951 2058.199951 2708700000
1 2015-01-05 2054.439941 2054.439941 2017.339966 2020.579956 2020.579956 3799120000
2 2015-01-06 2022.150024 2030.250000 1992.439941 2002.609985 2002.609985 4460110000
3 2015-01-07 2005.550049 2029.609985 2005.550049 2025.900024 2025.900024 3805480000
4 2015-01-08 2030.609985 2064.080078 2030.609985 2062.139893 2062.139893 3934010000 
df=SP_500[["Date","Close"]]

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()
ds y unique_id
0 2015-01-02 2058.199951 1
1 2015-01-05 2020.579956 1
2 2015-01-06 2002.609985 1
3 2015-01-07 2025.900024 1
4 2015-01-08 2062.139893 1 
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，这告诉我们原假设不能被拒绝；另一方面，我们的序列数据不是平稳的。

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

收益系列

自1970年代以来，金融行业伴随着计算机和互联网技术的进步而迅速繁荣。金融产品（包括各种衍生品）的交易产生了大量数据，这些数据形成了金融时间序列。在金融领域，金融产品的收益最为引人关注，因此我们将注意力集中在收益系列上。如果 \(P_t\) 是某一金融产品在时间 t 的收盘价，则该产品的收益为

\[X_t = \frac{(P_t − P_{t−1})}{P_{t−1}} ≈ \log(P_t) − \log(P_{t−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()
ds y unique_id return
1 2015-01-05 2020.579956 1 -1.827811
2 2015-01-06 2002.609985 1 -0.889347
3 2015-01-07 2025.900024 1 1.162984
4 2015-01-08 2062.139893 1 1.788828
5 2015-01-09 2044.810059 1 -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()
ds y unique_id return sq_return
1 2015-01-05 2020.579956 1 -1.827811 3.340891
2 2015-01-06 2002.609985 1 -0.889347 0.790938
3 2015-01-07 2025.900024 1 1.162984 1.352532
4 2015-01-08 2062.139893 1 1.788828 3.199906
5 2015-01-09 2044.810059 1 -0.840381 0.706240

收益与平方收益

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 检验，以寻找自相关，理想情况下我们的残差应该是白噪声。

  • Ho : 数据是独立分布的，没有自相关。
  • Ha : 数据不是独立分布的；它们表现出序列相关性。

带有 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_stat lb_pvalue bp_stat bp_pvalue
1 49.222273 2.285409e-12 49.155183 2.364927e-12
2 62.991348 2.097020e-14 62.899234 2.195861e-14
3 63.944944 8.433622e-14 63.850663 8.834380e-14
4 74.343652 2.742989e-15 74.221024 2.911751e-15
5 80.234862 7.494100e-16 80.093498 8.022242e-16

将数据分为训练集和测试集

让我们将数据分成以下两组： 1. 用于训练我们的 GARCH 模型的数据 2. 用于测试我们模型的数据

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

df=df[["ds","unique_id","return"]]
df.columns=["ds", "unique_id", "y"]
train = df[df.ds<='2023-05-31'] # Let's forecast the last 30 days
test = df[df.ds>'2023-05-31']
train.shape, test.shape
((2116, 3), (83, 3))

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

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

GARCH的实现与StatsForecast

要了解GARCH模型函数的参数，下面列出了相关内容。有关更多信息，请访问文档

p : int
    序列的滞后版本数量。
q: int
    波动性的滞后版本数量。
alias : str
    模型的自定义名称。
prediction_intervals : Optional[ConformalIntervals]
    计算符合预测区间的信息。
    默认情况下，模型将计算本地预测
    区间。

加载库

from statsforecast import StatsForecast 
from statsforecast.models import GARCH

实例化模型

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

season_length = 7 # 每日数据 
horizon = len(test) # 预测数量偏差调整=True，包含漂移=True，

models = [GARCH(1,1),
          GARCH(1,2),
          GARCH(2,2),
          GARCH(2,1),
          GARCH(3,1),
          GARCH(3,2),
          GARCH(3,3), 
          GARCH(1,3), 
          GARCH(2,3)]

我们通过实例化一个新的 StatsForecast 对象来拟合模型，使用以下参数：

models: 一个模型列表。选择您想要的模型并导入它们。

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

  • n_jobs: n_jobs: 整数，处理并行操作时使用的作业数量，使用 -1 表示所有核心。

  • fallback_model: 如果一个模型失败，则使用的备用模型。

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

sf = StatsForecast(df=train,
                   models=models,
                   freq='C', # 自定义营业日频率
                   n_jobs=-1)

交叉验证

我们构建了不同的GARCH模型，因此我们需要确定哪个是最佳模型，以便能够对其进行训练，从而进行预测。要知道哪个是最佳模型，我们需要进行交叉验证。

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

以下图表展示了这种交叉验证策略：

执行时间序列交叉验证

时间序列模型的交叉验证被认为是最佳实践，但大多数实现速度非常慢。statsforecast库将交叉验证实现为分布式操作，使得执行过程 less time-consuming。如果你有大型数据集，你还可以使用Ray、Dask或Spark在分布式集群中进行交叉验证。

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: 索引。如果您不喜欢使用索引，只需运行crossvalidation_df.resetindex()
  • ds: 日期戳或时间索引
  • cutoff: n_windows的最后日期戳或时间索引。
  • y: 真值
  • "model": 包含模型名称和拟合值的列。
crossvalidation_df
ds cutoff y GARCH(1,1) GARCH(1,2) GARCH(2,2) GARCH(2,1) GARCH(3,1) GARCH(3,2) GARCH(3,3) GARCH(1,3) GARCH(2,3)
unique_id
1 2023-01-03 2023-01-02 -0.404962 1.906913 1.886413 1.908357 1.930915 1.913208 1.904380 1.903494 1.885897 1.873233
1 2023-01-04 2023-01-02 -1.202064 -0.824837 -0.838750 -0.827394 -0.824686 -0.835927 -0.822887 -0.841358 -0.833009 -0.828155
1 2023-01-05 2023-01-02 1.746133 -0.665798 -0.653746 -0.666613 -0.680718 -0.663415 -0.677103 -0.649123 -0.661879 -0.665392
... ... ... ... ... ... ... ... ... ... ... ... ...
1 2023-05-29 2023-02-03 1.304909 -0.389434 -0.350367 -0.313831 -0.191451 -0.192649 -0.186990 -0.184408 -0.188771 -0.185092
1 2023-05-30 2023-02-03 0.001660 -0.340987 -0.296400 -0.278776 -0.158889 -0.157296 -0.160854 -0.154074 -0.156617 -0.158255
1 2023-05-31 2023-02-03 -0.610862 0.316639 0.274853 0.248995 0.135092 0.134911 0.132532 0.135533 0.134844 0.133779

415 rows × 12 columns

from datasetsforecast.losses import rmse

def compute_cv_rmse(crossvalidation_df):
    """计算所有生成模型的平均绝对误差（MAE）"""
    res = {}
    for mod in models: 
        res[mod] = rmse(crossvalidation_df['actual'], crossvalidation_df[str(mod)])
    return pd.Series(res)
crossvalidation_df.rename(columns = {'y' : 'actual'}, inplace = True) # 重命名实际值 
rmse_cv = crossvalidation_df.groupby(['unique_id', 'cutoff']).apply(compute_cv_rmse)

mae = rmse_cv.groupby('unique_id').mean()
mae.style.highlight_min(color = 'red', axis = 1)
  GARCH(1,1) GARCH(1,2) GARCH(2,2) GARCH(2,1) GARCH(3,1) GARCH(3,2) GARCH(3,3) GARCH(1,3) GARCH(2,3)
unique_id                  
1 1.598613 1.721620 1.547210 1.375646 1.442002 1.419195 1.373299 1.377253 1.374233

注意： 此结果可能会因您用于训练和测试模型的数据和时期以及您想要测试的模型而有所不同。这是一个示例，目的是能够教授使用 StatsForecast 的方法，特别是 GARCH 模型及在交叉验证中使用的参数，以确定此示例的最佳模型。

在之前的结果中，可以看到最佳模型是模型 \(\text{GARCH}(1,1)\)。

通过使用交叉验证找到的这个最佳模型结果，我们将继续训练我们的模型，然后进行预测。

拟合模型

season_length = 7 # 每日数据 
horizon = len(test) # 预测数量偏差调整=True，包含漂移=True，

models = [GARCH(1,1)]
sf = StatsForecast(df=train,
                   models=models,
                   freq='C', # 自定义工作日频率
                   n_jobs=-1)
sf.fit()
StatsForecast(models=[GARCH(1,1)])

让我们来看一下我们的Theta模型的结果。我们可以通过以下指令进行观察：

result=sf.fitted_[0,0].model_
result
{'p': 1,
 'q': 1,
 'coeff': array([0.11715665, 0.261296  , 0.64151786]),
 'message': 'Optimization terminated successfully',
 'y_vals': array([-0.61086244]),
 'sigma2_vals': array([0.81241364]),
 'fitted': array([        nan,  2.22688619, -0.75658963, ..., -0.20975516,
         0.83589047,  0.1360351 ]),
 'actual_residuals': array([        nan, -3.11623338,  1.91957389, ...,  1.51466384,
        -0.83423013, -0.74689754])}

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

正如我们所看到的，上面的结果以字典的形式输出，为了从字典中提取每个元素，我们将使用 .get() 函数提取元素，然后将其保存在 pd.DataFrame() 中。

residual=pd.DataFrame(result.get("actual_residuals"), columns=["residual Model"])
residual
residual Model
0 NaN
1 -3.116233
2 1.919574
... ...
2113 1.514664
2114 -0.834230
2115 -0.746898

2116 rows × 1 columns

from scipy import stats

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

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

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

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

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

plt.show();

预测方法

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

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

预测方法需要两个参数：预测接下来的 h（期望时间范围）和 level。

  • h (int): 表示预测未来 h 步。在这种情况下，表示提前 12 个月。

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

这里的预测对象是一个新的数据框，其中包括一个列，包含模型的名称和 y 预测值，以及不确定性区间的列。根据您的计算机，这一步骤大约需要 1 分钟。（如果您希望将时间缩短到几秒钟，请移除像 ARIMA 和 Theta 的自动模型。）

Y_hat = sf.forecast(horizon, fitted=True)

Y_hat.head()
ds GARCH(1,1)
unique_id
1 2023-06-01 1.393380
1 2023-06-02 -0.640599
1 2023-06-05 -0.508753
1 2023-06-06 -0.947621
1 2023-06-07 0.798889 
Y_hat = sf.forecast(horizon, fitted=True, level=[95])

Y_hat.head()
ds GARCH(1,1) GARCH(1,1)-lo-95 GARCH(1,1)-hi-95
unique_id
1 2023-06-01 1.393380 -0.048836 2.835595
1 2023-06-02 -0.640599 -2.789734 1.508536
1 2023-06-05 -0.508753 -2.327246 1.309740
1 2023-06-06 -0.947621 -2.476394 0.581153
1 2023-06-07 0.798889 -0.871354 2.469133 
values=sf.forecast_fitted_values()
values.head()
ds y GARCH(1,1) GARCH(1,1)-lo-95 GARCH(1,1)-hi-95
unique_id
1 2015-01-05 -1.827811 NaN NaN NaN
1 2015-01-06 -0.889347 2.226886 -0.858147 5.311919
1 2015-01-07 1.162984 -0.756590 -3.841623 2.328444
1 2015-01-08 1.788828 -0.636398 -3.721431 2.448635
1 2015-01-09 -0.840381 -1.472993 -4.558026 1.612040

添加95%的置信区间与预测方法

sf.forecast(h=horizon, level=[95])
ds GARCH(1,1) GARCH(1,1)-lo-95 GARCH(1,1)-hi-95
unique_id
1 2023-06-01 1.393380 -0.048836 2.835595
1 2023-06-02 -0.640599 -2.789734 1.508536
1 2023-06-05 -0.508753 -2.327246 1.309740
... ... ... ... ...
1 2023-09-21 -0.199650 -1.780158 1.380858
1 2023-09-22 -0.161219 -1.425179 1.102742
1 2023-09-25 0.136796 -0.916991 1.190583

83 rows × 4 columns

Y_hat=Y_hat.reset_index()
Y_hat.head()
unique_id ds GARCH(1,1) GARCH(1,1)-lo-95 GARCH(1,1)-hi-95
0 1 2023-06-01 1.393380 -0.048836 2.835595
1 1 2023-06-02 -0.640599 -2.789734 1.508536
2 1 2023-06-05 -0.508753 -2.327246 1.309740
3 1 2023-06-06 -0.947621 -2.476394 0.581153
4 1 2023-06-07 0.798889 -0.871354 2.469133 
# 将预测结果与真实值合并
test['unique_id'] = test['unique_id'].astype(int)
Y_hat1 = test.merge(Y_hat, how='left', on=['unique_id', 'ds'])
Y_hat1
ds unique_id y GARCH(1,1) GARCH(1,1)-lo-95 GARCH(1,1)-hi-95
0 2023-06-01 1 0.985445 1.393380 -0.048836 2.835595
1 2023-06-02 1 1.453442 -0.640599 -2.789734 1.508536
2 2023-06-05 1 -0.200358 -0.508753 -2.327246 1.309740
... ... ... ... ... ... ...
80 2023-09-26 1 -1.473453 NaN NaN NaN
81 2023-09-27 1 0.022931 NaN NaN NaN
82 2023-09-28 1 0.589317 NaN NaN NaN

83 rows × 6 columns

fig, ax = plt.subplots(1, 1)
plot_df = pd.concat([train, Y_hat1]).set_index('ds').tail(200)
plot_df[['y', "GARCH(1,1)"]].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()

带置信区间的预测方法

要生成预测，请使用预测方法。

预测方法接受两个参数：预测下一个 h（表示时间范围）和 level。

  • h (int)： 表示预测未来 h 步。在本例中，预测 30 天。

  • level (list of floats)： 该可选参数用于概率预测。设置预测区间的水平（或置信百分位）。例如，level=[95] 意味着模型预计真实值在该区间内的概率为 95%。

这里的预测对象是一个新的数据框，其中包含模型名称及 y hat 值的列，及不确定性区间的列。

此步骤应少于 1 秒。

sf.predict(h=horizon)
ds GARCH(1,1)
unique_id
1 2023-06-01 1.393380
1 2023-06-02 -0.640599
1 2023-06-05 -0.508753
... ... ...
1 2023-09-21 -0.199650
1 2023-09-22 -0.161219
1 2023-09-25 0.136796

83 rows × 2 columns

forecast_df = sf.predict(h=horizon, level=[80,95]) 

forecast_df.head(10)
ds GARCH(1,1) GARCH(1,1)-lo-95 GARCH(1,1)-lo-80 GARCH(1,1)-hi-80 GARCH(1,1)-hi-95
unique_id
1 2023-06-01 1.393380 -0.048836 0.450365 2.336394 2.835595
1 2023-06-02 -0.640599 -2.789734 -2.045843 0.764645 1.508536
1 2023-06-05 -0.508753 -2.327246 -1.697802 0.680296 1.309740
... ... ... ... ... ... ...
1 2023-06-12 -1.251548 -6.549860 -4.715928 2.212832 4.046764
1 2023-06-13 0.479689 -3.951083 -2.417437 3.376815 4.910460
1 2023-06-14 -0.318133 -3.508017 -2.403886 1.767620 2.871751

10 rows × 6 columns

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

df_plot=pd.concat([df, forecast_df]).set_index('ds').tail(220)
df_plot
unique_id y GARCH(1,1) GARCH(1,1)-lo-95 GARCH(1,1)-lo-80 GARCH(1,1)-hi-80 GARCH(1,1)-hi-95
ds
2023-03-15 1 -0.698088 NaN NaN NaN NaN NaN
2023-03-16 1 1.756201 NaN NaN NaN NaN NaN
2023-03-17 1 -1.101946 NaN NaN NaN NaN NaN
... ... ... ... ... ... ... ...
2023-09-21 NaN NaN -0.199650 -1.780158 -1.233088 0.833789 1.380858
2023-09-22 NaN NaN -0.161219 -1.425179 -0.987678 0.665241 1.102742
2023-09-25 NaN NaN 0.136796 -0.916991 -0.552238 0.825831 1.190583

220 rows × 7 columns

def plot_forecasts(y_hist, y_true, y_pred, models):
    _, ax = plt.subplots(1, 1, figsize = (20, 7))
    y_true = y_true.merge(y_pred, how='left', on=['unique_id', 'ds'])
    df_plot = pd.concat([y_hist, y_true]).set_index('ds').tail(12*10)
    df_plot[['y'] + models].plot(ax=ax, linewidth=2 )
    colors = ['green']
  # 指定图形特征：
    ax.fill_between(df_plot.index, 
                df_plot['GARCH(1,1)-lo-80'], 
                df_plot['GARCH(1,1)-hi-80'],
                alpha=.20,
                color='lime',
                label='GARCH(1,1)_level_80')
    ax.fill_between(df_plot.index, 
                df_plot['GARCH(1,1)-lo-95'], 
                df_plot['GARCH(1,1)-hi-95'],
                alpha=.2,
                color='white',
                label='GARCH(1,1)_level_95')
    ax.set_title('', fontsize=22)
    ax.set_ylabel("Return", fontsize=20)
    ax.set_xlabel('Month-Days', fontsize=20)
    ax.legend(prop={'size': 15})
    ax.grid(True)
    plt.show()
plot_forecasts(train, test, forecast_df, models=["GARCH(1,1)" ])

模型评估

现在我们可以使用适当的评估指标来计算预测的准确性。在这里，我们将使用均方根误差（RMSE）。为此，我们首先需要install datasetsforecast，这是一个由Nixtla开发的Python库，其中包含计算RMSE的函数。

%%capture
!pip install datasetsforecast

计算RMSE的函数接受两个参数：

  1. 实际值。
  2. 预测值，在本例中为 GARCH。

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

from datasetsforecast.losses import mae, mape, mase, rmse, smape
def evaluate_performace(y_hist, y_true, y_pred, model):
    y_true = y_true.merge(y_pred, how='left', on=['unique_id', 'ds'])
    evaluation = {}
    evaluation[model] = {}
    for metric in [mase, mae, mape, rmse, smape]:
        metric_name = metric.__name__
        if metric_name =='mase':
            evaluation[model][metric_name] = metric(y_true['y'].values, 
                                                y_true[model].values, 
                                                y_hist['y'].values, seasonality=7)
        else:
            evaluation[model][metric_name] = metric(y_true['y'].values, y_true[model].values)
    return pd.DataFrame(evaluation).T
evaluate_performace(train, test, Y_hat, model="GARCH(1,1)")
mae mape mase rmse smape
GARCH(1,1) 0.886834 372.562147 NaN 1.108626 136.528289

致谢

我们要感谢Naren Castellon撰写本教程。

参考文献

  1. Changquan Huang • Alla Petukhina. 施普林格系列（2022）。使用Python进行应用时间序列分析和预测。
  2. Bollerslev, T. (1986). 广义自回归条件异方差。计量经济学杂志，31(3)，307-327。
  3. Engle, R. F. (1982). 自回归条件异方差及对英国通货膨胀方差的估计。计量经济学：计量经济学会杂志，987-1007。
  4. James D. Hamilton. 时间序列分析 普林斯顿大学出版社，普林斯顿，新泽西，第一版，1994年。
  5. Nixtla 参数
  6. Pandas 可用频率
  7. Rob J. Hyndman 和 George Athanasopoulos (2018). “预测原则与实践，时间序列交叉验证”。
  8. 季节性周期 - Rob J Hyndman

Give us a ⭐ on Github