I had a difficult time to understand GARCH and its variants. In this post, I am going to show you what I have come across while learning and experimenting on this topic. If you are well-versed in this area, please do keep reading and point out the mistakes in this piece. Hopefully, it will help someone in the future.

We are going to use BTCUSD as an example for the rest of this article. Feel free to try the following techniques on other cryptocurrencies pairs. Shoot me a message on Patreon, Discord or Twitter and let me know what you have found.

You might also like to sign up for my newsletter and stay tuned with the latest articles.

Data

First of all, let us get BTCUSD minute data. Due to the sheer size of the data, you can scale down the time window if the code takes too long.

from catalyst.api import symbol , record from catalyst import run_algorithm import pandas as pd trading_pair = 'btc_usd' frequency = 'minute' exchange = 'bitfinex' start = '2017-4-1' end = '2019-4-1' capital_base = 1000 quote_currency = trading_pair . split ( '_' )[ 1 ] def initialize ( context ): context . asset = symbol ( trading_pair ) def handle_data ( context , data ): # The last known price and volume of current date/minute and the day/minute before if frequency == 'daily' : price = data . current ( context . asset , 'price' ) volume = data . current ( context . asset , 'volume' ) elif frequency == 'minute' : price = data . current ( context . asset , 'price' ) volume = data . current ( context . asset , 'volume' ) record ( price = price , volume = volume ) if __name__ == '__main__' : perf = run_algorithm ( capital_base = capital_base , data_frequency = frequency , initialize = initialize , handle_data = handle_data , exchange_name = exchange , quote_currency = quote_currency , start = pd . to_datetime ( start , utc = True ), end = pd . to_datetime ( end , utc = True ))

Once it is done, we will have 2-year of BTCUSD minute data from Bitfinex.

Now, it is time to resample our data and perform a few simple calculations.

import numpy as np # Resample to minute data m_df = pd . DataFrame ( perf . loc [:, [ 'price' , 'volume' ]], index = perf . index ) . resample ( 'T' , closed = 'left' , label = 'left' ) . mean () . copy () m_df . index . name = 'timestamp' m_df [ 'log_price' ] = np . log ( m_df . price ) m_df [ 'return' ] = m_df . price . pct_change () . dropna () m_df [ 'log_return' ] = m_df [ 'log_price' ] - m_df [ 'log_price' ] . shift ( 1 ) # Calculate squared log return m_df [ 'squared_log_return' ] = np . power ( m_df [ 'log_return' ], 2 ) # Scale up 100x m_df [ 'return_100x' ] = np . multiply ( m_df [ 'return' ], 100 ) m_df [ 'log_return_100x' ] = np . multiply ( m_df [ 'log_return' ], 100 ) m_df . head ()

price volume log_price return log_return squared_log_return return_100x log_return_100x timestamp 2017-04-01 00:01:00+00:00 1081.5 6.752945 6.986104 NaN NaN NaN NaN NaN 2017-04-01 00:02:00+00:00 1081.4 3.928597 6.986012 -0.000092 -0.000092 8.550413e-09 -0.009246 -0.009247 2017-04-01 00:03:00+00:00 1081.7 2.540370 6.986289 0.000277 0.000277 7.693949e-08 0.027742 0.027738 2017-04-01 00:04:00+00:00 1081.6 16.532780 6.986197 -0.000092 -0.000092 8.547253e-09 -0.009245 -0.009245 2017-04-01 00:05:00+00:00 1082.2 15.992445 6.986751 0.000555 0.000555 3.075589e-07 0.055473 0.055458

You may wonder why we need scale up 100 times for simple returns and log returns. This is because ARCH Python package might run into convergence issues. Do not worry, and I will show you another way around without scaling up.

Volatility

Since there is a general consense that the volatility is more sticky than the price in financial market, much of the effort has been allocated to study the models of volatility. Before we can jump into any fancy models and math, let us define what is volatility and how to caculate it.

The following is cited from Wikipedia.

