Using GARCH to Forecast Stock Volatility and Assess Rolling Forecasts

# Is this the correct way to forecast stock price volatility using GARCH


# Is this the correct way to forecast stock price volatility using GARCH












I am attempting to make a forecast of a stock's volatility some time into the future (say 90 days). It seems that GARCH is a traditionally used model for this.

I have implemented this below using Python's `arch` library. Everything I do is explained in the comments, the only thing that needs to be changed to run the code is to provide your own daily prices, rather than where I retrieve them from my own API.

```
import utils
import numpy as np
import pandas as pd
import arch
import matplotlib.pyplot as plt

ticker = 'AAPL'         # Ticker to retrieve data for
forecast_horizon = 90   # Number of days to forecast

# Retrive prices from IEX API
prices = utils.dw.get(filename=ticker, source='iex', iex_range='5y')
df = prices[['date', 'close']]

df['daily_returns'] = np.log(df['close']).diff()            # Daily log returns
df['monthly_std'] = df['daily_returns'].rolling(21).std()   # Standard deviation across trading month
df['annual_vol'] = df['monthly_std'] * np.sqrt(252)         # Annualize monthly standard devation
df = df.dropna().reset_index(drop=True)

# Convert decimal returns to %
returns = df['daily_returns'] * 100

# Fit GARCH model
am = arch.arch_model(returns[:-forecast_horizon])
res = am.fit(disp='off')

# Calculate fitted variance values from model parameters
# Convert variance to standard deviation (volatility)
# Revert previous multiplication by 100
fitted = 0.1 * np.sqrt(
    res.params['omega'] +
    res.params['alpha[1]'] *
    res.resid**2 +
    res.conditional_volatility**2 *
    res.params['beta[1]']
)

# Make forecast
# Convert variance to standard deviation (volatility)
# Revert previous multiplication by 100
forecast = 0.1 * np.sqrt(res.forecast(horizon=forecast_horizon).variance.values[-1])

# Store actual, fitted, and forecasted results
vol = pd.DataFrame({
    'actual': df['annual_vol'],
    'model': np.append(fitted, forecast)
})

# Plot Actual vs Fitted/Forecasted
plt.plot(vol['actual'][:-forecast_horizon], label='Train')
plt.plot(vol['actual'][-forecast_horizon - 1:], label='Test')
plt.plot(vol['model'][:-forecast_horizon], label='Fitted')
plt.plot(vol['model'][-forecast_horizon - 1:], label='Forecast')
plt.legend()
plt.show()
```

For Apple, this produces the following plot:

Clearly, the fitted values are constantly far lower than the actual values, and this results in the forecast being a huge underestimation, too (This is a poor example given that Apple's volatility was unusually high in this test period, but with all companies I try, the model is always underestimating the fitted values).

Am I doing everything correct, and the GARCH model just isn't very powerful, or modelling volatility is very difficult? Or is there some error I am making?

## Answer by KOB (score 3, accepted)

https://quant.stackexchange.com/a/43713

I have solved this problem by performing a rolling forecast the following way. I am unsure (1) if this rolling forecast is correct, and (2) how to then perform a rolling forecast 30 days into the future.

```
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from rpy2.robjects.packages import importr
import rpy2.robjects as robjects
from rpy2.robjects import numpy2ri

ticker = 'AAPL'
forecast_horizon = 30

prices = utils.dw.get(filename=ticker, source='iex', iex_range='5y')
df = prices[['date', 'close']]

df['daily_returns'] = np.log(df['close']).diff()            # Daily log returns
df['monthly_std'] = df['daily_returns'].rolling(21).std()   # Standard deviation across trading month
df['annual_vol'] = df['monthly_std'] * np.sqrt(252)         # Convert monthly standard devation to annualized volatility
df = df.dropna().reset_index(drop=True)

# Initialize R GARCH model
rugarch = importr('rugarch')
garch_spec = rugarch.ugarchspec(
    mean_model=robjects.r('list(armaOrder = c(0,0))'),
    variance_model=robjects.r('list(garchOrder=c(1,1))'),
    distribution_model='std'
)

# Used to convert training set to R list for model input
numpy2ri.activate()

# Train R GARCH model on returns as %
garch_fitted = rugarch.ugarchfit(
    spec=garch_spec,
    data=df['daily_returns'].values * 100,
    out_sample=forecast_horizon
)

numpy2ri.deactivate()

# Model's fitted standard deviation values
# Revert previous multiplication by 100
# Convert to annualized volatility
fitted = 0.01 * np.sqrt(252) * np.array(garch_fitted.slots['fit'].rx2('sigma')).flatten()

# Forecast using R GACRH model
garch_forecast = rugarch.ugarchforecast(
    garch_fitted,
    n_ahead=1,
    n_roll=forecast_horizon - 1
)

# Model's forecasted standard deviation values
# Revert previous multiplication by 100
# Convert to annualized volatility
forecast = 0.01 * np.sqrt(252) * np.array(garch_forecast.slots['forecast'].rx2('sigmaFor')).flatten()

volatility = pd.DataFrame({
    'actual': df['annual_vol'].values,
    'model': np.append(fitted, forecast),
})

plt.plot(volatility['actual'][:-forecast_horizon], label='Train')
plt.plot(volatility['actual'][-forecast_horizon - 1:], label='Test')
plt.plot(volatility['model'][:-forecast_horizon], label='Fitted')
plt.plot(volatility['model'][-forecast_horizon - 1:], label='Forecasted')
plt.legend()
plt.show()
```

Which produces this plot: