GARCH 模型

引言

GARCH（广义自回归条件异方差）模型是一种统计技术，用于模拟和预测金融和经济时间序列中的波动性。它由 Robert Engle 于 1982 年开发，是 Andrew Lo 和 Craig MacKinlay 于 1988 年提出的 ARCH（自回归条件异方差）模型的扩展。

GARCH 模型可以捕捉时间序列数据中存在的条件异方差性，即时间序列方差随时间波动的现象。这在金融数据分析中特别有用，因为波动性是衡量风险的重要指标。

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

GARCH 模型定义

定义 1. 阶为 $(p≥1,q≥0)$ 的 $\text{GARCH}(p,q)$ 模型形式如下

\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 )$ 过程。

实际上，对于某些时间序列，由 (1) 定义的 $\text{ARCH}(p)$ 模型只有在阶 $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}$

其中 $\alpha_i =0$ 对于 $i >p,βj =0$ 对于 $j >q,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{for} \ any \ \ 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) 反映了厚尾和波动性聚类。

尽管资产收益序列通常可以视为白噪声，但存在这样的收益序列可能具有自相关性。更重要的是，给定的原始时间序列不一定是收益序列，同时其值可能为负。如果时间序列具有自相关性，我们必须首先为该序列建立一个合适的模型（例如，一个 ARMA 模型），以便消除其中的任何自相关性。然后检查残差序列是否具有 ARCH 效应，如果具有，则进一步对残差建模。换句话说，如果时间序列 $Y_t$ 具有自相关性并具有 ARCH 效应，那么可以捕捉 $Y_t$ t 特征的 GARCH 模型应具有以下形式：

其中等式 (2) 被称为均值方程（模型），等式 (3) 被称为波动性（方差）方程（模型）， $Z_t$ 代表外生回归变量。如果 $Y_t$ 是收益序列，那么通常 $Y_t = r + X_t$ ，其中 $r$ 是一个常数，表示预期收益是固定的。

广义自回归条件异方差 (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 模型可以应用的一些领域包括：

金融市场：GARCH 模型广泛用于模拟股票、债券、货币等金融资产收益的波动性（风险）。它允许捕捉波动性的变化性质。
大宗商品价格：石油、黄金、谷物等原材料的价格表现出可以用 GARCH 建模的条件波动性。
信用风险：贷款和债券的违约风险随时间推移也呈现出波动性，GARCH 模型很适合对此进行建模。
经济时间序列：通货膨胀、GDP、失业率等宏观经济指标具有可以用 GARCH 建模的条件波动性。
隐含波动率：GARCH 模型可以估计金融期权中的隐含波动率。
预测：GARCH 允许对任何时间序列进行条件波动性预测。
风险分析：GARCH 有助于衡量和管理投资组合和资产的风险。
金融：GARCH 模型在金融领域广泛用于模拟股票、债券和货币等金融资产的价格波动性。
经济学：GARCH 模型在经济学中用于模拟商品和服务价格、通货膨胀和其他经济指标的波动性。
环境科学：GARCH 模型应用于环境科学，模拟温度、降水量和空气质量等变量的波动性。
社会科学：GARCH 模型在社会科学中用于模拟犯罪、移民和就业等变量的波动性。
工程：GARCH 模型应用于工程领域，模拟电力需求、工业生产和车辆交通等变量的波动性。
健康科学：GARCH 模型在健康科学中用于模拟传染病病例数和药品价格等变量的波动性。

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

加载库和数据

提示

需要用到 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
0	2015-01-02 00:00:00+00:00	2058.199951	2058.199951	2072.360107	2046.040039	2058.899902	2708700000
1	2015-01-05 00:00:00+00:00	2020.579956	2020.579956	2054.439941	2017.339966	2054.439941	3799120000
2	2015-01-06 00:00:00+00:00	2002.609985	2002.609985	2030.250000	1992.439941	2022.150024	4460110000
3	2015-01-07 00:00:00+00:00	2025.900024	2025.900024	2029.609985	2005.550049	2005.550049	3805480000
4	2015-01-08 00:00:00+00:00	2062.139893	2062.139893	2064.080078	2030.609985	2030.609985	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 00:00:00+00:00	2058.199951	1
1	2015-01-05 00:00:00+00:00	2020.579956	1
2	2015-01-06 00:00:00+00:00	2002.609985	1
3	2015-01-07 00:00:00+00:00	2025.900024	1
4	2015-01-08 00:00:00+00:00	2062.139893	1

print(df.dtypes)

ds           datetime64[ns, UTC]
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 的收盘价，那么该产品的收益为：

$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	收益率
1	2015-01-05 00:00:00+00:00	2020.579956	1	-1.827811
2	2015-01-06 00:00:00+00:00	2002.609985	1	-0.889347
3	2015-01-07 00:00:00+00:00	2025.900024	1	1.162984
4	2015-01-08 00:00:00+00:00	2062.139893	1	1.788828
5	2015-01-09 00:00:00+00:00	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	收益率	平方收益率
1	2015-01-05 00:00:00+00:00	2020.579956	1	-1.827811	3.340891
2	2015-01-06 00:00:00+00:00	2002.609985	1	-0.889347	0.790938
3	2015-01-07 00:00:00+00:00	2025.900024	1	1.162984	1.352532
4	2015-01-08 00:00:00+00:00	2062.139893	1	1.788828	3.199906
5	2015-01-09 00:00:00+00:00	2044.810059	1	-0.840381	0.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 检验来寻找自相关性，理想情况下，我们的残差应该是白噪声。

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), (78, 3))

