Comparison of Volatility Forecasts from Machine Learning Models and Conditional Heteroskedasticity Models

1. Abstract

The paper compares the accuracy of volatility forecasts of GARCH(1,1), EGARCH(1,1), GJR-GARCH(1,1) models and XGBoost, LightGBM machine learning models. Machine learning models are taken from the dquant library, and conditional heteroskedasticity from the arch library. Data is taken from 01.01.2014 to 31.12.2024 from the yfinance library. Tickers: BTC-USD (bitcoin/dollar), ETH-USD (ether/dollar), EURUSD=X (euro/dollar), BZ=F (oil), GC=F (gold), SI=F (silver), SPY. Training was carried out on data from 01.01.2014 to 01.01.2023, validation from 01.01.2023 to 31.12.2024. Quality metrics of forecasts: QLIKE, MAE. The Diebold-Mariano test is used to check the statistical significance of the results. Machine learning models showed a statistically significant advantage on all assets in terms of MAE and on 4 out of 7 assets by QLIKE.

2. Introduction

Volatility is a statistical financial indicator that characterizes the variability of the price of a financial instrument. Volatility in finance is usually defined as the standard deviation or variance of the return of an asset over a certain period of time.

Quantifying the degree of price variability is important in trading; volatility is used everywhere, from risk management to trading strategies based on volatility forecasting.

A revolutionary breakthrough in modeling volatility was made by Robert Engle in 1982, when he proposed the ARCH (Autoregressive Conditional Heteroskedasticity) model. In the ARCH model, the conditional variance at the current time is represented as a linear function of the squares of past values of the process.

Mathematically, the ARCH(q) model can be represented as follows:

where:

• σ_t² — conditional variance at time t;
• ω — constant (usually positive);
• α_i — model parameters, which must be positive to ensure that the variance is positive;
• ε_t−i — innovations of the process (usually assumed to be independent and identically distributed with zero mean).

For the development of the ARCH model, Robert Engle was awarded the Nobel Prize in Economics in 2003, which highlights the significance of this contribution to financial econometrics.

In 1986, Tim Bollerslev proposed a generalization of the ARCH model, which became known as GARCH (Generalized Autoregressive Conditional Heteroskedasticity). The key idea was to include past conditional variances in the equation for the current conditional variance:

where β_j — parameters reflecting the influence of past conditional variances.

This formulation turned out to be more parsimonious in terms of parameters, and often the GARCH(1,1) model was sufficient to adequately describe volatility, whereas an equivalent ARCH model would require an infinite number of parameters.

Machine learning models allow you to process a huge amount of data; they are used in trading, including in predicting volatility. For this, in most cases, boosting machine learning models are used: they do not require as much data as neural networks, but at the same time their forecasts do not lag behind the forecasts of neural network models.

The purpose of this work is to compare the accuracy of 1-step-ahead volatility forecasts on 7 different assets using GARCH family models and gradient boosting models. GARCH models will be taken from the arch library, gradient boosting models from the author's library dquant. The accuracy of forecasts will be compared on the MAE and QLIKE metrics. The statistical significance of the most accurate models will be checked on the Diebold-Mariano test.

3. Methodology

3.1 Data

The data is taken from the python library yfinance. For training, a period of 9 years of OHLC data from the daily time frame was taken. For validation, 2 years. The following assets are used: BTC-USD (bitcoin/dollar), ETH-USD (ether/dollar), EURUSD=X (euro/dollar), BZ=F (oil), GC=F (gold), SI=F (silver), SPY.

It is important to use different markets with different price movements. For example, if you already have a euro/dollar pair, you do not need to use a pound/dollar pair. Due to the strong correlation between these two instruments, the results will be similar.

After obtaining the data, it is necessary to divide it into training and validation data. Our models will be trained on the training data. The validation data will only be used after training. The python code is provided:

START_DATE = '2014-01-01'
SPLIT_DATE = '2023-01-01'
END_DATE   = '2024-12-31'

def get_data(ticker, start_date, end_date):
    print(f"Loading {ticker} from {start_date} to {end_date}...")
    raw = yf.download(ticker, start=start_date, end=end_date, auto_adjust=True)
    df = pd.DataFrame({
        'open':   raw[('Open', ticker)].values,
        'high':   raw[('High', ticker)].values,
        'low':    raw[('Low', ticker)].values,
        'close':  raw[('Close', ticker)].values,
        'volume': raw[('Volume', ticker)].values
    }, index=raw.index)
    return df

df = get_data('BTC-USD', START_DATE, END_DATE)

train_mask = df.index < SPLIT_DATE
test_mask = df.index >= SPLIT_DATE

df_train = df[train_mask].copy()
df_test = df[test_mask].copy()

3.2 Models used, their training and validation

3.2.1 Machine learning models

XGBoost and LightGBM machine learning models from the DQuant python library will be used. This library provides a convenient interface for training ML models for volatility prediction, and it was created specifically for this purpose.

We will use returns as input data and Parkinson volatility as targets. The Parkinson volatility formula is as follows:

