Prerequesites

This tutorial assumes basic familiarity with StatsForecast. For a minimal example visit the Quick Start

Introduction

The Generalized Autoregressive Conditional Heteroskedasticity (GARCH) model is used for time series that exhibit non-constant volatility over time. Here volatility refers to the conditional standard deviation. The GARCH(p,q) model is given by

where vtv_t is independent and identically distributed with zero mean and unit variance, and σt\sigma_t evolves according to

The coefficients in the equation above must satisfy the following conditions:

  1. w>0w>0, αi0\alpha_i \geq 0 for all ii, and βj0\beta_j \geq 0 for all jj
  2. k=1max(p,q)αk+βk<1\sum_{k=1}^{max(p,q)} \alpha_k + \beta_k < 1. Here it is assumed that αi=0\alpha_i=0 for i>pi>p and βj=0\beta_j=0 for j>qj>q.

A particular case of the GARCH model is the ARCH model, in which q=0q=0. Both models are commonly used in finance to model the volatility of stock prices, exchange rates, interest rates, and other financial instruments. They’re also used in risk management to estimate the probability of large variations in the price of financial assets.

By the end of this tutorial, you’ll have a good understanding of how to implement a GARCH or an ARCH model in StatsForecast and how they can be used to analyze and predict financial time series data.

Outline:

  1. Install libraries
  2. Load and explore the data
  3. Train models
  4. Perform time series cross-validation
  5. Evaluate results
  6. Forecast volatility

Tip

You can use Colab to run this Notebook interactively

Install libraries

We assume that you have StatsForecast already installed. If not, check this guide for instructions on how to install StatsForecast

Install the necessary packages using pip install statsforecast

pip install statsforecast -U

Load and explore the data

In this tutorial, we’ll use the last 5 years of prices from the S&P 500 and several publicly traded companies. The data can be downloaded from Yahoo! Finance using yfinance. To install it, use pip install yfinance.

pip install yfinance

We’ll also need pandas to deal with the dataframes.

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') # use monthly prices
df.head()
[*********************100%***********************]  10 of 10 completed
Adj CloseVolume
AAPLAMZNGOOGMETAMSFTNFLXNKENVDASPYTSLAAAPLAMZNGOOGMETAMSFTNFLXNKENVDASPYTSLA
Date
2018-01-0139.74160472.54450258.497002186.88999989.248772270.29998864.92978760.830006258.82168623.62066726387176001927424000574768000495655700574258400238377600157812200114562160019855067001864072500
2018-02-0142.27900775.62249855.236500178.32000788.083969291.38000563.79719259.889591249.41081222.87066737115772002755680000847640000516251600725663300184585800160317000149155280029237220001637850000
2018-03-0139.98705372.36699751.589500159.78999386.138298295.35000663.23564957.348976241.60675017.74200128549108002608002000907066000996201700750754800263449400174066700141184400023235618002359027500
2018-04-0139.38645678.30650350.866501172.00000088.261810312.45999165.28846755.692314243.82801819.59333226646172002598392000834318000750072700668130700262006000158981900111440080019984665002854662000
2018-05-0144.53677781.48100354.249500191.77999993.282692351.60000668.54384662.450180249.75526418.98200024839052001432310000636988000401144100509417900142050800129566300119782400016063972002333671500

The data downloaded includes different prices. We’ll use the adjusted closing price, which is the closing price after accounting for any corporate actions like stock splits or dividend distributions. It is also the price that is used to examine historical returns.

Notice that the dataframe that yfinance returns has a MultiIndex, so we need to select both the adjusted price and the tickers.

df = df.loc[:, (['Adj Close'], tickers)]
df.columns = df.columns.droplevel() # drop MultiIndex
df = df.reset_index()
df.head()
DateSPYMSFTAAPLGOOGAMZNTSLANVDAMETANKENFLX
02018-01-01258.82168689.24877239.74160458.49700272.54450223.62066760.830006186.88999964.929787270.299988
12018-02-01249.41081288.08396942.27900755.23650075.62249822.87066759.889591178.32000763.797192291.380005
22018-03-01241.60675086.13829839.98705351.58950072.36699717.74200157.348976159.78999363.235649295.350006
32018-04-01243.82801888.26181039.38645650.86650178.30650319.59333255.692314172.00000065.288467312.459991
42018-05-01249.75526493.28269244.53677754.24950081.48100318.98200062.450180191.77999968.543846351.600006