使用 StatsForecast 实现 GARCH

加载库

from statsforecast import StatsForecast 
from statsforecast.models import GARCH

实例化模型

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

season_length = 7 # Dayly data 
horizon = len(test) # number of predictions biasadj=True, include_drift=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: int，并行处理中使用的作业数，使用 -1 表示所有核心。
fallback_model: 如果模型失败时使用的备用模型。

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

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

交叉验证

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

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

下图描述了这样的交叉验证策略

执行时间序列交叉验证

时间序列模型的交叉验证被认为是一种最佳实践，但大多数实现非常慢。statsforecast 库将交叉验证实现为分布式操作，从而使该过程花费的时间更少。如果您有大型数据集，您还可以使用 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: 序列标识符
ds: 日期戳或时间索引
cutoff: n_windows 的最后一个日期戳或时间索引。
y: 真实值
"model": 包含模型名称和拟合值的列。

crossvalidation_df

	unique_id	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)
0	1	2023-01-04 00:00:00+00:00	2023-01-03 00:00:00+00:00	0.753897	1.678755	1.678412	1.680475	1.686649	1.719494	2.210902	1.702743	1.647114	1.637795
1	1	2023-01-05 00:00:00+00:00	2023-01-03 00:00:00+00:00	-1.164553	-0.728069	-0.745487	-0.730648	-0.722156	-0.738119	-0.824748	-0.755277	-0.740976	-0.744150
2	1	2023-01-06 00:00:00+00:00	2023-01-03 00:00:00+00:00	2.284078	-0.589733	-0.582982	-0.590078	-0.598076	-0.587109	-0.866347	-0.571160	-0.587807	-0.584692
…	…	…	…	…	…	…	…	…	…	…	…	…	…
387	1	2023-05-26 00:00:00+00:00	2023-02-07 00:00:00+00:00	1.304909	-1.697814	-1.694747	-1.702537	-1.735631	-1.729903	-1.712997	-1.663399	-1.702160	-1.687723
388	1	2023-05-30 00:00:00+00:00	2023-02-07 00:00:00+00:00	0.001660	-0.326945	-0.337504	-0.329686	-0.330120	-0.334717	-0.327583	-0.330260	-0.338245	-0.332412
389	1	2023-05-31 00:00:00+00:00	2023-02-07 00:00:00+00:00	-0.610862	0.807625	0.787054	0.807819	0.841536	0.811702	0.836159	0.772193	0.801933	0.804526

from utilsforecast.evaluation import evaluate
from utilsforecast.losses import rmse

evals = evaluate(crossvalidation_df.drop(columns='cutoff'), metrics=[rmse], agg_fn='mean')
evals

	metric	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)
0	rmse	1.383143	1.526258	1.481056	1.389969	1.453538	1.539906	1.392352	1.515796	1.389061

evals.drop(columns='metric').loc[0].idxmin()

'GARCH(1,1)'

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

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

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

拟合模型

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

models = [GARCH(1,1)]

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

sf.fit(df=train)

StatsForecast(models=[GARCH(1,1)])

让我们看看 Theta 模型的结果。我们可以使用以下指令查看它

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

