harlstmgarch

A Python implementation of the three-stage hybrid HAR-LSTM-GARCH framework for realized-volatility forecasting, with Bayesian-optimized LSTM tuning, a full suite of econometric/ML benchmarks, and the Model Confidence Set test for statistically rigorous model comparison.

Author: Belhimer Hocine — Maître de conférences A (Lecturer A), Higher School of Accounting and Finance of Constantine (ESCF Constantine) — hbelhimer@escf.constantine.dz

Method reference: Ben Romdhane, W., & Boubaker, H. (2026). A Hybrid HAR-LSTM-GARCH Model for Forecasting Volatility in Energy Markets. Journal of Risk and Financial Management, 19(2), 77. https://doi.org/10.3390/jrfm19010077

---

Table of contents

• The model
• Installation
• Quick start
• API reference (detailed syntax)
  • harlstmgarch.data
  • harlstmgarch.har
  • harlstmgarch.lstm
  • harlstmgarch.tune
  • harlstmgarch.garch
  • harlstmgarch.hybrid
  • harlstmgarch.benchmarks
  • harlstmgarch.metrics
  • harlstmgarch.plots
• End-to-end example
• Reproducing the paper
• Methodology ↔ paper equations
• Limitations
• Citation
• Changelog
• License

---

The model

The framework follows a sequential decomposition principle — each stage models what the previous stage leaves behind:

Stage	Component	Role	Eqs.
1	HAR (Corsi, 2009)	Linear filter for persistent, multi-scale dynamics (daily/weekly/monthly RV, OLS)	6, 24–25
2	LSTM	Recurrent net learning the non-linear patterns left in the HAR residuals; tuned by Bayesian Optimization	17–23, 27–28
3	GARCH(1,1)	Conditional variance of the hybrid forecast errors → time-varying prediction intervals	30–33

            RV_t, RV_t^w, RV_t^m
                     │
        ┌────────────▼────────────┐   e_t = RV_t − L_t
        │  Stage 1: HAR (OLS)      ├──────────────┐
        └────────────┬────────────┘              │
              L_{t+1}│                            ▼
                     │              ┌──────────────────────────┐
                     │              │ Stage 2: LSTM (Bayes-opt) │
                     │              │   ê_{t+1}=f(e_t,…)        │
                     │              └────────────┬─────────────┘
                     ▼                           │ ê_{t+1}
   H_{t+1}=L_{t+1}+ê_{t+1}  ◄────────────────────┘
                     │   z_t = RV_t − H_t
        ┌────────────▼────────────┐
        │  Stage 3: GARCH(1,1)     │  σ̂²_{z,t+1}
        └────────────┬────────────┘
                     ▼
   Point: H_{t+1}   Interval: H_{t+1} ± 1.96·σ̂_{z,t+1}

---

Installation

pip install harlstmgarch

With development tools (build / twine / pytest):

pip install "harlstmgarch[dev]"

Requirements (installed automatically): Python ≥ 3.10, numpy, pandas, scipy, statsmodels, arch, tensorflow ≥ 2.13, scikit-learn, matplotlib.

---

Quick start

import pandas as pd
from harlstmgarch import (
    realized_vol_from_returns, har_components, chronological_split,
    HARLSTMGARCH, metrics,
)

# 1. Realized volatility from a daily price series (paper Eq. 2)
prices = pd.read_csv("BrentCrude.csv", parse_dates=["DATE"]).set_index("DATE")["Close"]
rv = realized_vol_from_returns(prices, window=22)["RV"]

# 2. HAR features (RV in levels) + chronological 70/15/15 split
df = har_components(rv)
train, val, test = chronological_split(df, train=0.70, val=0.15)

# 3. Tune the residual LSTM by Bayesian Optimization, then fit all 3 stages
model = HARLSTMGARCH(lstm_kwargs={"lookback": 20})
model.tune_lstm(train, val, n_calls=20)      # GP + Matérn 5/2
model.fit(train, val)

# 4. Strictly one-step-ahead out-of-sample forecast (point + 95% interval)
fc = model.forecast(test)

print(metrics.scoreboard(test["RV"],
      {"HAR": fc.har, "HAR-LSTM-GARCH": fc.point}))
print("95% coverage:",
      metrics.interval_coverage(test["RV"], fc.lower, fc.upper))