The library provides a convenient interface for creating features and targets, but for clarity, we implement these functions manually. The library allows you to embed your own functions for creating features and targets.

def parkinson_func(df):
    df = df.copy()
    df = df.iloc[1:]
    return np.array(np.sqrt(
        (1 / (4 * np.log(2))) * (np.log(df['high'] / df['low'])) ** 2
    ))

def return_func(df):
    return np.array(df['close'].pct_change().dropna())

We embed these functions into the model training code:

from dquant.models import VolClustXGB, VolClustLightGBM

def dquant_xgb_train():
    model_dq = VolClustXGB({}, early_stopping=True, output=False, loss='QLIKE')
    model_dq.fit(
        df_train,
        feature_list=FEATURES,
        input_bars=INPUT_BARS,
        horizon=HORIZON,
        trees_count=TREES_COUNT,
        show_results=False,
        feature_func=return_func,
        target_func=parkinson_func
    )
    return model_dq

The training will be based on the QLIKE loss function, which was created specifically for training to predict volatility, and MAE.

Training was conducted with the default hyperparameters. For XGBoost, these are:

{
    'learning_rate': 0.1,
    'max_depth': 6,
    'min_child_weight': 5,
    'gamma': 0.1,
    'subsample': 0.8,
    'colsample_bytree': 0.8,
    'reg_alpha': 0.0,
    'reg_lambda': 1.0,
    'random_state': 42,
    'tree_method': 'hist',
    'device': 'cpu'
}

For LightGBM:

{
    'learning_rate': 0.1,
    'max_depth': 6,
    'num_leaves': 2 ** 6 - 1,
    'min_child_samples': 5,
    'min_split_gain': 0.1,
    'bagging_fraction': 0.8,
    'bagging_freq': 1,
    'feature_fraction': 0.8,
    'reg_alpha': 0.0,
    'reg_lambda': 1.0,
    'random_state': 42,
    'verbosity': -1,
    'boosting_type': 'gbdt'
}

The optimal number of trees is selected during the training process.

Code for obtaining model predictions:

def dquant_xgb_forecasting():
    model_dq = dquant_xgb_train()
    dquant_preds = []
    dquant_dates = []

    for i in range(INPUT_BARS, len(df_test)):
        test_date = df_test.index[i]
        available_data = df.loc[:test_date].iloc[-INPUT_BARS:].copy()

        if len(available_data) < INPUT_BARS:
            continue

        try:
            forecast_vals = model_dq.forecast(available_data, show=False)
            pred_vol = forecast_vals[-1] if HORIZON > 1 else forecast_vals[0]
        except Exception as e:
            print(f"DQuant forecast error on {test_date.date()}: {e}")
            continue

        actual_vol = df_test.loc[test_date, 'parkinson_vol']
        dquant_dates.append(test_date)
        dquant_preds.append({
            'date': test_date,
            'actual': actual_vol,
            'predicted': pred_vol
        })

    df_dquant = pd.DataFrame(dquant_preds).set_index('date')
    return df_dquant

The same should be done with VolClustLightGBM.

3.2.2 Conditional Heteroskedasticity Models

Conditional heteroskedasticity models will be used: GARCH(1,1), EGARCH(1,1), GJR-GARCH(1,1) from the python library arch.

GARCH (Generalized Autoregressive Conditional Heteroskedasticity) — conditional heteroskedasticity model. The formula looks like this:

Where:

• σ_t² — conditional variance at time t;
• ε_t−1 — standardized residual at time t−1;
• ω — constant;
• β — persistence coefficient;
• α — impact of the shock size.

Implementation in python:

def garch_forecasting():
    garch_preds = []
    returns_full = df['returns'].dropna() * 100  # scaling for GARCH

    for test_date in df_test.index:
        hist_returns = returns_full[:test_date].iloc[:-1]

        if len(hist_returns) < 500:
            continue

        try:
            model_garch = arch_model(hist_returns, vol='GARCH', p=1, q=1, dist='normal')
            fitted = model_garch.fit(disp='off')
            forec = fitted.forecast(horizon=HORIZON, reindex=False)
            cond_var = forec.variance.values[-1, 0]
            pred_vol = np.sqrt(cond_var) / 100.0  # reverse scaling
        except Exception as e:
            print(f"GARCH forecast error on {test_date.date()}: {e}")
            continue

        actual_vol = df.loc[test_date, 'parkinson_vol']
        garch_preds.append({
            'date': test_date,
            'actual': actual_vol,
            'predicted': pred_vol
        })

    df_garch = pd.DataFrame(garch_preds).set_index('date')
    return df_garch

EGARCH (Exponential GARCH) is a modification of GARCH that takes into account the asymmetric effects of volatility.

Where:

• z_t−1 = ε_t−1 / σ_t−1 — standardized residual;
• ω — constant;
• β — persistence coefficient;
• α — impact of shock magnitude;
• γ — asymmetry parameter.

The code will be the same, except for the code of the model itself:

model_garch = arch_model(hist_returns, vol='EGARCH', p=1, q=1, dist='normal')

