scivol 0.2.1

Volatility modeling in Python with fast GARCH-family C extensions.

scivol

Fast, accurate volatility models for time series in Python.

scivol is a Python library for likelihood-based volatility modeling with GARCH-family models. It provides composable model specifications such as GARCH(1, 1) + StudentT() and ARMA(1, 1) + GARCH(1, 1) + SkewT(). scivol uses analytical gradients and Hessians for maximum speed and accuracy, and verifies those derivatives against automatic differentiation. It also provides MLE and QMLE estimation, diagnostic tests and automatic model selection.

Supported models

• Univariate volatility models: GARCH(p, q), GJRGARCH(p, q), EGARCH(p, q) (Normal, StudentT, SkewT, and GED surfaces)
• Conditional densities: Normal(), StudentT(), SkewT(), GED()
• Mean model: ARMA(p, q) (standalone with Normal() and GED())
• Composite univariate specifications: volatility-only models and ARMA(p, q) + GARCH(P, Q) with any shipped linked density, plus ARMA(p, q) + EGARCH(P, Q) for Normal(), StudentT(), SkewT(), and GED()
• Exogenous-mean specifications: standalone ARX/HARX with Normal(), StudentT(), SkewT(), and GED(), plus linked ARX/HARX + GARCH(P, Q), ARX/HARX + GJRGARCH(P, Q), and ARX/HARX + EGARCH(P, Q) with the same shipped densities
• Multivariate correlation model: DCC(1, 1) with Gaussian correlation dynamics
• Automatic model selection helpers: AutoVol() and AutoDensity()

Install

pip install scivol

Prebuilt wheels are the default installation path on supported platforms. If no compatible wheel is available, pip falls back to building from the source distribution, which requires a working C toolchain. tqdm ships as a normal dependency, so progress bars are available out of the box.

Quick start

import numpy as np
from scivol import GARCH, Normal, StudentT

np.random.seed(42)
returns = np.random.randn(1000) * 0.01

spec = GARCH(1, 1) + Normal()
result = spec.fit(returns)
result.summary()

For workflows that iterate over many fits, progress bars are enabled by default. You can disable them globally:

import scivol

scivol.settings.show_progress = False

Output:

══════════════════════════════════════════════════════════════════════
                    GARCH Model Estimation Results                    
══════════════════════════════════════════════════════════════════════
Model:       GARCH(1,1)+Normal
Method:      MLE
Date:        2026-01-30 15:42:38
──────────────────────────────────────────────────────────────────────
No. Observations:    1000            Converged:      Yes
No. Parameters:      3               Iterations:     42
Time Elapsed:        0.030s
──────────────────────────────────────────────────────────────────────
Log-Likelihood:            3456.7890
AIC:                      -6907.5780
BIC:                      -6892.8560
──────────────────────────────────────────────────────────────────────

                         Parameter Estimates                          
──────────────────────────────────────────────────────────────────────
Parameter            Coef      Std Err     t-stat      P>|t|
──────────────────────────────────────────────────────────────────────
omega          1.2345e-06    2.34e-07      5.27    <0.001
alpha[1]           0.0523      0.0084      6.23    <0.001
beta[1]            0.9412      0.0092    102.30    <0.001
──────────────────────────────────────────────────────────────────────

                          Model Diagnostics                           
──────────────────────────────────────────────────────────────────────
Persistence (alpha + beta): 0.993500
Stationary:               Yes
Unconditional Variance:   1.900000e-04
Half-life (periods):      106.3
══════════════════════════════════════════════════════════════════════

---

Contents

1. Model specification
2. Components
3. Estimation
4. Automatic model selection
5. Results
6. Display settings
7. Diagnostics
8. API reference

---

Model specification

Build models by combining components with +. scivol orders them automatically: MEAN, then VOLATILITY, then DENSITY.

from scivol import EGARCH, GARCH, GJRGARCH, ARMA, DCC, Normal, StudentT, SkewT

spec = GARCH(1, 1)                          # Normal density by default
spec = GARCH(1, 1) + StudentT()             # Explicit density
spec = GJRGARCH(1, 1) + StudentT()          # Asymmetric volatility
spec = EGARCH(2, 1) + Normal()              # Log-variance volatility model
spec = ARMA(1, 1) + GARCH(1, 1) + SkewT()  # Mean + volatility + density
spec = ARMA(1, 1) + EGARCH(2, 1) + Normal() # Linked mean + EGARCH volatility
dcc = DCC(1, 1)                             # Dynamic correlations (multivariate)