---

API reference (detailed syntax)

All public names are importable directly from the top level, e.g. from harlstmgarch import HARLSTMGARCH, metrics.

harlstmgarch.data

load_realized_measures(path, date_col="DATE", rv_col="RV") -> pd.DataFrame

Load a realized-measures CSV (dates as YYYYMMDD). The input rv_col is treated as a realized variance; its square root is returned as a RV column (volatility scale), indexed by date.

realized_vol_from_returns(prices, window=22, annualize=True) -> pd.DataFrame

Build daily RV from a daily closing-price pd.Series (paper Eq. 2): RV_t = sqrt( (252/window) · Σ_{i=0}^{window-1} r_{t-i}² ) with r_t = ln(P_t/P_{t-1}). Returns columns RV and returns.

har_components(rv, weekly=5, monthly=22) -> pd.DataFrame

HAR design matrix. For each day t: target RV (=RV_t) and the lagged regressors RV_d = RV_{t-1}, RV_w = mean(last 5), RV_m = mean(last 22) — all computed from information available at t-1 (no look-ahead).

chronological_split(df, train=0.70, val=0.15) -> (train, val, test)

Time-ordered split (never shuffled). Test fraction is 1 − train − val.

harlstmgarch.har

class HARModel
    .fit(df) -> self                 # OLS on RV ~ RV_d + RV_w + RV_m
    .predict(df) -> pd.Series        # linear forecast
    .residuals(df) -> pd.Series      # RV − prediction
    .coefficients() -> pd.DataFrame  # Estimate / Std.Error / t / p  (cf. Table 6)
    .rsquared -> float

harlstmgarch.lstm

class ResidualLSTM(lookback=20, units=32, learning_rate=5.5e-4,
                   batch_size=32, epochs=300, patience=25, seed=42)
    .fit(train_resid, val_resid=None) -> self
    .predict(resid_history, resid_future_index) -> pd.Series

A single-layer LSTM (+ linear dense head) that forecasts the next residual from its own history; inputs are min/max-scaled to [-1, 1], loss is MSE (Adam), early stopping on validation loss. predict is rolling and strictly one-step-ahead over resid_future_index.

harlstmgarch.tune

class BayesianLSTMTuner(lookback=20, n_calls=20, n_init=6, patience=20,
                        max_eval_epochs=120, seed=42)
    .optimize(e_train, e_val) -> self
    .best_params_     # dict: {units, learning_rate, batch_size, epochs}
    .best_score_      # best validation MSE
    .history_         # pd.DataFrame of every trial

SEARCH_SPACE   # dict of the paper's ranges (units, learning_rate, batch_size, epochs)

Gaussian-Process Bayesian Optimization (Matérn 5/2 kernel + Expected Improvement) of the residual-LSTM hyperparameters, minimizing one-step validation MSE — reproduces paper §3.4.2 using only scikit-learn. max_eval_epochs caps epochs during search for speed (set None to honour the sampled value exactly).

harlstmgarch.garch

class ResidualGARCH(rescale_factor=100.0)
    .fit(z) -> self                       # zero-mean GARCH(1,1), Normal innovations
    .conditional_sigma(z) -> pd.Series    # recursive one-step σ_t (OOS)
    .parameters() -> pd.DataFrame         # ω, α, β with SE / t / p
    .persistence -> float                 # α + β

Stage 3: models the conditional variance of the hybrid forecast errors z_t (σ_t² = ω + α·z_{t-1}² + β·σ_{t-1}²). The conditional σ gives the prediction intervals.

harlstmgarch.hybrid

class HARLSTMGARCH(lstm_kwargs={}, z_level=1.96)
    .tune_lstm(train, val, n_calls=20, n_init=6,
               max_eval_epochs=120, seed=42) -> dict   # runs Bayesian Opt., rebuilds the LSTM
    .fit(train, val=None) -> self                      # fits all three stages
    .forecast(test) -> HybridForecast
    .har, .lstm, .garch, .tuner_                        # fitted components

@dataclass HybridForecast
    .point   # H_t = HAR + LSTM-residual forecast
    .har     # stage-1 linear forecast
    .sigma   # stage-3 conditional std of the forecast error
    .lower   # point − z_level·sigma
    .upper   # point + z_level·sigma