In finance, volatility (symbol σ) is the degree of variation of a trading price series over time as measured by the standard deviation of logarithmic returns.

We can arrive at the mathematical definiton with ease.

Where,

$\sigma$ is volatility,

$\bar{r}$ is mean log return,

$r_i$ is the log return at time i,

$n$ is the number of log return observations.

Daily Volatility & Variance

A simple volatility calculated from a range of the sample time series is not that useful when we need to implement more advanced models. We are more inclined to estimate the volatility/variance for a certain time interval.

Equally as important as the forecast methodology is the proxy of variance that one measures performance against.

One common way is to use daily squared log return $r_d^2$ as a proxy of eastimating the true variance. Let us resample the data by day and calculate daily squared log returns.

# Resample to daily data d_df = pd . DataFrame ( m_df . loc [:, [ 'price' , 'volume' ]], index = m_df . index ) . resample ( 'D' , closed = 'left' , label = 'left' ) . mean () . copy () d_df [ 'log_price' ] = np . log ( d_df . price ) d_df [ 'return' ] = d_df . price . pct_change () . dropna () d_df [ 'log_return' ] = d_df [ 'log_price' ] - d_df [ 'log_price' ] . shift ( 1 ) d_df [ 'squared_log_return' ] = np . power ( d_df [ 'log_return' ], 2 ) # Scale up 100x d_df [ 'return_100x' ] = np . multiply ( d_df [ 'return' ], 100 ) d_df [ 'log_return_100x' ] = np . multiply ( d_df [ 'log_return' ], 100 ) d_df . head ()

price volume log_price return log_return squared_log_return return_100x log_return_100x timestamp 2017-04-01 00:00:00+00:00 1086.633148 11.808373 6.990839 NaN NaN NaN NaN NaN 2017-04-02 00:00:00+00:00 1095.521806 14.183194 6.998986 0.008180 0.008147 0.000066 0.818000 0.814672 2017-04-03 00:00:00+00:00 1138.100972 21.106648 7.037116 0.038867 0.038130 0.001454 3.886656 3.813028 2017-04-04 00:00:00+00:00 1143.858542 14.609319 7.042163 0.005059 0.005046 0.000025 0.505893 0.504617 2017-04-05 00:00:00+00:00 1130.486875 15.630853 7.030404 -0.011690 -0.011759 0.000138 -1.168997 -1.175883

However, the downside of this method is incredibly noisy.

Another commonly accepted approach is called realized volatility (RV), coined by Bollershev. Basically, the summation of squared intraday returns aggregated to the daily level can be used to measure that day’s variance. For instance, if we have BTCUSD minute data, its daily variance can be estimated as follows:

Here is how we implement it in Python.

d_df [ 'realized_variance_1min' ] = pd . Series ( m_df . loc [:, 'squared_log_return' ], index = perf . index ) . resample ( 'D' , closed = 'left' , label = 'left' ) . sum () . copy () d_df [ 'realized_volatility_1min' ] = np . sqrt ( d_df [ 'realized_variance_1min' ]) d_df . head ()

price volume log_price return log_return squared_log_return return_100x log_return_100x realized_variance_1min realized_volatility_1min timestamp 2017-04-01 00:00:00+00:00 1086.633148 11.808373 6.990839 NaN NaN NaN NaN NaN 0.000645 0.025388 2017-04-02 00:00:00+00:00 1095.521806 14.183194 6.998986 0.008180 0.008147 0.000066 0.818000 0.814672 0.000757 0.027519 2017-04-03 00:00:00+00:00 1138.100972 21.106648 7.037116 0.038867 0.038130 0.001454 3.886656 3.813028 0.000751 0.027406 2017-04-04 00:00:00+00:00 1143.858542 14.609319 7.042163 0.005059 0.005046 0.000025 0.505893 0.504617 0.000510 0.022578 2017-04-05 00:00:00+00:00 1130.486875 15.630853 7.030404 -0.011690 -0.011759 0.000138 -1.168997 -1.175883 0.000493 0.022203

