波动率预测 (GARCH & ARCH)

在这个例子中，我们将使用 GARCH 和 ARCH 模型预测 S&P 500 以及几家上市公司的波动率。

先决条件

本教程假设您对StatsForecast有基本的了解。有关最小示例，请访问快速入门

介绍

广义自回归条件异方差（GARCH）模型用于具有随时间变化的非恒定波动性的时间序列中。这里的波动性指的是条件标准差。GARCH(p,q)模型由以下公式给出：

\[\begin{equation} y_t = v_t \sigma_t \end{equation}\]

其中 \(v_t\) 服从独立同分布，均值为零，方差为单位，而 \(\sigma_t\) 根据以下公式演变：

\[\begin{equation} \sigma_t^2 = w + \sum_{i=1}^p \alpha_i y^2_{t-i} + \sum_{j=1}^q \beta_j \sigma_{t-j}^2 \end{equation}\]

上述方程中的系数必须满足以下条件：

\(w>0\)，对于所有 \(i\) 的 \(\alpha_i \geq 0\)，对于所有 \(j\) 的 \(\beta_j \geq 0\)
\(\sum_{k=1}^{max(p,q)} \alpha_k + \beta_k < 1\)。这里假设对于 \(i>p\) 时 \(\alpha_i=0\)，对于 \(j>q\) 时 \(\beta_j=0\)。

GARCH模型的一个特例是ARCH模型，其中特殊情况为\(q=0\)。这两种模型在金融领域常用于建模股票价格、汇率、利率和其他金融工具的波动性。它们还被用于风险管理，以估计金融资产价格大幅波动的概率。

在本教程结束时，您将对如何在 StatsForecast 中实现GARCH或ARCH模型有一个良好的理解，并了解它们如何用于分析和预测金融时间序列数据。

大纲：

安装库
加载和探索数据
训练模型
执行时间序列交叉验证
评估结果
预测波动性

Tip

您可以使用 Colab 互动运行此笔记本

安装库

我们假设您已经安装了StatsForecast。如果没有，请查看本指南以获取如何安装StatsForecast 的说明。

安装必要的包，使用 pip install statsforecast

%%capture
pip install statsforecast -U

加载并探索数据

在本教程中，我们将使用过去5年的S&P 500和几家上市公司的价格数据。数据可以通过yfinance从Yahoo! Finance下载。要安装它，请使用pip install yfinance。

%%capture
pip install yfinance

我们还需要 pandas 来处理数据框。

import yfinance as yf
import pandas as pd 

tickers = ['SPY', 'MSFT', 'AAPL', 'GOOG', 'AMZN', 'TSLA', 'NVDA', 'META', 'NKE', 'NFLX'] 
df = yf.download(tickers, start = '2018-01-01', end = '2022-12-31', interval='1mo') # 使用月度价格
df.head()

[*********************100%***********************]  10 of 10 completed

	Adj Close										...	Volume
	AAPL	AMZN	GOOG	META	MSFT	NFLX	NKE	NVDA	SPY	TSLA	...	AAPL	AMZN	GOOG	META	MSFT	NFLX	NKE	NVDA	SPY	TSLA
Date
2018-01-01	39.741604	72.544502	58.497002	186.889999	89.248772	270.299988	64.929787	60.830006	258.821686	23.620667	...	2638717600	1927424000	574768000	495655700	574258400	238377600	157812200	1145621600	1985506700	1864072500
2018-02-01	42.279007	75.622498	55.236500	178.320007	88.083969	291.380005	63.797192	59.889591	249.410812	22.870667	...	3711577200	2755680000	847640000	516251600	725663300	184585800	160317000	1491552800	2923722000	1637850000
2018-03-01	39.987053	72.366997	51.589500	159.789993	86.138298	295.350006	63.235649	57.348976	241.606750	17.742001	...	2854910800	2608002000	907066000	996201700	750754800	263449400	174066700	1411844000	2323561800	2359027500
2018-04-01	39.386456	78.306503	50.866501	172.000000	88.261810	312.459991	65.288467	55.692314	243.828018	19.593332	...	2664617200	2598392000	834318000	750072700	668130700	262006000	158981900	1114400800	1998466500	2854662000
2018-05-01	44.536777	81.481003	54.249500	191.779999	93.282692	351.600006	68.543846	62.450180	249.755264	18.982000	...	2483905200	1432310000	636988000	401144100	509417900	142050800	129566300	1197824000	1606397200	2333671500

5 rows × 60 columns

下载的数据包含不同的价格。我们将使用调整后的收盘价，即在考虑到股票分拆或股息分配等公司行动后得出的收盘价。它也是用于研究历史回报的价格。