train/val/test are frames produced by har_components. Set z_level for other interval widths (1.96 → 95%, 2.576 → 99%).

harlstmgarch.benchmarks

class GARCHBenchmark(returns, rescale_factor=1.0)      # GARCH(1,1), Skewed-Student
class GJRGARCHBenchmark(returns, rescale_factor=1.0)   # GJR-GARCH(1,1), leverage
    .fit() -> self
    .forecast_sigma(index) -> pd.Series                # one-step conditional volatility

class StandaloneLSTM(lookback=20, units=97, learning_rate=5.81e-3,
                     batch_size=12, epochs=198, patience=25, seed=42)
    .fit(train_rv, val_rv=None) -> self
    .forecast(rv_history, index) -> pd.Series          # LSTM on raw RV (no HAR)

The competitors used in the paper's Table 4 comparison. StandaloneLSTM defaults to the paper's standalone hyperparameters (97 units, lr 5.81e-3, …).

harlstmgarch.metrics

# point-forecast losses
rmse(actual, forecast) -> float
mae(actual, forecast) -> float
mape(actual, forecast) -> float
r2(actual, forecast) -> float
qlike(actual, forecast) -> float            # mean( RV/RV̂ − log(RV/RV̂) − 1 )
qlike_loss_series(actual, forecast) -> np.ndarray   # per-obs QLIKE (for the MCS)

scoreboard(actual, forecasts: dict[str, Series]) -> pd.DataFrame
    # one row per model: RMSE, MAE, MAPE, R2, QLIKE

# statistical tests
mcleod_li(residuals, lags=20) -> (Q, p_value)        # ARCH-type non-linearity (Eq. 26)
diebold_mariano(actual, f1, f2, power=2) -> (DM, p)  # equal predictive accuracy
interval_coverage(actual, lower, upper) -> float     # empirical PI coverage
model_confidence_set(losses: dict[str, ndarray], alpha=0.05,
                     n_boot=5000, block=None, seed=42) -> pd.DataFrame
    # Hansen-Lunde-Nason (2011) MCS; columns: MCS_p, eliminated_order
    # models with MCS_p >= alpha are in the confidence set

harlstmgarch.plots

Publication-quality matplotlib figures (each returns a Figure; pass path=... to save a PNG):

plots.timeseries_with_split(rv, splits, title, path=None)
plots.forecast_comparison(actual, forecasts, title, path=None)
plots.prediction_band(actual, point, lower, upper, title, path=None)
plots.scatter_fit(actual, forecasts, title, path=None)
plots.residual_diagnostics(har_resid, hybrid_resid, title, lags=30, path=None)
plots.metric_bars(scoreboard, metric, title, path=None, lower_is_better=True)
plots.training_history(history, title, path=None)

---

End-to-end example

import pandas as pd
from harlstmgarch import (
    realized_vol_from_returns, har_components, chronological_split,
    HARModel, HARLSTMGARCH, GARCHBenchmark, GJRGARCHBenchmark,
    StandaloneLSTM, metrics,
)

# --- data ---------------------------------------------------------------
raw = pd.read_csv("BrentCrude.csv", parse_dates=["DATE"]).set_index("DATE")
out = realized_vol_from_returns(raw["Close"], window=22)
rv, ret = out["RV"], out["returns"]
df = har_components(rv)
train, val, test = chronological_split(df, 0.70, 0.15)

# --- stage 1: HAR + justification for the LSTM --------------------------
har = HARModel().fit(train)
print(har.coefficients().round(5))                 # cf. paper Table 6
Q, p = metrics.mcleod_li(har.residuals(train), lags=20)
print(f"McLeod-Li Q(20) = {Q:.1f}, p = {p:.3g}")   # non-linearity remains

# --- stages 1-3: Bayesian-optimized hybrid ------------------------------
model = HARLSTMGARCH(lstm_kwargs={"lookback": 20})
print("best LSTM:", model.tune_lstm(train, val, n_calls=20))
model.fit(train, val)
fc = model.forecast(test)
print(model.garch.parameters().round(4), "persistence:", model.garch.persistence)