Another common interval to calculate RV is 5-min data interval. Of course, the formula mentioned above will be re-written like this.

# Resample to 5-min data five_min_df = pd . DataFrame ( m_df . loc [:, [ 'price' , 'volume' ]], index = m_df . index ) . resample ( '5T' , closed = 'left' , label = 'left' ) . mean () . copy () five_min_df [ 'log_price' ] = np . log ( five_min_df . price ) five_min_df [ 'log_return' ] = five_min_df [ 'log_price' ] - five_min_df [ 'log_price' ] . shift ( 1 ) five_min_df [ 'squared_log_return' ] = np . power ( five_min_df [ 'log_return' ], 2 ) five_min_df . head ()

price volume log_price log_return squared_log_return timestamp 2017-04-01 00:00:00+00:00 1081.55 7.438673 6.986150 NaN NaN 2017-04-01 00:05:00+00:00 1081.98 19.154190 6.986548 0.000397 1.580051e-07 2017-04-01 00:10:00+00:00 1082.96 2.362552 6.987453 0.000905 8.196350e-07 2017-04-01 00:15:00+00:00 1083.78 2.617299 6.988210 0.000757 5.728938e-07 2017-04-01 00:20:00+00:00 1083.86 2.761562 6.988284 0.000074 5.448358e-09

d_df [ 'realized_variance_5min' ] = pd . Series ( five_min_df . loc [:, 'squared_log_return' ], index = perf . index ) . resample ( 'D' , closed = 'left' , label = 'left' ) . sum () . copy () d_df [ 'realized_volatility_5min' ] = np . sqrt ( d_df [ 'realized_variance_5min' ]) d_df . head ()

price volume log_price return log_return squared_log_return return_100x log_return_100x realized_variance_5min realized_volatility_5min timestamp 2017-04-01 00:00:00+00:00 1086.633148 11.808373 6.990839 NaN NaN NaN NaN NaN 0.000435 0.020851 2017-04-02 00:00:00+00:00 1095.521806 14.183194 6.998986 0.008180 0.008147 0.000066 0.818000 0.814672 0.000434 0.020839 2017-04-03 00:00:00+00:00 1138.100972 21.106648 7.037116 0.038867 0.038130 0.001454 3.886656 3.813028 0.000704 0.026524 2017-04-04 00:00:00+00:00 1143.858542 14.609319 7.042163 0.005059 0.005046 0.000025 0.505893 0.504617 0.000458 0.021391 2017-04-05 00:00:00+00:00 1130.486875 15.630853 7.030404 -0.011690 -0.011759 0.000138 -1.168997 -1.175883 0.000380 0.019500

Some simple configurations for our charts.

import warnings warnings . simplefilter ( 'ignore' ) % matplotlib inline # Increase chart resolution % config InlineBackend . figure_format = 'retina' import seaborn as sns sns . set_style ( 'darkgrid' ) sns . mpl . rcParams [ 'figure.figsize' ] = ( 20.0 , 4.0 ) sns . mpl . rcParams [ 'savefig.dpi' ] = 150 sns . mpl . rcParams [ 'font.family' ] = 'serif' sns . mpl . rcParams [ 'font.size' ] = 12 sns . mpl . rcParams [ 'lines.linewidth' ] = 1.3

While we are at it, let us plot those three and take a look.

ax1 = d_df . loc [:, [ 'squared_log_return' , 'realized_variance_1min' , 'realized_variance_5min' ]] . plot () ax1 . xaxis . label . set_visible ( False )

BTCUSD estimated daily variance

As expected, the realized variance plots are visually smoother than squared log return plot, especially around January 2018.

Evaluation

There are three commonly accepted criteria to evalute a variety of models in financial area - Mean Absolute Error (MAE), Mean Absolute Percentage Error (MAPE) and Root Mean Square Error (RMSE).

Mean Absolute Error (MAE)

Mean Absolute Percentage Error (MAPE)

Root Mean Square Error (RMSE)