{'p': 1,
 'q': 1,
 'coeff': array([0.03745049, 0.18399111, 0.7890637 ]),
 'message': 'Optimization terminated successfully',
 'y_vals': array([-0.61086242]),
 'sigma2_vals': array([0.76298402]),
 'fitted': array([        nan,  2.14638896, -0.76426268, ..., -0.19747638,
         0.76993462,  0.13183178]),
 'actual_residuals': array([        nan, -3.03573613,  1.92724695, ...,  1.50238505,
        -0.7682743 , -0.7426942 ])}

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

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

residual=pd.DataFrame(result.get("actual_residuals"), columns=["residual Model"])
residual

	残差模型
0	NaN
1	-3.035736
2	1.927247
…	…
2113	1.502385
2114	-0.768274
2115	-0.742694

from scipy import 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（范围）和 level。

h (int): 表示未来 h 步的预测。在此示例中，提前 12 个月。
level (浮点数列表): 此可选参数用于概率预测。设置预测区间的水平（或置信百分位数）。例如，level=[90] 表示模型期望真实值有 90% 的时间落在此区间内。

此处的 forecast 对象是一个新的数据框，其中包含一列模型名称和 y 的预测值，以及不确定性区间的列。根据您的计算机性能，此步骤大约需要 1 分钟。（如果您想将速度提高到几秒钟，请移除 ARIMA 和 Theta 等 AutoModels）

Y_hat = sf.forecast(df=train, h=horizon, fitted=True)
Y_hat.head()

	unique_id	ds	GARCH(1,1)
0	1	2023-06-01 00:00:00+00:00	1.366914
1	1	2023-06-02 00:00:00+00:00	-0.593121
2	1	2023-06-05 00:00:00+00:00	-0.485200
3	1	2023-06-06 00:00:00+00:00	-0.927145
4	1	2023-06-07 00:00:00+00:00	0.766640

Y_hat = sf.forecast(df=train, h=horizon, fitted=True, level=[95])
Y_hat.head()

	unique_id	ds	GARCH(1,1)	GARCH(1,1)-lo-95	GARCH(1,1)-hi-95
0	1	2023-06-01 00:00:00+00:00	1.366914	-0.021035	2.754863
1	1	2023-06-02 00:00:00+00:00	-0.593121	-2.435497	1.249254
2	1	2023-06-05 00:00:00+00:00	-0.485200	-2.139216	1.168815
3	1	2023-06-06 00:00:00+00:00	-0.927145	-2.390566	0.536276
4	1	2023-06-07 00:00:00+00:00	0.766640	-0.771479	2.304759

values=sf.forecast_fitted_values()
values.head()

	unique_id	ds	y	GARCH(1,1)	GARCH(1,1)-lo-95	GARCH(1,1)-hi-95
0	1	2015-01-05 00:00:00+00:00	-1.827811	NaN	NaN	NaN
1	1	2015-01-06 00:00:00+00:00	-0.889347	2.146389	-0.972874	5.265652
2	1	2015-01-07 00:00:00+00:00	1.162984	-0.764263	-3.883526	2.355000
3	1	2015-01-08 00:00:00+00:00	1.788828	-0.650707	-3.769970	2.468556
4	1	2015-01-09 00:00:00+00:00	-0.840381	-1.449049	-4.568312	1.670214

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

sf.forecast(df=train, h=horizon, level=[95])

	unique_id	ds	GARCH(1,1)	GARCH(1,1)-lo-95	GARCH(1,1)-hi-95
0	1	2023-06-01 00:00:00+00:00	1.366914	-0.021035	2.754863
1	1	2023-06-02 00:00:00+00:00	-0.593121	-2.435497	1.249254
2	1	2023-06-05 00:00:00+00:00	-0.485200	-2.139216	1.168815
…	…	…	…	…	…
75	1	2023-09-14 00:00:00+00:00	-1.686546	-3.049859	-0.323233
76	1	2023-09-15 00:00:00+00:00	-0.322556	-2.497448	1.852335
77	1	2023-09-18 00:00:00+00:00	0.799407	-1.027642	2.626457

sf.plot(train, Y_hat.merge(test), max_insample_length=200)

使用置信区间的 predict 方法

使用 predict 方法生成预测。

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

h (int): 表示未来 h 步的预测。在此示例中，提前 30 天。
level (浮点数列表): 此可选参数用于概率预测。设置预测区间的水平（或置信百分位数）。例如，level=[95] 表示模型期望真实值有 95% 的时间落在此区间内。

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

此步骤应小于 1 秒钟。

sf.predict(h=horizon)

	unique_id	ds	GARCH(1,1)