请注意，yfinance 返回的 dataframe 具有多重索引，因此我们需要同时选择调整后的价格和股票代码。

df = df.loc[:, (['Adj Close'], tickers)]
df.columns = df.columns.droplevel() # 删除多重索引
df = df.reset_index()
df.head()

	Date	SPY	MSFT	AAPL	GOOG	AMZN	TSLA	NVDA	META	NKE	NFLX
0	2018-01-01	258.821686	89.248772	39.741604	58.497002	72.544502	23.620667	60.830006	186.889999	64.929787	270.299988
1	2018-02-01	249.410812	88.083969	42.279007	55.236500	75.622498	22.870667	59.889591	178.320007	63.797192	291.380005
2	2018-03-01	241.606750	86.138298	39.987053	51.589500	72.366997	17.742001	57.348976	159.789993	63.235649	295.350006
3	2018-04-01	243.828018	88.261810	39.386456	50.866501	78.306503	19.593332	55.692314	172.000000	65.288467	312.459991
4	2018-05-01	249.755264	93.282692	44.536777	54.249500	81.481003	18.982000	62.450180	191.779999	68.543846	351.600006

StatsForecast的输入是一个长格式的数据框，包含三列：unique_id、ds和y：

unique_id：（字符串、整数或类别）系列的唯一标识符。
ds：（日期戳或整数）格式为YYYY-MM-DD或YYYY-MM-DD HH:MM:SS的日期戳，或一个用于索引时间的整数。
y：（数值型）我们希望进行预测的测量值。

因此，我们需要重新构造数据。我们将通过创建一个新的数据框来实现这一点，命名为price。

prices = df.melt(id_vars = 'Date')
prices = prices.rename(columns={'Date': 'ds', 'variable': 'unique_id', 'value': 'y'})
prices = prices[['unique_id', 'ds', 'y']]
prices

	unique_id	ds	y
0	SPY	2018-01-01	258.821686
1	SPY	2018-02-01	249.410812
2	SPY	2018-03-01	241.606750
3	SPY	2018-04-01	243.828018
4	SPY	2018-05-01	249.755264
...	...	...	...
595	NFLX	2022-08-01	223.559998
596	NFLX	2022-09-01	235.440002
597	NFLX	2022-10-01	291.880005
598	NFLX	2022-11-01	305.529999
599	NFLX	2022-12-01	294.880005

600 rows × 3 columns

我们可以使用 StatsForecast 类的 plot 方法绘制该系列。

from statsforecast import StatsForecast

/home/ubuntu/statsforecast/statsforecast/core.py:21: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)
  from tqdm.autonotebook import tqdm

StatsForecast.plot(prices)

通过价格，我们可以计算标普500和上市公司的对数收益。这是我们感兴趣的变量，因为它可能与GARCH框架相结合得很好。对数收益由以下公式给出：

\(return_t = log \big( \frac{price_t}{price_{t-1}} \big)\)

我们将在价格数据框上计算收益，然后将使用StatsForecast的格式创建一个收益数据框。为此，我们需要numpy。

import numpy as np 
prices['rt'] = prices['y'].div(prices.groupby('unique_id')['y'].shift(1))
prices['rt'] = np.log(prices['rt'])

returns = prices[['unique_id', 'ds', 'rt']]
returns = returns.rename(columns={'rt':'y'})
returns

	unique_id	ds	y
0	SPY	2018-01-01	NaN
1	SPY	2018-02-01	-0.037038
2	SPY	2018-03-01	-0.031790
3	SPY	2018-04-01	0.009152
4	SPY	2018-05-01	0.024018
...	...	...	...
595	NFLX	2022-08-01	-0.005976
596	NFLX	2022-09-01	0.051776
597	NFLX	2022-10-01	0.214887
598	NFLX	2022-11-01	0.045705
599	NFLX	2022-12-01	-0.035479

600 rows × 3 columns

Warning

如果数据的秩序非常小（例如 \(<1e-5\)），scipy.optimize.minimize 可能无法成功终止。在这种情况下，请对数据进行重缩放，然后生成 GARCH 或 ARCH 模型。

StatsForecast.plot(returns)

从这个图中，我们可以看到，收益似乎适合GARCH框架，因为大的冲击_往往_会随之而来其他大的冲击。这并不意味着在每次大的冲击之后我们都应期待另一次；只是大方差的概率大于小方差的概率。

训练模型

首先，我们需要从 statsforecast.models 导入 GARCH 和 ARCH 模型，然后通过实例化一个新的 StatsForecast 对象来拟合它们。请注意，我们将使用不同的 \(p\) 和 \(q\) 值。在接下来的部分中，我们将通过交叉验证确定哪些值产生最准确的模型。我们还将导入 Naive 模型，因为我们将其用作基准。