One component per role. If you omit the density, Normal() is added for you.

Alternative operators produce the same result:

spec = garch + normal      # __add__
spec = garch < normal      # __lt__
spec = garch << normal     # __lshift__
spec = normal >> garch     # __rlshift__

---

Components

GARCH(p, q)

Conditional variance:

$$ \sigma_t^2 = \omega + \sum_i \alpha_i \varepsilon_{t-i}^2 + \sum_j \beta_j \sigma_{t-j}^2 $$

Stationary when $\sum \alpha + \sum \beta < 1$.

spec = GARCH(1, 1)  # most common
spec = GARCH(2, 1)

Parameters: omega ($\omega > 0$), alpha[1:p] (ARCH terms), beta[1:q] (GARCH terms).

GJRGARCH(p, q)

Adds a leverage term so that negative shocks raise volatility more than positive shocks of the same size:

$$ \sigma_t^2 = \omega + \sum_i \left(\alpha_i + \gamma_i I(\varepsilon_{t-i} < 0)\right)\varepsilon_{t-i}^2 + \sum_j \beta_j \sigma_{t-j}^2 $$

A negative shock contributes $(\alpha + \gamma)\varepsilon^2$ to the next period's variance; a positive shock contributes $\alpha \varepsilon^2$.

Stationarity:

• Symmetric densities (Normal, Student-t): $\alpha + 0.5\gamma + \beta < 1$
• Asymmetric densities (Skew-t): $\alpha + \gamma P(z < 0) + \beta < 1$

spec = GJRGARCH(1, 1) + StudentT()

EGARCH(p, q)

Log-variance volatility dynamics:

$$ \log h_t = \omega + \sum_{i=1}^{p}\alpha_i(|z_{t-i}| - \mathbb{E}|z|) + \sum_{i=1}^{p}\gamma_i z_{t-i} + \sum_{j=1}^{q}\beta_j \log h_{t-j} $$

spec = EGARCH(1, 1) + Normal()
spec = EGARCH(2, 1) + StudentT()

EGARCH is currently shipped for Normal(), StudentT(), SkewT(), and GED().

Parameters: omega, alpha[1:p], gamma[1:p] (leverage), beta[1:q].

DCC(p, q)

Dynamic Conditional Correlation for multivariate return series.

DCC.fit() uses a two-step workflow:

1. Fit a univariate volatility model to each series.
2. Fit Gaussian DCC dynamics to the resulting standardised residuals.

from scivol import DCC, GARCH, StudentT

dcc = DCC(1, 1)
result = dcc.fit(returns_df, univariate_spec=GARCH(1, 1) + StudentT())

Result access focuses on the economically useful correlation objects:

• result.Rt for the full time-varying correlation path
• result.corr(i, j) for a single pair
• result.unconditional_corr for the long-run correlation matrix

ARMA(p, q)

Conditional mean. Standalone ARMA(p, q) is supported with Normal() and GED(). Linked composite mean-volatility fits currently ship for ARMA(p, q) + GARCH(P, Q) with Normal(), StudentT(), SkewT(), and GED(), and for ARMA(p, q) + GJR-GARCH(P, Q) with Normal(), StudentT(), SkewT(), and GED(), and for ARMA(p, q) + EGARCH(P, Q) with Normal(), StudentT(), SkewT(), and GED().

spec = ARMA(2, 1)                 # standalone ARMA with Normal errors
spec = ARMA(1, 1) + GARCH(1, 1)
spec = ARMA(1, 1) + EGARCH(1, 1) + StudentT()

ARX / HARX

Array-based mean models with explicit regressors via x=. Standalone ARX and HARX now ship with public fit() under the default Normal() surface, and standalone StudentT(), SkewT(), and GED() fits now ship as well. Linked ARX/HARX + GARCH(P, Q), ARX/HARX + GJRGARCH(P, Q), and ARX/HARX + EGARCH(P, Q) fits ship for Normal(), StudentT(), SkewT(), and GED().

from scivol import ARX, HARX, EGARCH, GARCH, GED, StudentT

spec = ARX(1)
result = spec.fit(y, x=x)

spec = HARX((1, 5))
result = spec.fit(y, x=x)

spec = ARX(1) + GARCH(1, 1) + StudentT()
result = spec.fit(y, x=x)

spec = HARX((1, 5)) + GARCH(1, 1) + GED()
result = spec.fit(y, x=x)