By definition, it is quite obvious that the smaller values calculated by those methods, the better performance of our model of interest. It is fairly easy to implement those criteria in Python as well.

from sklearn.metrics import mean_absolute_error , mean_squared_error from math import sqrt def mae ( observation , forecast ): error = mean_absolute_error ( observation , forecast ) print ( 'Mean Absolute Error (MAE): {:.3g}' . format ( error )) return error def mape ( observation , forecast ): observation , forecast = np . array ( observation ), np . array ( forecast ) # Might encounter division by zero error when observation is zero error = np . mean ( np . abs (( observation - forecast ) / observation )) * 100 print ( 'Mean Absolute Percentage Error (MAPE): {:.3g}' . format ( error )) return error def rmse ( observation , forecast ): error = sqrt ( mean_squared_error ( observation , forecast )) print ( 'Root Mean Square Error (RMSE): {:.3g}' . format ( error )) return error def evaluate ( pd_dataframe , observation , forecast ): """ :params :pd_dataframe pandas dataframe :observation column name of expected values :forecast column name of forecasted values :return the results of MAE, MAPE and RMSE, respectively """ first_valid_date = pd_dataframe [ forecast ] . first_valid_index () mae_error = mae ( pd_dataframe [ observation ] . loc [ first_valid_date :, ], pd_dataframe [ forecast ] . loc [ first_valid_date :, ]) mape_error = mape ( pd_dataframe [ observation ] . loc [ first_valid_date :, ], pd_dataframe [ forecast ] . loc [ first_valid_date :, ]) rmse_error = rmse ( pd_dataframe [ observation ] . loc [ first_valid_date :, ], pd_dataframe [ forecast ] . loc [ first_valid_date :, ]) ax = pd_dataframe . loc [:, [ observation , forecast ]] . plot () ax . xaxis . label . set_visible ( False ) return mae_error , mape_error , rmse_error

Exponentially Weighted Moving Average (EWMA)

Exponentially Weighted Moving Average (EWMA) is simply introduced here to help understanding ARCH. However, we not going for a deep dive other than deriving a few equations. You may skip this session.

One simple way of estimating future volatility is to use historical weighted volatilities. The assumption is that the further we go back in time, the less impact the previous volatility has on the future. Since this method simply chooses a constant rate over time, it is hard to miss that exponential decreasing pattern on the chart.

Before we can move forward, we need to simplify the equation from the volatility definition. When we have a larger number of observations and relatively smaller time window (i.e. daily data), it is safe to assume $n-1 \approx n$ and $(\bar{r} - r_i)^2 \approx r_i^2$. If we take it a bit further by squaring both sides, we get this.

The added subscript $t$ means today's volatility.

In other words, the variance of log return is the average of all squared log returns within the observation window $n$. However, this form simply allows the even contribution from all previous log returns. The next step is to implement exponentially weighted approach by introducing lambda $\lambda$.

The above equation can also be written like this.

Autoregressive Conditional Heteroskedasticity (ARCH)

In our previous posts, we have developed autoregressive models for univariate time series data that is stationary (AR), has a trend (ARIMA), and has a seasonal component (SARIMA). All of these models assume the data has a constant variance over time. However, financial markets and crypto trading data are exposed to consistent changes of the variance in time series - volatility.

Moderate changes in variance can be easily dealt with the power transformation, such as natural log and Box-Cox transformation. When the variance is largely time-dependent, the autoregressive process is introduced to treat the data.

In a nutshell, ARCH models the variance at a time step as a function of the residual errors from a mean process (e.g. a zero mean).

In other words, AR model is applied to the variance of a time series.

ARCH is only applicable where the time series is stationary, other than the change in variance, meaning it does not have a trend or seasonal component.

Similar to the equation mentioned in EWMA, $\alpha_0$, $\alpha_1$ … and $\alpha_q$ are considered as the parameters of the model instead of the weights.

$p$: The number of lag residual errors to include in the ARCH model.

Normally, researchers will start with ARIMA model to select the best combination of the parameters based on the lowest score of AIC, BIC or HQIC as we have discussed in a previous post.