from statsforecast.models import (
    GARCH, 
    ARCH, 
    Naive
)

models = [ARCH(1), 
          ARCH(2), 
          GARCH(1,1),
          GARCH(1,2),
          GARCH(2,2),
          GARCH(2,1),
          Naive()
]

要实例化一个新的StatsForecast对象，我们需要以下参数：

df：包含训练数据的数据框。
models：在上一步中定义的模型列表。
freq：一个字符串，表示数据的频率。在这里我们将使用MS，即代表每个月的开始。您可以在此处查看pandas可用频率的列表。
n_jobs：一个整数，表示并行处理所使用的作业数量。使用-1以选择所有核心。

sf = StatsForecast(
    df = returns, 
    models = models, 
    freq = 'MS',
    n_jobs = -1
)

执行时间序列交叉验证

时间序列交叉验证是一种评估模型在过去表现的方法。它通过在历史数据上定义一个滑动窗口，并预测之后的时间段来实现。这里我们将使用StatsForecast的cross-validation方法来确定S&P 500和所选公司的最准确模型。

该方法接受以下参数：

df：包含训练数据的数据框。
h（整数）：表示未来h步的预测。
step_size（整数）：每个窗口之间的步长，意味着您希望多频繁运行预测过程。
n_windows（整数）：用于交叉验证的窗口数量，意味着您希望评估过去的预测过程数量。

对于这个特定的例子，我们将使用4个窗口，每个窗口3个月，或者一年中的所有季度。

crossvalidation_df = sf.cross_validation(
    df = returns,
    h = 3,
    step_size = 3,
    n_windows = 4
  )

crossvalidation_df 对象是一个包含以下列的数据框：

unique_id: 索引。
ds: 日期戳或时间索引。
cutoff: n_windows 的最后一个日期戳或时间索引。
y: 真实值。
"model": 包含模型名称和拟合值的列。

crossvalidation_df = crossvalidation_df.reset_index()
crossvalidation_df.rename(columns = {'y' : 'actual'}, inplace = True)
crossvalidation_df.head()

	unique_id	ds	cutoff	actual	ARCH(1)	ARCH(2)	GARCH(1,1)	GARCH(1,2)	GARCH(2,2)	GARCH(2,1)	Naive
0	AAPL	2022-01-01	2021-12-01	-0.015837	0.142416	0.144013	0.142951	0.226098	0.141690	0.144018	0.073061
1	AAPL	2022-02-01	2021-12-01	-0.056855	-0.056896	-0.057158	-0.056387	-0.087001	-0.058787	-0.057161	0.073061
2	AAPL	2022-03-01	2021-12-01	0.057156	-0.045899	-0.046478	-0.047512	-0.073625	-0.045714	-0.046479	0.073061
3	AAPL	2022-04-01	2022-03-01	-0.102178	0.138661	0.140211	0.136213	0.136124	0.136127	0.136546	0.057156
4	AAPL	2022-05-01	2022-03-01	-0.057505	-0.056013	-0.056268	-0.054599	-0.057080	-0.057085	-0.053791	0.057156

StatsForecast.plot(returns, crossvalidation_df.drop(['cutoff', 'actual'], axis=1))

关于交叉验证的教程可以在这里找到。

评估结果

为了计算预测的准确性，我们将使用平均绝对误差（mae），它是绝对误差的总和除以预测的数量。我们在datasetsforecast上有MAE的实现，因此我们将安装它，然后导入mae函数。

%%capture
pip install datasetsforecast -U

from datasetsforecast.losses import mae

需要计算每个窗口的MAE，然后对所有窗口的结果进行平均。为此，我们将创建以下函数。

def compute_cv_mae(crossvalidation_df):
    """计算所有生成模型的平均绝对误差（MAE）"""
    res = {}
    for mod in models: 
        res[mod] = mae(crossvalidation_df['actual'], crossvalidation_df[str(mod)])
    return pd.Series(res)

mae_cv = crossvalidation_df.groupby(['unique_id', 'cutoff']).apply(compute_cv_mae)

mae = mae_cv.groupby('unique_id').mean()
mae.style.highlight_min(color = 'lightblue', axis = 1)

	ARCH(1)	ARCH(2)	GARCH(1,1)	GARCH(1,2)	GARCH(2,2)	GARCH(2,1)	Naive