spec = ARX(1) + StudentT()
result = spec.fit(y, x=x)

spec = HARX((1, 5)) + EGARCH(1, 1) + GED()
result = spec.fit(y, x=x)

Normal()

Gaussian density. No extra parameters. Default when none is specified.

StudentT()

Heavier tails than Normal. One extra parameter: nu ($\nu > 2$). Lower $\nu$ means fatter tails; as $\nu$ grows, the distribution approaches Normal.

SkewT()

Hansen's skewed Student-t. Two extra parameters: nu ($\nu > 2$) and lambda ($-1 < \lambda < 1$). $\lambda = 0$ gives a symmetric Student-t; $\lambda < 0$ shifts weight to the left tail.

---

Estimation

MLE

The default method.

result = spec.fit(data)

result = spec.fit(
    data,
    solver="trust",
    log_mode=True,
    verbose=True,
)

Solvers:

Solver	Method	Notes
"slsqp"	Sequential quadratic programming	Default; fast and reliable
"nelder-mead"	Derivative-free simplex	Slow but dependable
"trust"	Trust-region (gradient + Hessian)	Fast when it converges
"trust-exact"	Trust-region in log-space	Most stable for difficult data

Log-mode (log_mode=True) transforms constrained parameters into unconstrained space before optimization:

Parameter	Constraint	Transform
$\omega$	$\omega > 0$	softplus(z)
$\alpha, \beta$	$> 0$, sum $< 1$	softmax
$\gamma$ (GJR)	$\gamma > 0$	4-class softmax
$\nu$	$\nu > 2$	2 + softplus(z)
$\lambda$	$-1 < \lambda < 1$	tanh(z)

This guarantees stationarity by construction and avoids boundary problems during optimization.

Default log_mode is benchmark-driven and may differ across model families and densities. If you do not pass log_mode, the effective path is exposed on the fitted result:

result = spec.fit(data)
result.fit_info.to_dict()
# solver='slsqp', log_mode=True/False, optimization_space='z-space'/'theta-space'

Current default path policy:

Family	Default path
GARCH + Normal / StudentT / SkewT	z-space
GJRGARCH + Normal	z-space
GJRGARCH + StudentT / SkewT	theta-space
EGARCH + Normal	theta-space
EGARCH + StudentT	theta-space
EGARCH + SkewT	theta-space
EGARCH + GED	theta-space
ARMA + Normal / GED	theta-space
ARX/HARX + Normal / StudentT / SkewT / GED	theta-space
ARMA + GARCH + Normal / StudentT / SkewT / GED	theta-space
ARMA + EGARCH + Normal / StudentT / SkewT / GED	theta-space
ARX/HARX + GARCH + Normal / StudentT / SkewT / GED	theta-space
ARX/HARX + EGARCH + Normal / StudentT / SkewT / GED	theta-space

For EGARCH + Normal, the current SLSQP policy is based on targeted checks of the shipped EGARCH(1,1) and EGARCH(2,1) surfaces. Those runs showed theta-space matching z-space on fit quality while usually converging in fewer iterations and less wall time.

QMLE

Quasi-maximum likelihood: fit under Normal likelihood, then compute sandwich standard errors valid under distributional misspecification. Pass method='qmle':

spec = GARCH(1, 1) + Normal()
result = spec.fit(data, method='qmle')

result.std_errors        # MLE standard errors
result.std_errors_robust # sandwich (robust) standard errors

This works for standalone GARCH, standalone GJRGARCH, and joint ARMA + GARCH specifications. For joint fits, the robust covariance is computed for the full mean-plus-volatility parameter vector under the Normal QMLE step.

For Student-t or Skew-t, QMLE runs a two-step procedure: first it estimates the mean/volatility parameters under Normal likelihood with sandwich errors, then it fixes those parameters and estimates the distribution parameters by MLE.

from scivol import ARMA

spec = GARCH(1, 1) + StudentT()
result = spec.fit(data, method='qmle')

spec = GJRGARCH(1, 1) + Normal()
result = spec.fit(data, method='qmle')

spec = ARMA(2, 1) + GARCH(1, 2) + Normal()
result = spec.fit(data, method='qmle')

---

Automatic model selection

By GARCH order

Search over lag orders with auto=True:

spec = GARCH(auto=True) + Normal()      # p, q in [1, 3]
spec = GJRGARCH(auto=True) + Normal()

spec = GARCH(auto={'max_p': 2, 'max_q': 2}) + Normal()  # narrower grid
spec = GJRGARCH(p=1, q='auto') + StudentT()              # fix p, search q