The input to StatsForecast is a dataframe in long format with three columns: unique_id, ds and y:

  • unique_id: (string, int or category) A unique identifier for the series.
  • ds: (datestamp or int) A datestamp in format YYYY-MM-DD or YYYY-MM-DD HH:MM:SS or an integer indexing time.
  • y: (numeric) The measurement we wish to forecast.

Hence, we need to reshape the data. We’ll do this by creating a new dataframe called 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_iddsy
0SPY2018-01-01258.821686
1SPY2018-02-01249.410812
2SPY2018-03-01241.606750
3SPY2018-04-01243.828018
4SPY2018-05-01249.755264
595NFLX2022-08-01223.559998
596NFLX2022-09-01235.440002
597NFLX2022-10-01291.880005
598NFLX2022-11-01305.529999
599NFLX2022-12-01294.880005

We can plot this series using the plot method of the StatsForecast class.

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)

With the prices, we can compute the logarithmic returns of the S&P 500 and the publicly traded companies. This is the variable we’re interested in since it’s likely to work well with the GARCH framework. The logarithmic return is given by

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

We’ll compute the returns on the price dataframe and then we’ll create a return dataframe with StatsForecast’s format. To do this, we’ll need 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_iddsy
0SPY2018-01-01NaN
1SPY2018-02-01-0.037038
2SPY2018-03-01-0.031790
3SPY2018-04-010.009152
4SPY2018-05-010.024018
595NFLX2022-08-01-0.005976
596NFLX2022-09-010.051776
597NFLX2022-10-010.214887
598NFLX2022-11-010.045705
599NFLX2022-12-01-0.035479

Warning

If the order of the data is very small (say <1e5<1e-5), scipy.optimize.minimize might not terminate successfully. In this case, rescale the data and then generate the GARCH or ARCH model.

StatsForecast.plot(returns)

From this plot, we can see that the returns seem suited for the GARCH framework, since large shocks tend to be followed by other large shocks. This doesn’t mean that after every large shock we should expect another one; merely that the probability of a large variance is greater than the probability of a small one.

Train models

We first need to import the GARCH and the ARCH models from statsforecast.models, and then we need to fit them by instantiating a new StatsForecast object. Notice that we’ll be using different values of pp and qq. In the next section, we’ll determine which ones produce the most accurate model using cross-validation. We’ll also import the Naive model since we’ll use it as a baseline.

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()
]

To instantiate a new StatsForecast object, we need the following parameters:

  • df: The dataframe with the training data.
  • models: The list of models defined in the previous step.
  • freq: A string indicating the frequency of the data. Here we’ll use MS, which correspond to the start of the month. You can see the list of panda’s available frequencies here.
  • n_jobs: An integer that indicates the number of jobs used in parallel processing. Use -1 to select all cores.
sf = StatsForecast(
    df = returns, 
    models = models, 
    freq = 'MS',
    n_jobs = -1
)

Perform time series cross-validation

Time series cross-validation is a method for evaluating how a model would have performed in the past. It works by defining a sliding window across the historical data and predicting the period following it. Here we’ll use StatsForercast’s cross-validation method to determine the most accurate model for the S&P 500 and the companies selected.

This method takes the following arguments:

  • df: The dataframe with the training data.
  • h (int): represents the h steps into the future that will be forecasted.
  • step_size (int): step size between each window, meaning how often do you want to run the forecasting process.
  • n_windows (int): number of windows used for cross-validation, meaning the number of forecasting processes in the past you want to evaluate.

For this particular example, we’ll use 4 windows of 3 months, or all the quarters in a year.

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

The crossvalidation_df object ia a dataframe with the following columns:

  • unique_id: index.
  • ds: datestamp or temporal index
  • cutoff: the last datestamp or temporal index for the n_windows.
  • y: true value
  • "model": columns with the model’s name and fitted value.