GJR-GARCH is another model for asymmetric effects, named after its authors (Glosten, Jagannathan, and Runkle).

Where:

• σ_t² — conditional variance at time t;
• ε_t−1 — standardized residual at time t−1;
• ω — constant;
• β — persistence coefficient;
• α — impact of shock magnitude;
• γ — asymmetry parameter (leverage effect), which enhances the impact of negative shocks;
• I_t−1 — indicator function, which is equal to 1 if ε_t−1 is negative, otherwise 0.

The code will be the same, except for the code of the model itself:

model_garch = arch_model(hist_returns, vol='GARCH', p=1, q=1, o=1, dist='normal')

o=1 adds the GJR-GARCH effect. The library handles all the mathematical implementation.

3.3 Metrics

The Mean Absolute Error (MAE) and Quasi-Likelihood (QLIKE) quality metrics were used.

Mean Absolute Error (MAE) is a regression model quality metric calculated as the arithmetic mean of the absolute values of the errors (the differences between the predicted and actual values).

Formula:

Where:

• y — true value;
• ŷ — predicted value.

It will be taken from the python library scikit-learn: from sklearn.metrics import mean_absolute_error.

Quasi-Likelihood (QLIKE) is a metric created specifically for evaluating volatility forecasts.

Formula:

Where:

• y — true value;
• ŷ — predicted value.

There is no such function in scikit-learn, so it will be implemented manually:

def qlike(y_true, y_pred):
    sigma2_true = y_true ** 2
    sigma2_pred = np.maximum(y_pred ** 2, 1e-10)
    return np.mean(np.log(sigma2_pred) + sigma2_true / sigma2_pred)

3.4 Test

After identifying the most accurate machine learning models and conditional heteroskedasticity models, we will compare the results of their predictions on the Diebold-Mariano test, to identify the statistical superiority in the predictions of any of the models.

The Diebold-Mariano test is used to determine how much the 2 predictions differ.

First, we calculate a series of errors. For MAE:

For QLIKE:

Where:

• e — error series;
• y — true values series;
• f — predictions series.

Next, we calculate the difference in loss between the two series:

If the difference in losses d < 0, then the forecast of the e₁ series is more accurate than the forecast of e₂.

Since we have a 1-step-ahead forecast, it is not necessary to consider autocorrelation. The test formula without considering autocorrelation is:

We calculate the p-value:

Where Φ is the distribution function of the standard normal law:

The code:

def dm_test(d):
    n = len(d)
    mean_d = np.mean(d)
    var_d = np.var(d)
    dm_stat = np.sqrt(n) * mean_d / np.sqrt(var_d)
    p_value = 2 * (1 - stats.norm.cdf(np.abs(dm_stat)))
    return dm_stat, p_value

err_dq = abs(the_best_dquant_aligned['actual'] - the_best_dquant_aligned['predicted'])
err_garch = abs(the_best_garch_aligned['actual'] - the_best_garch_aligned['predicted'])
d = err_dq - err_garch
dm_stat, p_value = dm_test(d)

3.5 The Code

All the source code is available on GitHub: github.com/artrdon/DQuant_Research.

4. Results

4.1 Results when using the QLIKE loss function

Diebold-Mariano test results for the models with the best accuracy on each asset:

4.2 Results when using the MAE loss function

Diebold-Mariano test results for the models with the best accuracy on each asset:

5. Conclusion

In this study, we compared conditional heteroskedasticity models from the arch library and machine learning models from the dquant library on 7 different financial instruments.

On the Diebold-Mariano test, when using MAE as the loss function, gradient boosting models showed a statistically significant advantage on all assets. This shows that machine learning models, particularly gradient boosting models, can capture nonlinear dependencies that GARCH models do not detect.

When using QLIKE as the loss function when training ML models, a statistical advantage was found for the MAE metric on all assets. For the QLIKE metric, a statistical advantage was found on 4 out of 7 assets. No statistical advantage was found on EURUSD=X, BZ=F, SI=F.

When using MAE as the loss function when training ML models, a statistical advantage was found for the MAE metric on all assets. For the QLIKE metric, a statistical advantage was found on 4 out of 7 assets. No statistical advantage was found for BTC-USD, EURUSD=X, GC=F. However, a statistical advantage was found for silver in favor of GARCH models.

The comparison was conducted on daily data, and the results may differ if intraday data is used.

Gradient boosting hyperparameters were not optimized, and the accuracy of machine learning models may be higher if optimal hyperparameters are used.

Only returns were used as a feature for machine learning models, and the accuracy of predictions may be higher if additional features are used. There was no walk forward validation; instead, the data was statically split into training and validation sets.

6. References

1. github.com/artrdon/dquant
2. mlgu.ru/1778
3. real-statistics.com — Diebold-Mariano test
4. elma365.com — MSE
5. Patton — Volatility forecast comparison (JoE, 2011)
6. mlgu.ru/6807
7. quantinsti.com — GARCH & GJR-GARCH volatility forecasting in Python
8. Wikipedia — Волатильность
9. Bollerslev — Generalized ARCH (JoE, 1986)