By volatility model

AutoVol searches across both GARCH and GJRGARCH families:

from scivol import AutoVol

spec = AutoVol() + Normal()
result = spec.fit(returns)

spec = AutoVol(candidates=['GJRGARCH'], max_p=2, max_q=2) + StudentT()

By distribution

AutoDensity searches across Normal, StudentT, and SkewT:

from scivol import AutoDensity

spec = GARCH(1, 1) + AutoDensity()
spec = GARCH(1, 1) + AutoDensity(candidates=['Normal', 'StudentT'])

Full search

Combine them to search volatility model, order, and distribution at once:

from scivol import AutoVol, AutoDensity

spec = AutoVol() + AutoDensity()
result = spec.fit(returns, verbose_selection=True)

print(result.spec)
result.selection_summary()

Selection criterion

The default score is:

$$ \mathrm{Score} = \mathrm{AIC} + \mathrm{diagnostic_weight} \times n_{\mathrm{failed_tests}} $$

where n_failed_tests counts failures of the DGT and Ljung-Box tests. Default diagnostic_weight is 50.

# Heavier diagnostic penalty
result = spec.fit(returns, diagnostic_weight=100.0)

# AIC only
result = spec.fit(returns, diagnostic_weight=0.0)

For full control, pass a callable:

def my_criterion(result, diagnostics):
    score = result.bic
    if diagnostics is not None:
        if diagnostics['dgt']['p_value'] < 0.01:
            score += 200
        if diagnostics['ljung_box'][1]['reject']:
            score += 100
    return score

spec = AutoVol() + AutoDensity()
result = spec.fit(returns, criterion=my_criterion)

The diagnostics dict matches what result.diagnostic_tests() returns:

{
    'distribution': 'StudentT',
    'dist_params': {'nu': 7.42, 'lam': None},
    'n_obs': 1000,
    'alpha': 0.05,
    'dgt': {
        'n_cells': 40, 'chi2_stat': 34.2,
        'df': 39, 'p_value': 0.689, 'reject': False,
    },
    'ljung_box': {
        1: {'lags': 10, 'q_stat': 8.42, 'p_value': 0.588, 'reject': False},
        2: {'lags': 10, 'q_stat': 7.91, 'p_value': 0.637, 'reject': False},
        3: {'lags': 10, 'q_stat': 11.2, 'p_value': 0.340, 'reject': False},
        4: {'lags': 10, 'q_stat': 9.87, 'p_value': 0.452, 'reject': False},
    },
    'pit': np.ndarray,
}

Pass diagnostic_kwargs to tune test settings:

result = spec.fit(returns, diagnostic_kwargs={'lags': 20, 'n_cells': 50})

When criterion is provided, diagnostic_weight is ignored.

Parallel fitting

Auto-selection fits candidates in parallel by default:

result = spec.fit(returns)           # all cores
result = spec.fit(returns, n_jobs=4) # 4 workers
result = spec.fit(returns, n_jobs=1) # sequential

Multi-series fitting

Fit one specification to many series at once:

spec = GARCH(1, 1) + Normal()
results = spec.fit_multiple([returns1, returns2, returns3], n_jobs=4)

for i, r in enumerate(results):
    print(f"Series {i}: persistence = {r.garch_params.persistence:.4f}")

Auto-selection works here too -- each series gets its own best model:

spec = GARCH(auto=True) + AutoDensity()
results = spec.fit_multiple(returns_list, n_jobs=4)

Inspecting candidates

result.selection_summary()

for c in result._selection_candidates[:5]:
    print(f"{c.spec}: AIC={c.aic:.2f}, Score={c.score:.2f}")

QMLE with AutoDensity is redundant (QMLE always uses Normal likelihood). scivol warns and fits Normal only.

---

Results

spec.fit() returns an EstimationResult.

Parameters

result.params            # flat array: [omega, alpha_1, ..., beta_q, nu?, lambda?]

gp = result.garch_params
gp.omega                 # constant
gp.alpha                 # ARCH coefficients (array)
gp.gamma                 # leverage coefficients (GJR-GARCH only)
gp.beta                  # GARCH coefficients (array)
gp.persistence           # alpha+beta (GARCH) or alpha+0.5*gamma+beta (GJR-GARCH)

dp = result.dist_params
dp.nu                    # degrees of freedom (StudentT, SkewT)
dp.lam                   # skewness (SkewT)

Fit statistics