crossvalidation_df = crossvalidation_df.reset_index()
crossvalidation_df.rename(columns = {'y' : 'actual'}, inplace = True)
crossvalidation_df.head()
unique_iddscutoffactualARCH(1)ARCH(2)GARCH(1,1)GARCH(1,2)GARCH(2,2)GARCH(2,1)Naive
0AAPL2022-01-012021-12-01-0.0158370.1424160.1440130.1429510.2260980.1416900.1440180.073061
1AAPL2022-02-012021-12-01-0.056855-0.056896-0.057158-0.056387-0.087001-0.058787-0.0571610.073061
2AAPL2022-03-012021-12-010.057156-0.045899-0.046478-0.047512-0.073625-0.045714-0.0464790.073061
3AAPL2022-04-012022-03-01-0.1021780.1386610.1402110.1362130.1361240.1361270.1365460.057156
4AAPL2022-05-012022-03-01-0.057505-0.056013-0.056268-0.054599-0.057080-0.057085-0.0537910.057156
StatsForecast.plot(returns, crossvalidation_df.drop(['cutoff', 'actual'], axis=1))

A tutorial on cross-validation can be found here.

Evaluate results

To compute the accuracy of the forecasts, we’ll use the mean average error (mae), which is the sum of the absolute errors divided by the number of forecasts. There’s an implementation of MAE on datasetsforecast, so we’ll install it and then import the mae function.

pip install datasetsforecast -U
from datasetsforecast.losses import mae

The MAE needs to be computed for every window and then it needs to be averaged across all of them. To do this, we’ll create the following function.

def compute_cv_mae(crossvalidation_df):
    """Compute MAE for all models generated"""
    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       
AAPL0.0685370.0689270.0689290.0856300.0725190.0685560.110426
AMZN0.1186120.1261820.1188580.1254700.1099130.1099120.115189
GOOG0.0938490.0937520.0995930.1151360.0946480.1136450.083233
META0.1983330.1988910.1996170.1997120.1997080.1988900.185346
MSFT0.0800220.0973010.0821830.0727650.0730060.0804940.086951
NFLX0.1593840.1595230.2196580.2317980.2300770.2241030.167421
NKE0.1078420.1142630.1030970.1071800.1071790.1070190.160405
NVDA0.1894620.2078750.1990040.1961720.2119280.2119280.215289
SPY0.0585130.0654980.0587000.0570510.0570510.0585260.089012
TSLA0.1920030.1926200.1902250.1923530.1916200.1914180.218857

Hence, the most accurate model to describe the logarithmic returns of Apple’s stock is an ARCH(1), for Amazon’s stock is a GARCH(2,1), and so on.

Forecast volatility

We can now generate a forecast for the next quarter. To do this, we’ll use the forecast method, which requieres the following arguments:

  • h: (int) The forecasting horizon.
  • level: (list[float]) The confidence levels of the prediction intervals
  • fitted : (bool = False) Returns insample predictions.
levels = [80, 95] # confidence levels for the prediction intervals 

forecasts = sf.forecast(h=3, level=levels)
forecasts = forecasts.reset_index()
forecasts.head()
unique_iddsARCH(1)ARCH(1)-lo-95ARCH(1)-lo-80ARCH(1)-hi-80ARCH(1)-hi-95ARCH(2)ARCH(2)-lo-95ARCH(2)-lo-80GARCH(2,1)GARCH(2,1)-lo-95GARCH(2,1)-lo-80GARCH(2,1)-hi-80GARCH(2,1)-hi-95NaiveNaive-lo-80Naive-lo-95Naive-hi-80Naive-hi-95
0AAPL2023-01-010.1504570.1336410.1394620.1614520.1672730.1501660.1334150.1392130.1476100.1314240.1370270.1581930.163795-0.128762-0.284463-0.3668860.0269390.109362
1AAPL2023-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.4655200.0914330.207997
2AAPL2023-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.5412050.1409200.283681
3AMZN2023-01-010.1521580.1349600.1409130.1634040.1693570.1486590.1322430.1379250.1485970.1321950.1378720.1593220.165000-0.139141-0.315716-0.4091900.0374350.130909
4AMZN2023-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.5210480.1105750.242767

With the results of the previous section, we can choose the best model for the S&P 500 and the companies selected. Some of the plots are shown below. Notice that we’re using somo additional arguments in the plot method:

  • level: (list[int]) The confidence levels for the prediction intervals (this was already defined).
  • unique_ids: (list[str, int or category]) The ids to plot.
  • models: (list(str)). The model to plot. In this case, is the model selected by cross-validation.
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)'])

References