The configuration for an ARCH model is best understood in the context of ACF and PACF plots of the variance of the time series. This can be achieved by subtracting the mean from each observation in the series and squaring the result, or just squaring the observation if you’re already working with white noise residuals from another model. - Jason Brownlee at Machine Learning Mastery

For the sake of this post, let us assume the mean log return is zero over time and run autocorrelations on those squared observations.

from statsmodels.graphics.tsaplots import plot_acf , plot_pacf _ = plot_acf ( d_df [ 'squared_log_return' ][ 1 :], lags = 100 , title = 'Squared Daily Log Return ACF' ) _ = plot_acf ( d_df [ 'realized_volatility_1min' ][ 1 :], lags = 100 , title = '1-min Daily Realized Variance ACF' ) _ = plot_acf ( d_df [ 'realized_volatility_5min' ][ 1 :], lags = 100 , title = '5-min Daily Realized Variance ACF' )

BTCUSD squared daily log return ACF

BTCUSD 1-min daily realized variance ACF

BTCUSD 5-min daily realized variance ACF

In case of squared log return, we have a positive correlation up to 62 lag time steps. On the other hand, the realized variance (both 1-min and 5-min) ACF correlograms have showed a positive correlation till 46 lag time steps. Since realized variance is generally considered as a more accurate estimate of the true variance, we are going to set the lag as 46.

Or is it? I have tried this lag. However, the evaluation errors are much larger than those using ARCH package default value. I could not figure out why. If you are having the same issue, please do let me know. And of course, if you know the solution, I am all ears!

Generalized Autoregressive Conditional Heteroskedasticity (GARCH)

Notice that we have only used squared log returns to estimate the variance so far. If we take one step further by incorporating previous variances into the model, we will get the following.

$p$: The number of lag variances to include in the GARCH model.

$q$: The number of lag residual errors to include in the GARCH model.

PAY ATTENTION TO NAME CONVENTIONS. $p$ in ARCH(p) specifies the number of lag residual errors, while $p$ in GARCH(p, q) is for the number of lag variances.

Exponential Generalized Autoregressive Conditional Heteroscedasticity (EGARCH)

Although there are thousands of academic papers and published books which focus on studying different volatility models in financial industry, we can find few studies to discuss the correlations between the existing forecasting models and cryptocurrencies markets. Interestingly, EGARCH seems to perform better than GARCH and its other variants over all, according to Naimy and Zheng.

Exponential Generalized Autoregressive Conditional Heteroscedasticity (EGARCH) is a variant of GARCH. It integrated the asymmetric influence on volatility from both positive and negative impacts by introducing natural log on the variable of interest. There was an empirical belief that the negative factors can perform stronger volatility than positive factors, even with the same absolute value. The practice of asymmetry volatility was further confirmed by the famous Leverage Effect Theory proposed by Black in 1976 “Studies in Stock Price Volatility Change”. Negative values are allowed in EGARCH, while conditional variance is positive.

In ARCH package, EGARCH variance dynamics are

where $e_{t}=\epsilon_{t}/\sigma_{t}$

ARCH Python Package

from arch import arch_model

ARCH is a package written in Python (with Cython and/or Numba used to improve performance). As discussed in ARCH official document, there are four types of forecasting schemes (basic, fixed window, rolling window and recursive forecasting) at our disposal. I am going to demonstrate how to implement 3 out of those 4 schemes via ARCH, GARCH and EGARCH.

model_in_sample = arch_model ( d_df [ 'log_return' ][ 1 :], mean = 'Zero' , vol = 'ARCH' , p = 1 , dist = 't' ) res_in_sample = model_in_sample . fit ( disp = 'off' )

ARCH covergence warning

You might encounter an error message like the one above - “ConvergenceWarning”. Recall we scale up 100x of the log return? This is because ARCH package cannot handle small numbers by defaut. We either scale up our data or re-adjust the default value.