result.loglikelihood
result.aic
result.bic
result.hqic

Conditional variances and residuals

result.sigma2        # sigma^2_t
result.volatility    # sigma_t
result.std_resid     # epsilon_t / sigma_t

Standard errors

result.std_errors        # MLE (from inverse Hessian)
result.std_errors_robust # sandwich (QMLE only)
result.cov_matrix        # H⁻¹
result.cov_robust        # H⁻¹ @ OPG @ H⁻¹

Output

result.fit_info             # effective solver / theta-vs-z path
result.summary()              # full table
result.summary(robust=True)   # with sandwich SEs
print(result)                 # compact
result.to_dict()              # for programmatic use

---

Display settings

Override how parameter names appear in summaries, print output, and to_dict() keys:

import scivol

scivol.settings.names.gamma = "leverage"
scivol.settings.names.nu = "df"
scivol.settings.names.alpha = "a"

# Now result.summary() shows "leverage[1]" instead of "gamma[1]"
# and result.to_dict() uses "leverage" as a key

Overrides apply to indexed variants too: renaming alpha to a turns alpha[1], alpha[2] into a[1], a[2].

Reset with:

scivol.settings.names.reset()

Available names: omega, alpha, gamma, beta, nu, lambda, const, ar, ma.

Internal attribute names (gp.omega, dp.nu, etc.) never change.

---

Diagnostics

DGT and Ljung-Box tests

Check whether the fitted distribution captures the data:

spec = GARCH(1, 1) + StudentT()
result = spec.fit(returns)
result.diagnostic_tests()

======================================================================
                     Model Diagnostic Tests
======================================================================
Distribution:  StudentT (nu=7.42)
Observations:  1000
Alpha:         0.05

DGT Test (Diebold-Gunther-Tay)
----------------------------------------------------------------------
  Cells:       40         df:          39
  Chi2 stat:   34.20      p-value:     0.6891
  Reject H0:   No (uniform PIT)

Ljung-Box Tests on PIT Moments
----------------------------------------------------------------------
  Moment       Lags       Q-stat      p-value     Reject
  ---------- ------ ------------ ------------ ----------
  (u-0.5)^1     10         8.42       0.5880         No
  (u-0.5)^2     10         7.91       0.6371         No
  (u-0.5)^3     10        11.23       0.3396         No
  (u-0.5)^4     10         9.87       0.4518         No
======================================================================

The DGT test checks whether PIT residuals are uniform -- "No" rejection means the distribution fits. Ljung-Box tests check for serial correlation in PIT moments: moment 1 targets the mean, moment 2 the variance, moments 3--4 skewness and kurtosis.

result.diagnostic_tests(alpha=0.01, n_cells=50, lags=20)

diag = result.diagnostic_tests(print_results=False)
diag['dgt']['p_value']
diag['ljung_box'][2]['reject']

Auto-selection uses these tests internally to penalize poorly fitting candidates.

Analytical gradients and Hessians are verified in the internal development test suite against independent AD reference implementations. That validation machinery is intended for library development rather than end-user workflows.

---

API reference

Exports

from scivol import (
    EGARCH, GARCH, GJRGARCH, ARMA, DCC,
    Normal, StudentT, SkewT,
    AutoDensity, AutoVol,
    Component, CompositeSpec, Role,
    DCCParams, DCCResult,
    settings, __version__,
)

spec.fit()

result = spec.fit(
    data,                      # 1D array, Series, or DataFrame
    method="mle",              # "mle" or "qmle"
    solver="trust",
    log_mode=True,
    verbose=False,
    n_jobs=None,               # parallel workers (auto-selection)
    diagnostic_weight=50.0,    # AIC penalty per failed test
    criterion=None,            # custom scoring callable
    diagnostic_kwargs=None,    # forwarded to diagnostic_tests()
)

GARCH

g = GARCH(1, 1)
g = GARCH(auto=True)
g = GARCH(auto={'max_p': 2, 'max_q': 2})

g.p, g.q, g.n_params, g.signature
g.fitted_params   # {'omega': ..., 'alpha': [...], 'beta': [...]}
g.persistence()
g.is_stationary()
g.unconditional_variance()

GJRGARCH

gjr = GJRGARCH(1, 1)
gjr = GJRGARCH(auto=True)

gjr.n_params              # 1 + 2p + q
gjr.fitted_params         # includes 'gamma'
gjr.persistence()         # alpha + 0.5*gamma + beta
gjr.persistence(p_neg=0.6)

