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(p, q) with Gaussian or common-df Student-t copula dynamics
• Automatic model selection helpers: AutoVol() and AutoDensity()

Performance

Scivol consistently outperforms other libraries in both speed and optimality. Below are benchmarks on synthetic data. Unless noted, each row uses 1,000,000 observations.

GARCH(1,1)-Normal:

GJR-GARCH(5,5)+Student-t:

AR(1)+EGARCH(2,1)+Normal:

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. Performance
2. Model specification
3. Components
4. Estimation
5. Automatic model selection
6. Results
7. Display settings
8. Diagnostics
9. 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 or Student-t copula DCC dynamics to the standardized residuals.

from scivol import DCC, GARCH, StudentT

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

t_dcc = DCC(1, 1, distribution="studentt")
t_result = t_dcc.fit(returns_df, univariate_spec=GARCH(1, 1) + StudentT())
t_result.params.nu

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
• result.forecast_correlation(horizon) for the forecast correlation path
• result.forecast(horizon) for joint marginal and correlation forecasts

Draw many outcomes from one unchanged forecast origin with scenarios():

scenarios = result.scenarios(50_000, seed=7)
empirical = result.scenarios(50_000, seed=7, innovations="empirical")

The fitted DCC distribution supplies the Gaussian or Student-t copula. The parametric mode maps its draws through each fitted marginal density; the empirical mode uses standardized-residual quantiles. simulate() uses the same copula and marginals along an evolving time path.

Rolling DCC fits warm-start both the marginals and correlation dynamics:

rolling = dcc.fit_rolling(
    returns_df,
    window=1000,
    step=20,
    univariate_spec=GARCH(1, 1) + StudentT(),
    n_jobs=4,
)

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 zero-mean, unit-variance 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:

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

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'

Volatility orders do not implicitly remove observations. Mean models apply the lags they require; hold_back=k requests at least k leading observations be trimmed. Results preserve pandas alignment, and the effective sample is recorded in result.fit_info.

Standalone GARCH, GJR-GARCH, and EGARCH fits accept constrained-space warm starts and explicit optimizer budgets. The same start works in theta- or z-space:

first = spec.fit(returns, maxiter=500)
updated = spec.fit(new_returns, initial_params=first, maxiter=100)
updated.fit_info.warm_started

Current default path policy:

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) + SkewT()
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()
print(len(spec.preview_candidates()))

result = spec.fit(returns, verbose_selection=True)

print(result.spec)
result.selection_summary()
df = result.selection.to_dataframe()

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.

Use a built-in criterion with select_by:

result = spec.fit(returns)                 # "aic+diag" default
result = spec.fit(returns, select_by="aic")
result = spec.fit(returns, select_by="bic")
result = spec.fit(returns, select_by="hqic")

Use OOS selection to fit each candidate on the training portion, score the holdout, then refit the winner on the full sample:

result = spec.fit(returns, select_by="oos")  # log score, last 20%
result = spec.fit(returns, select_by="oos", oos_test_size=500)

oos_loss="qlike" and oos_loss="mse" score variance forecasts, not the full predictive density. They are blocked with AutoDensity unless you pass allow_density_blind_loss=True.

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:

import pandas as pd

spec = GARCH(1, 1) + Normal()
returns = pd.DataFrame(
    {"asset_1": returns1, "asset_2": returns2, "asset_3": returns3}
)
results = spec.fit(returns, n_jobs=4)

for name, r in results.items():
    print(f"{name}: persistence = {r.garch_params.persistence:.4f}")

For repeated portfolio work, fit_many() returns structured successes and failures and supports per-fit limits and prior results:

batch = spec.fit_many(
    returns,
    n_jobs=4,
    maxiter=500,
    timeout=10.0,
    on_error="record",
)

updated = spec.fit_many(new_returns, initial_results=batch, n_jobs=4)
batch.results
batch.failures
batch.warning_summary
batch.to_dataframe()

Rolling windows reuse each series' last successful parameters:

rolling = spec.fit_rolling(returns, window=1000, step=20, n_jobs=4)
latest = rolling[-1].batch

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

spec = GARCH(auto=True) + AutoDensity()
results = spec.fit(returns, n_jobs=4)

Inspecting candidates

result.selection_summary()

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

result.selection.best
result.selection.to_dataframe()

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⁻¹
result.warnings          # persistent fit warnings shown by result.summary()

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.

Every shipped univariate surface has full theta-space and log-space validation against independent AD reference implementations. The machine-checked contract and current verification status live in MODELS.md.

API reference

Exports

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

spec.fit()

result = spec.fit(
    data,                      # 1D array, Series, or DataFrame
    hold_back=0,               # explicit leading observations to trim
    scale=None,                # optional estimation scale
    common_sample=None,        # align auto-selection candidates
    method="mle",              # "mle" or "qmle"
    solver="trust",
    log_mode=True,
    initial_params=None,        # prior result or constrained parameter vector
    maxiter=5000,               # standalone volatility fits
    maxfev=50000,               # derivative-free evaluation budget
    verbose=False,
    n_jobs=None,               # parallel workers (selection or multi-series)
    diagnostic_weight=50.0,    # AIC penalty per failed test
    select_by="aic+diag",      # "aic+diag", "aic", "bic", "hqic", or "oos"
    criterion=None,            # custom scoring callable
    diagnostic_kwargs=None,    # forwarded to diagnostic_tests()
    oos_loss="logscore",       # for select_by="oos"
    oos_test_size=0.2,         # float fraction or int holdout count
)

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

t_dcc = DCC(1, 1, distribution="studentt")

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
result.selection.to_dataframe()

DCCResult

result.params, result.log_likelihood, result.aic, result.bic
result.params.nu              # Student-t DCC
result.Rt
result.corr(0, 1)
result.unconditional_corr
result.std_errors, result.std_errors_robust
result.scenarios(10_000, seed=7)
result.simulate(250, seed=7)
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}")