unique_id
AAPL	0.068537	0.068927	0.068929	0.085630	0.072519	0.068556	0.110426
AMZN	0.118612	0.126182	0.118858	0.125470	0.109913	0.109912	0.115189
GOOG	0.093849	0.093752	0.099593	0.115136	0.094648	0.113645	0.083233
META	0.198333	0.198891	0.199617	0.199712	0.199708	0.198890	0.185346
MSFT	0.080022	0.097301	0.082183	0.072765	0.073006	0.080494	0.086951
NFLX	0.159384	0.159523	0.219658	0.231798	0.230077	0.224103	0.167421
NKE	0.107842	0.114263	0.103097	0.107180	0.107179	0.107019	0.160405
NVDA	0.189462	0.207875	0.199004	0.196172	0.211928	0.211928	0.215289
SPY	0.058513	0.065498	0.058700	0.057051	0.057051	0.058526	0.089012
TSLA	0.192003	0.192620	0.190225	0.192353	0.191620	0.191418	0.218857

因此，描述苹果公司股票对数收益的最准确模型是ARCH(1)，而亚马逊股票的模型是GARCH(2,1)，等等。

预测波动性

我们现在可以生成下一个季度的预测。为此，我们将使用forecast方法，该方法需要以下参数：

h： (int) 预测范围。
level： (list[float]) 预测区间的置信水平。
fitted： (bool = False) 返回样本内预测值。

levels = [80, 95] # 预测区间的置信水平 

forecasts = sf.forecast(h=3, level=levels)
forecasts = forecasts.reset_index()
forecasts.head()

	unique_id	ds	ARCH(1)	ARCH(1)-lo-95	ARCH(1)-lo-80	ARCH(1)-hi-80	ARCH(1)-hi-95	ARCH(2)	ARCH(2)-lo-95	ARCH(2)-lo-80	...	GARCH(2,1)	GARCH(2,1)-lo-95	GARCH(2,1)-lo-80	GARCH(2,1)-hi-80	GARCH(2,1)-hi-95	Naive	Naive-lo-80	Naive-lo-95	Naive-hi-80	Naive-hi-95
0	AAPL	2023-01-01	0.150457	0.133641	0.139462	0.161452	0.167273	0.150166	0.133415	0.139213	...	0.147610	0.131424	0.137027	0.158193	0.163795	-0.128762	-0.284463	-0.366886	0.026939	0.109362
1	AAPL	2023-02-01	-0.056942	-0.073923	-0.068046	-0.045839	-0.039961	-0.057209	-0.074349	-0.068417	...	-0.059511	-0.078059	-0.071639	-0.047384	-0.040964	-0.128762	-0.348956	-0.465520	0.091433	0.207997
2	AAPL	2023-03-01	-0.048390	-0.064842	-0.059148	-0.037633	-0.031939	-0.049279	-0.066340	-0.060435	...	-0.054537	-0.075435	-0.068201	-0.040874	-0.033640	-0.128762	-0.398444	-0.541205	0.140920	0.283681
3	AMZN	2023-01-01	0.152158	0.134960	0.140913	0.163404	0.169357	0.148659	0.132243	0.137925	...	0.148597	0.132195	0.137872	0.159322	0.165000	-0.139141	-0.315716	-0.409190	0.037435	0.130909
4	AMZN	2023-02-01	-0.057306	-0.074504	-0.068551	-0.046060	-0.040107	-0.061187	-0.080794	-0.074007	...	-0.069302	-0.094455	-0.085749	-0.052856	-0.044150	-0.139141	-0.388856	-0.521048	0.110575	0.242767

5 rows × 37 columns

根据上一部分的结果，我们可以为标准普尔500指数及所选公司选择最佳模型。下面显示了一些图表。请注意，我们在plot方法中使用了一些额外的参数：

level: （list[int]）预测区间的置信水平（这已经定义过）。
unique_ids: （list[str, int 或类别]）要绘制的 ID。
models: （list(str)）。要绘制的模型。在这种情况下，是通过交叉验证选择的模型。

StatsForecast.plot(returns, forecasts, level=levels, unique_ids = ['AAPL'], models = ['GARCH(2,1)'])

StatsForecast.plot(returns, forecasts, level=levels, unique_ids = ['MSFT'], models = ['ARCH(2)'])

StatsForecast.plot(returns, forecasts, level=levels, unique_ids = ['NFLX'], models = ['ARCH(1)'])

参考文献

Engle, R. F. (1982). 自回归条件异方差及对英国通货膨胀方差的估计。经济计量学：经济计量学会期刊，987-1007。
Bollerslev, T. (1986). 广义自回归条件异方差。计量经济学杂志，31(3), 307-327。
Hamilton, J. D. (1994). 时间序列分析。普林斯顿大学出版社。
Tsay, R. S. (2005). 金融时间序列分析。约翰·威利父子公司。

Give us a ⭐ on Github