AutoVol

av = AutoVol()
av = AutoVol(candidates=['GJRGARCH'], max_p=2, max_q=2)
av.get_candidates()  # list of (model, p, q) tuples

AutoDensity

ad = AutoDensity()
ad = AutoDensity(candidates=['Normal', 'StudentT'])

DCC

dcc = DCC(1, 1)
dcc.p, dcc.q, dcc.n_params, dcc.signature

EstimationResult

result.params, result.garch_params, result.dist_params
result.loglikelihood, result.aic, result.bic, result.hqic
result.sigma2, result.volatility, result.std_resid
result.std_errors, result.std_errors_robust
result.cov_matrix, result.cov_robust
result.success, result.niter, result.time_elapsed
result.summary(), result.to_dict(), result.diagnostic_tests()
result.selection_summary()
result._selection_candidates  # list of all evaluated models

DCCResult

result.params, result.log_likelihood, result.aic, result.bic
result.Rt
result.corr(0, 1)
result.unconditional_corr
result.std_errors, result.std_errors_robust
result.summary()

---

Examples

Fit GARCH(1,1) and inspect results

import numpy as np
from scivol import GARCH, Normal

returns = np.random.randn(1000) * 0.01

spec = GARCH(1, 1) + Normal()
result = spec.fit(returns)

print(f"Persistence: {result.garch_params.persistence:.4f}")
print(f"Log-likelihood: {result.loglikelihood:.2f}")
result.summary()

Compare MLE and sandwich standard errors

from scivol import GARCH, Normal

spec = GARCH(1, 1) + Normal()
result = spec.fit(returns, method='qmle')

print("Parameter       MLE SE    Robust SE")
print("-" * 40)
for i, name in enumerate(['omega', 'alpha', 'beta']):
    print(f"{name:10}  {result.std_errors[i]:10.6f}  {result.std_errors_robust[i]:10.6f}")

GJR-GARCH leverage effect

from scivol import GJRGARCH, StudentT

spec = GJRGARCH(1, 1) + StudentT()
result = spec.fit(returns)

gp = result.garch_params
print(f"alpha: {gp.alpha[0]:.4f}")
print(f"gamma: {gp.gamma[0]:.4f}")
print(f"beta:  {gp.beta[0]:.4f}")
print(f"nu:    {result.dist_params.nu:.2f}")

if gp.gamma[0] > 0:
    ratio = (gp.alpha[0] + gp.gamma[0]) / gp.alpha[0]
    print(f"Negative shocks hit {ratio:.1f}x harder than positive")

Full automatic search

from scivol import AutoVol, AutoDensity

spec = AutoVol() + AutoDensity()
result = spec.fit(returns, verbose_selection=True)

print(f"Best: {result.spec}")
result.selection_summary()

Fit DCC and inspect correlations

import numpy as np
import pandas as pd
from scivol import DCC, GARCH, Normal

rng = np.random.default_rng(42)
returns = rng.standard_normal(1000) * 0.01
returns_df = pd.DataFrame(
    {
        "stock": returns,
        "bond": returns * 0.3 + rng.standard_normal(len(returns)) * 0.005,
        "commodity": returns * -0.2 + rng.standard_normal(len(returns)) * 0.008,
    }
)

dcc = DCC(1, 1)
result = dcc.fit(returns_df, univariate_spec=GARCH(1, 1) + Normal())

stock_bond = result.corr("stock", "bond")
print(stock_bond.tail())
print(result.unconditional_corr)

One-step-ahead variance forecast

from scivol import GARCH, GJRGARCH, Normal
import numpy as np

# GARCH(1,1)
spec = GARCH(1, 1) + Normal()
result = spec.fit(returns)

gp = result.garch_params
h_next = gp.omega + gp.alpha[0] * returns[-1]**2 + gp.beta[0] * result.sigma2[-1]
print(f"GARCH forecast sigma: {np.sqrt(h_next):.6f}")

# GJR-GARCH(1,1)
spec_gjr = GJRGARCH(1, 1) + Normal()
result_gjr = spec_gjr.fit(returns)

gp = result_gjr.garch_params
indicator = 1.0 if returns[-1] < 0 else 0.0
h_next = (gp.omega
    + gp.alpha[0] * returns[-1]**2
    + gp.gamma[0] * indicator * returns[-1]**2
    + gp.beta[0] * result_gjr.sigma2[-1])
print(f"GJR forecast sigma:   {np.sqrt(h_next):.6f}")