Metadata-Version: 2.4
Name: harlstmgarch
Version: 0.2.0
Summary: Three-stage hybrid HAR-LSTM-GARCH framework for realized volatility forecasting
Author-email: Belhimer Hocine <hbelhimer@escf.constantine.dz>
Maintainer-email: Belhimer Hocine <hbelhimer@escf.constantine.dz>
License: MIT
Project-URL: Reference Paper, https://doi.org/10.3390/jrfm19010077
Keywords: volatility forecasting,realized volatility,HAR,LSTM,GARCH,hybrid models,financial econometrics,energy markets,Bayesian optimization,model confidence set
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Intended Audience :: Financial and Insurance Industry
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering
Classifier: Topic :: Office/Business :: Financial :: Investment
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24
Requires-Dist: pandas>=2.0
Requires-Dist: scipy>=1.10
Requires-Dist: statsmodels>=0.14
Requires-Dist: arch>=6.0
Requires-Dist: tensorflow>=2.13
Requires-Dist: scikit-learn>=1.3
Requires-Dist: matplotlib>=3.7
Provides-Extra: dev
Requires-Dist: build>=1.0; extra == "dev"
Requires-Dist: twine>=5.0; extra == "dev"
Requires-Dist: pytest>=7.0; extra == "dev"
Dynamic: license-file

# 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.

[![PyPI](https://img.shields.io/pypi/v/harlstmgarch.svg)](https://pypi.org/project/harlstmgarch/)
[![Python](https://img.shields.io/pypi/pyversions/harlstmgarch.svg)](https://pypi.org/project/harlstmgarch/)
[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](LICENSE)

**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](#the-model)
- [Installation](#installation)
- [Quick start](#quick-start)
- [API reference (detailed syntax)](#api-reference-detailed-syntax)
  - [`harlstmgarch.data`](#harlstmgarchdata)
  - [`harlstmgarch.har`](#harlstmgarchhar)
  - [`harlstmgarch.lstm`](#harlstmgarchlstm)
  - [`harlstmgarch.tune`](#harlstmgarchtune)
  - [`harlstmgarch.garch`](#harlstmgarchgarch)
  - [`harlstmgarch.hybrid`](#harlstmgarchhybrid)
  - [`harlstmgarch.benchmarks`](#harlstmgarchbenchmarks)
  - [`harlstmgarch.metrics`](#harlstmgarchmetrics)
  - [`harlstmgarch.plots`](#harlstmgarchplots)
- [End-to-end example](#end-to-end-example)
- [Reproducing the paper](#reproducing-the-paper)
- [Methodology ↔ paper equations](#methodology--paper-equations)
- [Limitations](#limitations)
- [Citation](#citation)
- [Changelog](#changelog)
- [License](#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

```bash
pip install harlstmgarch
```

With development tools (build / twine / pytest):

```bash
pip install "harlstmgarch[dev]"
```

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

---

## Quick start

```python
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`

```python
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.

```python
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`.

```python
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).

```python
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`

```python
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`

```python
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`

```python
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`

```python
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`

```python
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`

```python
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`

```python
# 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):

```python
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

```python
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/`:

```bash
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:

```bibtex
@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.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](LICENSE). © 2026 Belhimer Hocine.