0	1	2023-06-01 00:00:00+00:00	1.366914
1	1	2023-06-02 00:00:00+00:00	-0.593121
2	1	2023-06-05 00:00:00+00:00	-0.485200
…	…	…	…
75	1	2023-09-14 00:00:00+00:00	-1.686546
76	1	2023-09-15 00:00:00+00:00	-0.322556
77	1	2023-09-18 00:00:00+00:00	0.799407

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

	unique_id	ds	GARCH(1,1)	GARCH(1,1)-lo-95	GARCH(1,1)-lo-80	GARCH(1,1)-hi-80	GARCH(1,1)-hi-95
0	1	2023-06-01 00:00:00+00:00	1.366914	-0.021035	0.459383	2.274445	2.754863
1	1	2023-06-02 00:00:00+00:00	-0.593121	-2.435497	-1.797786	0.611543	1.249254
2	1	2023-06-05 00:00:00+00:00	-0.485200	-2.139216	-1.566703	0.596303	1.168815
…	…	…	…	…	…	…	…
7	1	2023-06-12 00:00:00+00:00	-1.051435	-4.790880	-3.496526	1.393657	2.688010
8	1	2023-06-13 00:00:00+00:00	0.421605	-3.001123	-1.816396	2.659607	3.844333
9	1	2023-06-14 00:00:00+00:00	-0.300086	-3.138338	-2.155920	1.555747	2.538166

sf.plot(train, test.merge(forecast_df), level=[80, 95], max_insample_length=200)

模型评估

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

from functools import partial

import utilsforecast.losses as ufl
from utilsforecast.evaluation import evaluate

evaluate(
    test.merge(Y_hat),
    metrics=[ufl.mae, ufl.mape, partial(ufl.mase, seasonality=season_length), ufl.rmse, ufl.smape],
    train_df=train,
)

	unique_id	metric	GARCH(1,1)
0	1	mae	0.843296
1	1	mape	3.703305
2	1	mase	0.794905
3	1	rmse	1.048076
4	1	smape	0.709150

参考文献

动态标准 Theta 模型 Holt 模型

在此页面上

目录
引言
GARCH 模型定义
广义自回归条件异方差 (GARCH) 模型的优缺点
广义自回归条件异方差 (GARCH) 模型可应用于多个领域
加载库和数据
读取数据
使用 plot 方法探索数据
增广迪基-福勒检验
收益序列
创建平方收益率
收益率 vs 平方收益率
Ljung-Box 检验
将数据拆分为训练集和测试集
使用 StatsForecast 实现 GARCH
加载库
实例化模型
交叉验证
执行时间序列交叉验证
拟合模型
预测方法
使用置信区间的 predict 方法
模型评估
参考资料

入门指南

教程

操作指南

分布式

实验

模型参考

API 参考

贡献

目录

引言

GARCH 模型定义

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

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

加载库和数据

读取数据

使用 plot 方法探索数据

增广迪基-福勒检验

收益序列

创建平方收益率

收益率 vs 平方收益率

Ljung-Box 检验

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

使用 StatsForecast 实现 GARCH

加载库

实例化模型

交叉验证

执行时间序列交叉验证

拟合模型

预测方法

使用置信区间的 predict 方法

模型评估

参考文献

入门指南

教程

操作指南

分布式

实验

模型参考

API 参考

贡献

​目录

​引言

​GARCH 模型定义

​广义自回归条件异方差 (GARCH) 模型的优缺点

​广义自回归条件异方差 (GARCH) 模型可应用于多个领域

​加载库和数据

​读取数据

​使用 plot 方法探索数据

​增广迪基-福勒检验

​收益序列

​创建平方收益率

​收益率 vs 平方收益率

​Ljung-Box 检验

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

​使用 StatsForecast 实现 GARCH

​加载库

​实例化模型

​交叉验证

​执行时间序列交叉验证

​拟合模型

​预测方法

​使用置信区间的 predict 方法

​模型评估

​参考文献

目录

引言

GARCH 模型定义

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

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

加载库和数据

读取数据

使用 plot 方法探索数据

增广迪基-福勒检验

收益序列

创建平方收益率

收益率 vs 平方收益率

Ljung-Box 检验

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

使用 StatsForecast 实现 GARCH

加载库

实例化模型

交叉验证

执行时间序列交叉验证

拟合模型

预测方法

使用置信区间的 predict 方法

模型评估

参考文献