# Scale up 100x model_in_sample = arch_model ( d_df [ 'log_return_100x' ][ 1 :], mean = 'Zero' , vol = 'ARCH' , p = 1 , dist = 't' ) res_in_sample = model_in_sample . fit ( disp = 'off' )

Alternatively, we can pass in some options in the fit() function.

model_in_sample = arch_model ( d_df [ 'log_return' ][ 1 :], mean = 'Zero' , vol = 'ARCH' , p = 1 , dist = 't' ) # Pass in options # 'ftol' is 1e-06 by defaut # Re-adjust to 1e-01 for our data res_in_sample = model_in_sample . fit ( disp = 'off' , options = { 'ftol' : 1e-01 })

Let us just stay with passing in ftol value for the rest of this tutorial. Feel free to experiment with the scaling up option.

Fixed Window Forecasting

We can use a fixed window forecasting scheme in ARCH package and see what happens.

Fixed-windows forecasting uses data up to a specified date to generate all forecasts after that date.

The nice part of this scheme is that we can pass the entire data set when initializing the model. We just need to provide the last observation timestamp when calling fit() function. The following forecast() function will, by default, produce forecasts after that timestamp.

We can create a simple function to run ARCH, GARCH and EGARCH all together.

def fixed_window_forecast ( pd_dataframe ): models = [ 'ARCH' , 'GARCH' , 'EGARCH' ] for model in models : mod = arch_model ( pd_dataframe [ 'log_return' ][ 1 :], mean = 'Zero' , vol = model , dist = 't' ) res = mod . fit ( disp = 'off' , last_obs = 356 , options = { 'ftol' : 1e-01 }) forecasts = res . forecast () pd_dataframe [ 'fixed_window_{}' . format ( model )] = forecasts . variance [ 'h.1' ] print ( '------------------------------------------------' ) print ( 'Fixed Window Forcasting {}' . format ( model )) evaluate ( pd_dataframe , 'realized_variance_1min' , 'fixed_window_{}' . format ( model )) fixed_window_forecast ( d_df )

------------------------------------------------ Fixed Window Forcasting ARCH Mean Absolute Error (MAE): 0.00123 Mean Absolute Percentage Error (MAPE): 3.11e+06 Root Mean Square Error (RMSE): 0.00148 ------------------------------------------------ Fixed Window Forcasting GARCH Mean Absolute Error (MAE): 0.000637 Mean Absolute Percentage Error (MAPE): 8.48e+05 Root Mean Square Error (RMSE): 0.00122 ------------------------------------------------ Fixed Window Forcasting EGARCH Mean Absolute Error (MAE): 0.000582 Mean Absolute Percentage Error (MAPE): 3.88e+05 Root Mean Square Error (RMSE): 0.00125

Fixed window forecasting ARCH

Fixed window forecasting GARCH

Fixed window forecasting EGARCH

Rolling Window Forecasting

In addition, we can also adopt a rolling window scheme and check if there is any improvement.

Rolling window forecasts use a fixed sample length and then produce one-step from the final observation.

Basically, the model generates one-step ahead forecast and drops the oldest daily log return while including the newest daily return for a re-fit. Similarly, we will have a simple function to run ARCH, GARCH and EGARCH all together.

def rolling_window_forecast ( pd_dataframe ): # Set rolling window window = 365 models = [ 'ARCH' , 'GARCH' , 'EGARCH' ] for model in models : index = d_df [ 1 :] . index end_loc = np . where ( index >= d_df . index [ window ])[ 0 ] . min () forecasts = {} for i in range ( len ( d_df [ 1 :]) - window + 2 ): mod = arch_model ( pd_dataframe [ 'log_return' ][ 1 :], mean = 'Zero' , vol = model , dist = 't' ) res = mod . fit ( first_obs = i , last_obs = i + end_loc , disp = 'off' , options = { 'ftol' : 1e03 }) temp = res . forecast () . variance fcast = temp . iloc [ i + end_loc - 1 ] forecasts [ fcast . name ] = fcast forecasts = pd . DataFrame ( forecasts ) . T pd_dataframe [ 'rolling_window_{}' . format ( model )] = forecasts [ 'h.1' ] print ( '------------------------------------------------' ) print ( 'Rolling Window Forcasting {}' . format ( model )) evaluate ( pd_dataframe , 'realized_variance_1min' , 'rolling_window_{}' . format ( model )) rolling_window_forecast ( d_df )