# --- benchmarks + scoreboard --------------------------------------------
rv_lstm = StandaloneLSTM().fit(train["RV"], val["RV"])
hist = pd.concat([train["RV"], val["RV"], test["RV"]])
forecasts = {
    "HAR": fc.har,
    "LSTM": rv_lstm.forecast(hist, test.index),
    "GARCH": GARCHBenchmark(ret).fit().forecast_sigma(test.index),
    "GJR-GARCH": GJRGARCHBenchmark(ret).fit().forecast_sigma(test.index),
    "HAR-LSTM-GARCH": fc.point,
}
board = metrics.scoreboard(test["RV"], forecasts)
print(board.round(5))                              # cf. paper Table 4

# --- Model Confidence Set on QLIKE --------------------------------------
losses = {k: metrics.qlike_loss_series(test["RV"], f.reindex(test.index))
          for k, f in forecasts.items()}
print(metrics.model_confidence_set(losses).round(4))   # cf. paper Table 5

---

Reproducing the paper

examples/run_case_study.py runs the entire pipeline (data → HAR → McLeod-Li → Bayesian-optimized LSTM → GARCH → benchmarks → scoreboard → MCS → figures) and writes all tables/figures to assets/:

python examples/run_case_study.py --prices BrentCrude.csv --n-calls 20

If no price file is supplied it falls back to a bundled realized-measures CSV so the pipeline is runnable out of the box; supply the Brent series to reproduce the published results.

---

Methodology ↔ paper equations

Concept	Paper	Code
Realized volatility (22-day)	Eq. 2	realized_vol_from_returns
HAR-RV model	Eq. 6, Table 6	HARModel
LSTM gates / cell	Eqs. 17–23	ResidualLSTM
Bayesian Optimization (GP, Matérn 5/2)	§3.4.2, Table 3	BayesianLSTMTuner, HARLSTMGARCH.tune_lstm
HAR residuals → LSTM	Eqs. 24–28	HARLSTMGARCH.fit
McLeod-Li test	Eq. 26 (Q(20)=4869.18)	metrics.mcleod_li
GARCH(1,1) on hybrid errors	Eqs. 30–33	ResidualGARCH
Benchmarks (GARCH/GJR/LSTM)	Table 4	benchmarks
QLIKE	§3.5	metrics.qlike
Model Confidence Set	§3.5, Table 5	metrics.model_confidence_set

---

Limitations

Mirroring the reference study: results are demonstrated on a single asset class (crude oil); RV is a proxy from daily squared returns (not intraday/jump-robust); the pipeline (HAR + Bayesian optimization + GARCH) is computationally non-trivial; and the LSTM stage remains a "grey box" (attention/SHAP are natural extensions). The framework is univariate — exogenous drivers (OPEC, inventories, geopolitical indices) are not yet included.

---

Citation

If you use this package, please cite the method paper:

@article{BenRomdhane2026HARLSTMGARCH,
  author  = {Ben Romdhane, Wiem and Boubaker, Heni},
  title   = {A Hybrid HAR-LSTM-GARCH Model for Forecasting Volatility in Energy Markets},
  journal = {Journal of Risk and Financial Management},
  year    = {2026},
  volume  = {19},
  number  = {2},
  pages   = {77},
  doi     = {10.3390/jrfm19010077}
}

---

Changelog

0.2.1

• Fix ResidualGARCH.conditional_sigma: fall back to the sample variance when the GARCH fit is (near-)integrated (alpha + beta ≈ 1), instead of dividing by 1 - (alpha + beta) ≈ 0 — prevents an exploding prediction interval at the first observation.
• Sync __version__ with the package metadata.

0.2.0

• Bayesian Optimization of the LSTM (tune.BayesianLSTMTuner, HARLSTMGARCH.tune_lstm) — GP with Matérn 5/2 kernel, no extra dependency.
• Benchmark models (benchmarks): GARCH(1,1)-skewt, GJR-GARCH, standalone LSTM.
• Model Confidence Set test (metrics.model_confidence_set).
• QLIKE corrected to the paper's level-ratio definition; added qlike_loss_series.
• Case study rewritten for full paper reproduction (Brent RV, MCS, benchmarks).
• Added LICENSE, classifiers and project metadata.

0.1.0

• Initial three-stage HAR-LSTM-GARCH implementation.

---

License

MIT — see LICENSE. © 2026 Belhimer Hocine.