------------------------------------------------ Rolling Window Forcasting ARCH Mean Absolute Error (MAE): 0.000858 Mean Absolute Percentage Error (MAPE): 1.52e+06 Root Mean Square Error (RMSE): 0.00121 ------------------------------------------------ Rolling Window Forcasting GARCH Mean Absolute Error (MAE): 0.000616 Mean Absolute Percentage Error (MAPE): 6.41e+05 Root Mean Square Error (RMSE): 0.00117 ------------------------------------------------ Rolling Window Forcasting EGARCH Mean Absolute Error (MAE): 0.00058 Mean Absolute Percentage Error (MAPE): 3.5e+05 Root Mean Square Error (RMSE): 0.00126

Rolling window forecasting ARCH

Rolling window forecasting GARCH

Rolling window forecasting EGARCH

Recursive Forecast Generation

There is one more forecasting scheme provided by ARCH package.

Recursive is similar to rolling except that the initial observation doesn’t change. This can be easily implemented by dropping the first_obs input.

def recursive_forecast ( pd_dataframe ): # Set rolling window window = 365 models = [ 'ARCH' , 'GARCH' , 'EGARCH' ] for model in models : index = d_df [ 1 :] . index end_loc = np . where ( index >= d_df . index [ window ])[ 0 ] . min () forecasts = {} for i in range ( len ( d_df [ 1 :]) - window + 2 ): mod = arch_model ( pd_dataframe [ 'log_return' ][ 1 :], mean = 'Zero' , vol = model , dist = 't' ) res = mod . fit ( last_obs = i + end_loc , disp = 'off' , options = { 'ftol' : 1e03 }) temp = res . forecast () . variance fcast = temp . iloc [ i + end_loc - 1 ] forecasts [ fcast . name ] = fcast forecasts = pd . DataFrame ( forecasts ) . T pd_dataframe [ 'recursive_{}' . format ( model )] = forecasts [ 'h.1' ] print ( '------------------------------------------------' ) print ( 'Recursive Forcasting {}' . format ( model )) evaluate ( pd_dataframe , 'realized_variance_1min' , 'recursive_{}' . format ( model )) recursive_forecast ( d_df )

------------------------------------------------ Recursive Forcasting ARCH Mean Absolute Error (MAE): 0.000867 Mean Absolute Percentage Error (MAPE): 1.69e+06 Root Mean Square Error (RMSE): 0.00121 ------------------------------------------------ Recursive Forcasting GARCH Mean Absolute Error (MAE): 0.000583 Mean Absolute Percentage Error (MAPE): 6.49e+05 Root Mean Square Error (RMSE): 0.00117 ------------------------------------------------ Recursive Forcasting EGARCH Mean Absolute Error (MAE): 0.000546 Mean Absolute Percentage Error (MAPE): 3.15e+05 Root Mean Square Error (RMSE): 0.00116

Recursive forecasting ARCH

Recursive forecasting GARCH

Recursive forecasting EGARCH

EGARCH has outperformed the other two in all 3 forecasting schemes based on MAE, MAPE, and RMSE.

The End

I hope you have enjoyed this piece. My apologies for this long post. There are a couple of other directions we could have pursued further. For instance, simulation and bootstrap are two methods which could be implemented in EGARCH for out-of-sample forecasting. It would be very interesting to see the simulation paths from their results. I guess those exciting challenges will be left to our readers who are eager to investigate further.

Again, thank you for your time and patience. All source code will be uploaded to my github. Let me know if you have any questions/comments/proposals. My contact information can be found at the beginning of this post.

Stay calm and happy trading!