Time Series (wraquant.ts)

Time series analysis and forecasting for financial data: decomposition, seasonality detection, changepoint detection, stationarity tests and transformations, feature extraction, anomaly detection, and forecasting.

Submodules:

• Decomposition – STL, seasonal, SSA, EMD, wavelets, unobserved components

• Forecasting – auto ARIMA, exponential smoothing, Holt-Winters, theta, ensemble, GARCH residual

• Stationarity – ADF, KPSS, Phillips-Perron, fractional differencing, detrending

• Seasonality – detection, Fourier features, multi-seasonal decomposition

• Changepoint – CUSUM, Bayesian changepoint detection

• Anomaly – isolation forest, Prophet-based, Grubbs test

• Features – autocorrelation, spectral, complexity features

• Stochastic – OU, jump-diffusion, regime-switching, and VAR forecasts

• Advanced – tsfresh, stumpy, wavelets, sktime, statsforecast, tslearn, darts

Quick Example

from wraquant.ts import auto_arima, stl_decompose, adf_test

# Test stationarity
adf = adf_test(returns)
print(f"ADF statistic: {adf['statistic']:.4f}")
print(f"p-value: {adf['p_value']:.4f}")
print(f"Stationary: {adf['is_stationary']}")

# STL decomposition (trend + seasonal + residual)
decomp = stl_decompose(prices, period=252)
trend = decomp['trend']
seasonal = decomp['seasonal']
residual = decomp['residual']

# Auto ARIMA forecasting
forecast = auto_arima(returns, h=21)
print(f"21-day forecast mean: {forecast['mean'].mean():.6f}")
print(f"AIC: {forecast['aic']:.2f}")
print(f"Order: {forecast['order']}")

Ensemble Forecasting

from wraquant.ts import ensemble_forecast, forecast_evaluation

# Combine multiple forecasting methods
ensemble = ensemble_forecast(returns, h=21, methods=["arima", "ets", "theta"])
print(f"Ensemble forecast: {ensemble['mean']}")

# Evaluate forecast accuracy with proper metrics
eval_result = forecast_evaluation(actual, predicted)
print(f"MAE: {eval_result['mae']:.6f}")
print(f"RMSE: {eval_result['rmse']:.6f}")
print(f"MAPE: {eval_result['mape']:.2%}")

Stochastic Process Forecasting

from wraquant.ts import ornstein_uhlenbeck_forecast, regime_switching_forecast

# OU process forecast (mean-reverting -- good for spreads)
ou = ornstein_uhlenbeck_forecast(spread, h=21)
print(f"OU mean-reversion speed: {ou['theta']:.4f}")
print(f"Long-run mean: {ou['mu']:.4f}")

# Regime-switching forecast
rs = regime_switching_forecast(returns, n_regimes=2, h=21)
print(f"Forecast (accounting for regime probabilities): {rs['mean']}")

See also

• Statistics (wraquant.stats) – Statistical tests and distribution fitting

• Volatility Modeling (wraquant.vol) – GARCH volatility forecasting

• Regime Detection (wraquant.regimes) – Changepoint and regime detection

API Reference

Time series analysis and forecasting.

Provides a complete pipeline for financial time series modeling: from decomposing a series into trend, seasonal, and residual components, through stationarity testing and transformation, to forecasting with classical, stochastic-process, and machine-learning-backed methods. Designed for the challenges of financial data – non-stationarity, regime changes, fat tails, and low signal-to-noise ratios.

Key sub-modules:

• Decomposition (decomposition) – seasonal_decompose (additive/ multiplicative), stl_decompose (STL – robust to outliers), ssa_decompose (Singular Spectrum Analysis), emd_decompose (Empirical Mode Decomposition), wavelet_decompose, and unobserved_components (Harvey structural model).

• Seasonality (seasonality) – detect_seasonality (auto-detect period), fourier_features (generate sin/cos regressors), multi_fourier_features, seasonal_strength, and multi_seasonal_decompose.

• Change-point detection (changepoint) – detect_changepoints (Bayesian or PELT) and cusum for detecting structural breaks in mean, variance, or trend.

• Stationarity (stationarity) – adf_test, kpss_test, phillips_perron, variance_ratio_test, difference, fractional_difference (preserve memory while achieving stationarity), detrend, and optimal_differencing.

• Features (features) – autocorrelation_features, spectral_features, and complexity_features for summarizing time series structure in ML pipelines.

• Anomaly detection (anomaly) – isolation_forest_ts, prophet_anomaly, and grubbs_test_ts.

• Forecasting (forecasting) – auto_arima, auto_forecast (automatic model selection), exponential_smoothing, theta_forecast, ensemble_forecast, rolling_forecast, holt_winters_forecast, and forecast_evaluation.

• Stochastic processes (stochastic) – ornstein_uhlenbeck_forecast, jump_diffusion_forecast, regime_switching_forecast, and var_forecast.

• Advanced integrations (advanced) – Wrappers for tsfresh, stumpy (matrix profile), tslearn (DTW, k-means), darts, sktime, and statsforecast.

Example

>>> from wraquant.ts import stl_decompose, auto_forecast, adf_test
>>> decomp = stl_decompose(prices, period=252)
>>> stationarity = adf_test(returns)
>>> forecast = auto_forecast(returns, h=20)

Use wraquant.ts for time series modeling and forecasting. For statistical tests on returns (normality, autocorrelation), see wraquant.stats. For volatility-specific forecasting (GARCH), see wraquant.vol. For regime-switching models, see wraquant.regimes.

seasonal_decompose(data, period=None, model='additive')

Decompose a time series into trend, seasonal, and residual components.

Classical decomposition extracts the trend via a centred moving average of width equal to the seasonal period, then computes the seasonal component as the average deviation from the trend for each season.

When to use:

Use classical decomposition for quick EDA and when the seasonal pattern is approximately constant over time. For production forecasting or when outliers are present, prefer stl_decompose which is more robust.

Mathematical formulation:

Additive: y_t = T_t + S_t + R_t Multiplicative: y_t = T_t * S_t * R_t

where T is trend, S is seasonal, and R is residual.

How to interpret:

• result.trend: long-term direction of the series. NaN at endpoints (half the period on each side).

• result.seasonal: periodic pattern that repeats every period observations. Constant across years.

• result.resid: what remains after removing trend and seasonal. Should look like white noise if the decomposition is adequate.

Parameters:

• data (Series) – Time series to decompose.

• period (int | None, default: None) – Seasonal period (e.g., 12 for monthly data with annual seasonality, 252 for daily with annual). Inferred from the index frequency when None.

• model (str, default: 'additive') – "additive" (default) – use when seasonal amplitude is roughly constant. "multiplicative" – use when seasonal amplitude scales with the level.

Return type:

Any

Returns:

DecomposeResult with trend, seasonal, and resid attributes, each a pd.Series.

Example

>>> import pandas as pd, numpy as np
>>> idx = pd.date_range("2020-01-01", periods=120, freq="MS")
>>> data = pd.Series(np.arange(120.0) + 5 * np.sin(np.arange(120) * 2 * np.pi / 12), index=idx)
>>> result = seasonal_decompose(data, period=12)
>>> result.seasonal.iloc[0] != 0
True

See also

stl_decompose: Robust decomposition with time-varying seasonality. trend_filter: Hodrick-Prescott trend extraction.

stl_decompose(data, period=None)

STL (Seasonal and Trend decomposition using Loess) decomposition.

STL uses locally weighted regression (Loess) to extract trend and seasonal components. Unlike classical decomposition, STL allows the seasonal pattern to evolve over time and is robust to outliers via iterative re-weighting.

When to use:

STL is the recommended default decomposition method. Use it when: - The seasonal pattern changes strength or shape over time. - The data contains outliers or level shifts. - You need values at the endpoints (no NaN padding). Prefer classical seasonal_decompose only for quick EDA where simplicity matters more than accuracy.

Mathematical background:

STL applies two nested Loess smoothers iteratively: 1. Inner loop: extract seasonal component via Loess on

sub-series (one per season), then extract trend via Loess on the deseasonalised series.

2. Outer loop: compute robustness weights based on residual magnitude, then re-run the inner loop.

How to interpret:

Same as classical decomposition: trend, seasonal, resid. The key difference is that seasonal can vary over time (plot it to see evolution). The resid should be approximately stationary with no remaining seasonal pattern.

Parameters:

• data (Series) – Time series to decompose.

• period (int | None, default: None) – Seasonal period (e.g., 7 for daily data with weekly seasonality). Uses 7 when None and no index frequency is available.

Return type:

Any

Returns:

STLForecast result with trend, seasonal, and resid attributes, each a pd.Series.

Example

>>> import pandas as pd, numpy as np
>>> idx = pd.date_range("2020-01-01", periods=365, freq="D")
>>> data = pd.Series(np.arange(365.0) + 3 * np.sin(np.arange(365) * 2 * np.pi / 7), index=idx)
>>> result = stl_decompose(data, period=7)
>>> hasattr(result, "trend")
True

See also

seasonal_decompose: Simpler classical decomposition. trend_filter: Hodrick-Prescott trend extraction.

References

• Cleveland et al. (1990), “STL: A Seasonal-Trend Decomposition Procedure Based on Loess”

trend_filter(data, method='hp', lamb=1600)

Extract the trend component from a time series.

Currently supports the Hodrick-Prescott (HP) filter, which separates a time series into a smooth trend and a cyclical component by minimising the sum of squared deviations from the trend subject to a penalty on the trend’s second difference (acceleration).

When to use:

The HP filter is standard in macroeconomic analysis (GDP, unemployment, inflation) where convention dictates its use. For financial time series, prefer STL or moving averages – the HP filter can create spurious cycles and has well-documented endpoint instability (Hamilton 2018).

Mathematical formulation:

min_tau sum_{t=1}^T (y_t - tau_t)^2 +

lambda * sum_{t=3}^T ((tau_t - tau_{t-1}) - (tau_{t-1} - tau_{t-2}))^2

where tau is the trend and lambda controls smoothness. Higher lambda => smoother trend.

How to interpret:

The returned series is the trend component. The cyclical component is data - trend. Common lambda values: - Monthly data: lambda = 129600 (Ravn & Uhlig) - Quarterly data: lambda = 1600 (the Hodrick-Prescott default) - Annual data: lambda = 6.25

Parameters:

• data (Series) – Time series.

• method (str, default: 'hp') – Filter method – "hp" (Hodrick-Prescott, default).

• lamb (float, default: 1600) – Smoothing parameter. Default 1600 (appropriate for quarterly data).

Return type:

Series

Returns:

Trend component as a Series.

Raises:

ValueError – If method is not recognized.

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(100)) + np.arange(100) * 0.5)
>>> trend = trend_filter(data, lamb=1600)
>>> len(trend) == len(data)
True

See also

stl_decompose: Preferred for most decomposition tasks. seasonal_decompose: Classical moving-average decomposition.

References

• Hodrick & Prescott (1997), “Postwar U.S. Business Cycles”

• Hamilton (2018), “Why You Should Never Use the HP Filter”

ssa_decompose(data, window=None, n_components=None, groups=None)

Singular Spectrum Analysis (SSA) decomposition.

SSA embeds a time series into a trajectory (Hankel) matrix, applies SVD to extract orthogonal components, and reconstructs interpretable signal components (trend, oscillatory modes, noise).

SSA vs STL:

• STL assumes a fixed seasonal period that you specify up front. It decomposes into exactly three components: trend, seasonal, residual. Works best when the periodicity is known and stable.

• SSA is data-adaptive and makes no assumption about periodicity. It discovers the dominant oscillatory modes directly from the data via singular values. Better for non-stationary signals, signals with multiple or changing periodicities, or when you want to separate trend from multiple oscillatory components.

Algorithm:

1. Embedding: construct the L x K trajectory matrix (Hankel) where L = window length, K = N - L + 1.

2. SVD: decompose the trajectory matrix into singular triplets.

3. Grouping: assign singular triplets to interpretable groups (trend, oscillatory, noise) either automatically or via user specification.

4. Reconstruction: diagonal averaging (Hankelisation) of each group’s partial matrix to recover the time-domain component.

Parameters:

• data (Series) – Time series to decompose. NaN values are dropped.

• window (int | None, default: None) – Embedding window length L. Should be large enough to capture the longest period of interest but not exceed N/2. If None, defaults to N // 2.

• n_components (int | None, default: None) – Number of leading singular triplets to retain. If None, all are retained.

• groups (dict[str, list[int]] | None, default: None) – Optional grouping of singular triplet indices into named components. For example {"trend": [0], "seasonal": [1, 2]}. Each key maps to a list of 0-based component indices. If None, each singular triplet is returned as its own component ("component_0", "component_1", …).

Returns:

• components: dict mapping component name to reconstructed pd.Series.

• singular_values: 1-D numpy array of singular values.

• explained_variance: 1-D numpy array of explained variance ratio per singular value (sums to 1.0 over all values).

• window: embedding window length used.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> t = np.arange(200, dtype=float)
>>> trend = 0.05 * t
>>> osc = 3 * np.sin(2 * np.pi * t / 20)
>>> noise = np.random.default_rng(42).normal(0, 0.3, 200)
>>> data = pd.Series(trend + osc + noise)
>>> result = ssa_decompose(data, window=40, n_components=5)
>>> # Components sum back to approximately the original
>>> recon = sum(result['components'].values())
>>> np.allclose(recon.values, data.values, atol=1e-8)
True

References

• Golyandina, N. & Zhigljavsky, A. (2013), Singular Spectrum Analysis for Time Series. Springer.

• Vautard, R. & Ghil, M. (1989), “Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series”, Physica D.

emd_decompose(data, max_imfs=10, max_siftings=100, sifting_threshold=0.05)

Empirical Mode Decomposition (EMD) via the sifting algorithm.

EMD adaptively decomposes a non-stationary, non-linear signal into a finite set of Intrinsic Mode Functions (IMFs) and a monotonic residual. Each IMF captures a narrow-band oscillatory mode whose frequency and amplitude can vary over time.

Unlike Fourier or wavelet analysis, EMD makes no assumption about basis functions – the decomposition is entirely data-driven.

Algorithm (sifting):

1. Identify local maxima and minima of the signal.

2. Construct upper and lower envelopes by cubic spline interpolation through the extrema.

3. Compute the mean of the envelopes.

4. Subtract the mean from the signal to obtain a candidate IMF.

5. Repeat steps 1-4 until the candidate satisfies the IMF criteria (symmetric envelopes, near-zero mean).

6. Subtract the extracted IMF from the signal and repeat from step 1 to extract the next IMF.

Parameters:

• data (Series) – Time series to decompose. NaN values are dropped.

• max_imfs (int, default: 10) – Maximum number of IMFs to extract (default 10).

• max_siftings (int, default: 100) – Maximum sifting iterations per IMF (default 100).

• sifting_threshold (float, default: 0.05) – Convergence threshold for the normalised squared difference between successive sifting iterations (default 0.05).

Returns:

• imfs: 2-D numpy array of shape (n_imfs, N) where each row is an IMF ordered from highest to lowest frequency.

• residual: 1-D numpy array of the monotonic residual.

• n_imfs: number of IMFs extracted.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> t = np.linspace(0, 1, 500)
>>> sig = np.sin(2 * np.pi * 5 * t) + np.sin(2 * np.pi * 20 * t) + t
>>> data = pd.Series(sig)
>>> result = emd_decompose(data, max_imfs=5)
>>> result['imfs'].shape[0] >= 1
True
>>> np.all(np.isfinite(result['imfs']))
True

References

• Huang, N.E. et al. (1998), “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis”, Proceedings of the Royal Society A.

wavelet_decompose(data, wavelet='db4', level=None)

Multi-level discrete wavelet decomposition.

Decomposes a time series into approximation and detail coefficients at multiple resolution levels using the Discrete Wavelet Transform (DWT). At each level, the signal is split into a low-frequency approximation and a high-frequency detail component.

When to use:

• Analysing time-frequency structure at multiple scales.

• Denoising (threshold the detail coefficients).

• Feature extraction for ML (energy at each scale).

• When the signal has transient or localised frequency content that Fourier analysis would miss.

Parameters:

• data (Series) – Time series to decompose. NaN values are dropped.

• wavelet (str, default: 'db4') – Wavelet family and order. Common choices: "db4" (Daubechies-4, good general purpose), "haar" (simplest, piecewise constant), "sym5" (near-symmetric Daubechies), "coif3" (Coiflet, near-symmetric with vanishing moments).

• level (int | None, default: None) – Decomposition level. If None, the maximum useful level is computed automatically via pywt.dwt_max_level.

Returns:

• approximation: 1-D numpy array of the final-level approximation coefficients (lowest frequency content).

• details: list of 1-D numpy arrays, one per level, ordered from the finest (highest frequency, level 1) to the coarsest (level level). details[0] is level 1 (finest detail).

• wavelet: wavelet name used.

• level: decomposition level used.

• coeffs: raw coefficient list [cA_n, cD_n, ..., cD_1] as returned by pywt.wavedec for advanced use.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> data = pd.Series(np.random.default_rng(42).normal(0, 1, 256))
>>> result = wavelet_decompose(data, wavelet="db4")
>>> len(result['details']) == result['level']
True

References

• Mallat, S. (2009), A Wavelet Tour of Signal Processing. Academic Press.

unobserved_components(data, level=True, trend=False, cycle=False, seasonal=None, stochastic_level=True, stochastic_trend=True, stochastic_cycle=True, stochastic_seasonal=True)

Unobserved Components Model (UCM) decomposition.

Fits a structural time series model (Harvey 1989) that decomposes the series into interpretable latent components estimated via the Kalman filter and smoother:

• Level (local level / random walk): stochastic intercept.

• Trend (local linear trend): level + stochastic slope.

• Cycle: damped stochastic trigonometric cycle.

• Seasonal: time-varying seasonal pattern.

• Irregular: white noise observation error.

The model is:

y_t = mu_t + gamma_t + c_t + epsilon_t

where mu is the trend component (level + slope), gamma is seasonal, c is cycle, and epsilon is irregular.

When to use:

• Macro-economic or business time series where you want interpretable, probabilistic decomposition.

• When you need confidence intervals on components (STL does not provide these).

• Structural break analysis: sudden changes show up as jumps in the stochastic level.

Parameters:

• data (Series) – Time series to decompose. Should have a DatetimeIndex for best results, but integer index also works.

• level (bool, default: True) – Include a stochastic level component (default True).

• trend (bool, default: False) – Include a stochastic trend / slope component (default False). When True, the model is a local linear trend.

• cycle (bool, default: False) – Include a damped stochastic cycle (default False).

• seasonal (int | None, default: None) – Seasonal period. If None, no seasonal component is included.

• stochastic_level (bool, default: True) – Allow the level to vary over time (default True).

• stochastic_trend (bool, default: True) – Allow the slope to vary over time (default True).

• stochastic_cycle (bool, default: True) – Allow the cycle amplitude/phase to vary (default True).

• stochastic_seasonal (bool, default: True) – Allow the seasonal pattern to evolve (default True).

Returns:

• trend: pd.Series of the estimated trend (level + slope).

• trend_ci: pd.DataFrame with lower and upper columns (95% confidence interval for the trend).

• seasonal: pd.Series of the seasonal component (or None if no seasonal was specified).

• seasonal_ci: pd.DataFrame or None.

• cycle: pd.Series of the cycle component (or None).

• cycle_ci: pd.DataFrame or None.

• irregular: pd.Series of the irregular / residual component.

• model: the fitted statsmodels UnobservedComponentsResults object for further inspection.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> idx = pd.date_range("2015-01-01", periods=120, freq="MS")
>>> t = np.arange(120, dtype=float)
>>> y = 0.5 * t + 3 * np.sin(2 * np.pi * t / 12) + np.random.default_rng(42).normal(0, 1, 120)
>>> data = pd.Series(y, index=idx)
>>> result = unobserved_components(data, level=True, trend=True, seasonal=12)
>>> result['trend'] is not None
True

References

• Harvey, A.C. (1989), Forecasting, Structural Time Series Models and the Kalman Filter. Cambridge University Press.

• Durbin, J. & Koopman, S.J. (2012), Time Series Analysis by State Space Methods. Oxford University Press.

detect_seasonality(data, max_period=365)

Detect dominant seasonal periods via spectral analysis.

Uses Welch’s periodogram to identify significant frequency peaks.

Parameters:

• data (Series) – Time series to analyse.

• max_period (int, default: 365) – Maximum period to consider.

Return type:

list[int]

Returns:

List of detected seasonal periods (in number of observations), sorted by spectral power (strongest first).

fourier_features(index, period, n_harmonics)

Generate Fourier sine/cosine features for a datetime index.

Useful for encoding seasonality as regression features.

Parameters:

• index (DatetimeIndex) – Datetime index.

• period (int) – Seasonal period (in the same time unit as the index).

• n_harmonics (int) – Number of Fourier harmonics to generate.

Return type:

DataFrame

Returns:

DataFrame with sin_k and cos_k columns for each harmonic.

multi_fourier_features(index, periods, n_harmonics=3)

Generate Fourier terms for multiple seasonal periods.

Creates sin/cos pairs for each period-harmonic combination, producing a feature matrix suitable for regression-based seasonal modelling (e.g., linear regression, Prophet-style models, or as exogenous regressors for ARIMA).

For a period P and harmonic k, the features are:

sin(2 * pi * k * t / P) and cos(2 * pi * k * t / P)

This captures up to k-th order seasonal variation within each period.

Parameters:

• index (DatetimeIndex) – Datetime index of the time series.

• periods (list[int]) – List of seasonal periods (in the same units as the index). For example, [7, 365] for daily data with weekly and yearly seasonality.

• n_harmonics (int | list[int], default: 3) – Number of harmonics per period. Can be a single int (applied to all periods) or a list of ints (one per period). Default 3.

Return type:

DataFrame

Returns:

DataFrame with columns named sin_P{period}_H{harmonic} and cos_P{period}_H{harmonic} for each period-harmonic pair. Shape is (len(index), 2 * sum(harmonics)).

Example

>>> import pandas as pd
>>> idx = pd.date_range("2020-01-01", periods=365, freq="D")
>>> df = multi_fourier_features(idx, periods=[7, 365], n_harmonics=3)
>>> df.shape[1]
12

seasonal_strength(data, period=None)

Quantify the strength of seasonality in a time series.

Uses the Wang, Smith, and Hyndman (2006) measure:

F_s = max(0, 1 - Var(R_t) / Var(S_t + R_t))

where S_t is the seasonal component and R_t is the remainder from an STL decomposition.

Returns a float in [0, 1]:

• 0.0: no detectable seasonality (remainder dominates).

• 1.0: perfectly seasonal (no noise).

• > 0.9 is typically considered “strong” seasonality.

• < 0.4 is typically “weak”.

Parameters:

• data (Series) – Time series to evaluate. NaN values are dropped.

• period (int | None, default: None) – Seasonal period for the STL decomposition. If None, defaults to 7.

Return type:

float

Returns:

Strength of seasonality as a float in [0, 1].

Example

>>> import numpy as np, pandas as pd
>>> t = np.arange(200, dtype=float)
>>> pure_seasonal = pd.Series(10 * np.sin(2 * np.pi * t / 20))
>>> strength = seasonal_strength(pure_seasonal, period=20)
>>> strength > 0.9
True

References

• Wang, X., Smith, K. & Hyndman, R. (2006), “Characteristic- Based Clustering for Time Series Data”, Data Mining and Knowledge Discovery.

multi_seasonal_decompose(data, periods)

Decompose a time series with multiple seasonal periods.

Handles complex seasonality (e.g., daily data with both weekly and yearly patterns) by using MSTL (Multiple Seasonal-Trend decomposition using Loess) from statsmodels when available, or an iterative STL approach as fallback.

Parameters:

• data (Series) – Time series to decompose. NaN values are dropped.

• periods (list[int]) – List of seasonal periods, ordered from shortest to longest. For example, [7, 365] for weekly and yearly seasonality in daily data.

Returns:

• trend: pd.Series of the trend component.

• seasonal: dict mapping each period to its pd.Series seasonal component.

• residual: pd.Series of the residual after removing trend and all seasonal components.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> t = np.arange(730, dtype=float)
>>> weekly = 3 * np.sin(2 * np.pi * t / 7)
>>> yearly = 5 * np.sin(2 * np.pi * t / 365)
>>> trend = 0.01 * t
>>> data = pd.Series(trend + weekly + yearly)
>>> result = multi_seasonal_decompose(data, periods=[7, 365])
>>> sorted(result['seasonal'].keys())
[7, 365]

References

• Bandara, K. et al. (2021), “MSTL: A Seasonal-Trend Decomposition Algorithm for Time Series with Multiple Seasonal Patterns”, arXiv:2107.13462.

cusum(data, threshold=1.0)

Cumulative sum (CUSUM) control chart for change detection.

Returns a CUSUM statistic series. Values exceeding threshold standard deviations indicate a potential shift.

Parameters:

• data (Series) – Time series.

• threshold (float, default: 1.0) – Detection threshold in standard deviation units.

Return type:

Series

Returns:

CUSUM statistic series.

detect_changepoints(data, method='pelt', penalty=None)

Detect change-points using the ruptures library.

Parameters:

• data (Series) – Time series.

• method (str, default: 'pelt') – Algorithm — "pelt" (default), "binseg", or "window".

• penalty (float | None, default: None) – Penalty value for the algorithm. When None, a sensible default is chosen.

Return type:

list[int]

Returns:

List of change-point indices (positions in the array, excluding the final point).

Raises:

ValueError – If method is not recognized.

difference(data, order=1)

Apply integer differencing to a time series.

Parameters:

• data (Series) – Time series.

• order (int, default: 1) – Differencing order (1 = first difference).

Return type:

Series

Returns:

Differenced series with NaN values dropped.

fractional_difference(data, d=0.5, threshold=1e-05)

Apply fractional differencing to preserve long-memory information.

Implements the fixed-width window fracdiff method from Advances in Financial Machine Learning (Lopez de Prado).

Parameters:

• data (Series) – Time series.

• d (float, default: 0.5) – Fractional differencing parameter (0 < d < 1).

• threshold (float, default: 1e-05) – Weight cutoff threshold for the window.

Return type:

Series

Returns:

Fractionally differenced series.

detrend(data, method='linear')

Remove trend from a time series.

Parameters:

• data (Series) – Time series.

• method (str, default: 'linear') – Detrending method — "linear" (default) or "constant" (demean).

Return type:

Series

Returns:

Detrended series.

Raises:

ValueError – If method is not recognized.

adf_test(data, max_lags=None, regression='c', significance=0.05)

Augmented Dickey-Fuller (ADF) unit root test.

Tests the null hypothesis that a unit root is present (i.e., the series is non-stationary). A low p-value leads to rejection of the null, concluding stationarity.

The test regression is:

Delta y_t = alpha + beta*t + gamma*y_{t-1} + sum_i delta_i * Delta y_{t-i} + e_t

where alpha is a constant (if regression='c'), beta*t is a time trend (if regression='ct'), and the number of lagged differences is chosen to remove serial correlation in the residuals.

ADF vs KPSS:

• ADF: null = unit root (non-stationary). Rejecting H0 supports stationarity.

• KPSS: null = stationary. Rejecting H0 supports non-stationarity.

• Best practice: run both. If ADF rejects and KPSS does not reject, strong evidence of stationarity. If both reject or both fail to reject, the series may be trend-stationary or difference-stationary.

Parameters:

• data (Series) – Time series to test. NaN values are dropped.

• max_lags (int | None, default: None) – Maximum number of lags for the lagged differences. If None, the optimal lag is selected automatically using the AIC criterion.

• regression (str, default: 'c') – Deterministic terms to include: "c" – constant only (default, most common). "ct" – constant and linear trend. "n" – no constant, no trend.

• significance (float, default: 0.05) – Significance level for the is_stationary convenience flag (default 0.05).

Returns:

• test_statistic: float, the ADF t-statistic.

• p_value: float, MacKinnon approximate p-value.

• optimal_lag: int, number of lags used.

• n_obs: int, number of observations used in the regression.

• critical_values: dict mapping significance levels ("1%", "5%", "10%") to critical values.

• is_stationary: bool, True if p-value < significance.

• interpretation: str, human-readable conclusion.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> # White noise is stationary
>>> stationary = pd.Series(rng.normal(0, 1, 500))
>>> result = adf_test(stationary)
>>> result['is_stationary']
True

References

• Dickey, D.A. & Fuller, W.A. (1979), “Distribution of the Estimators for Autoregressive Time Series With a Unit Root”, JASA.

• Said, S.E. & Dickey, D.A. (1984), “Testing for Unit Roots in Autoregressive-Moving Average Models of Unknown Order”, Biometrika.

kpss_test(data, regression='c', n_lags='auto', significance=0.05)

Kwiatkowski-Phillips-Schmidt-Shin (KPSS) stationarity test.

Tests the null hypothesis that the series is stationary around a deterministic trend. This is the opposite of the ADF test: here, rejecting H0 implies non-stationarity.

When to use ADF vs KPSS vs both:

• ADF alone: quick check; but has low power against near-unit-root alternatives.

• KPSS alone: confirms stationarity; but may over-reject for long-memory processes.

• Both (recommended): a confirmatory strategy. Four cases: 1. ADF rejects, KPSS does not reject -> stationary. 2. ADF does not reject, KPSS rejects -> non-stationary. 3. Both reject -> trend-stationary (difference to remove trend). 4. Neither rejects -> inconclusive, may need more data.

Parameters:

• data (Series) – Time series to test. NaN values are dropped.

• regression (str, default: 'c') – Deterministic component: "c" – test for level stationarity (default). "ct" – test for trend stationarity.

• n_lags (int | str, default: 'auto') – Number of lags for the Newey-West estimator. "auto" (default) uses the data-dependent rule. An integer value fixes the lag truncation.

• significance (float, default: 0.05) – Significance level for the is_stationary convenience flag (default 0.05).

Returns:

• test_statistic: float, the KPSS statistic.

• p_value: float, interpolated p-value (may be bounded at 0.01 or 0.10 depending on the table).

• n_lags: int, number of lags used.

• critical_values: dict mapping significance levels to critical values.

• is_stationary: bool, True if p-value >= significance (i.e., cannot reject the null of stationarity).

• interpretation: str, human-readable conclusion.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> stationary = pd.Series(rng.normal(0, 1, 500))
>>> result = kpss_test(stationary)
>>> result['is_stationary']
True

References

• Kwiatkowski, D. et al. (1992), “Testing the null hypothesis of stationarity against the alternative of a unit root”, Journal of Econometrics.

phillips_perron(data, regression='c', significance=0.05)

Phillips-Perron (PP) unit root test.

Like the ADF test, the PP test has the null hypothesis of a unit root (non-stationarity). However, the PP test uses a non-parametric correction to the t-statistic that is robust to heteroskedasticity and serial correlation in the error term, without needing to specify a lag order.

When to prefer PP over ADF:

• When the error process exhibits heteroskedasticity (e.g., financial returns with volatility clustering).

• When you are uncertain about the appropriate lag length for ADF.

• As a robustness check alongside ADF.

The PP test statistic modifies the ADF statistic with a correction factor based on the Newey-West long-run variance estimate.

Parameters:

• data (Series) – Time series to test. NaN values are dropped.

• regression (str, default: 'c') – Deterministic terms: "c" – constant only (default). "ct" – constant and trend. "n" – no constant, no trend.

• significance (float, default: 0.05) – Significance level for the is_stationary convenience flag (default 0.05).

Returns:

• test_statistic: float, the PP t-statistic.

• p_value: float, MacKinnon approximate p-value.

• n_lags: int, truncation lag for Newey-West.

• critical_values: dict mapping significance levels to critical values.

• is_stationary: bool, True if p-value < significance.

• interpretation: str, human-readable conclusion.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> stationary = pd.Series(rng.normal(0, 1, 500))
>>> result = phillips_perron(stationary)
>>> result['is_stationary']
True

References

• Phillips, P.C.B. & Perron, P. (1988), “Testing for a Unit Root in Time Series Regression”, Biometrika.

optimal_differencing(data, max_d=2, significance=0.05)

Automatically determine the optimal differencing order for stationarity.

Sequentially applies integer differencing (d = 0, 1, 2, …) and runs the ADF test at each order. Returns the smallest d for which the ADF test rejects the unit root null.

Parameters:

• data (Series) – Time series. NaN values are dropped.

• max_d (int, default: 2) – Maximum differencing order to try (default 2). Higher orders are rarely needed in practice.

• significance (float, default: 0.05) – Significance level for the ADF test (default 0.05).

Returns:

• optimal_d: int, the smallest differencing order that achieves stationarity (or max_d if stationarity is not achieved).

• test_results: dict mapping each d to its adf_test result dictionary.

• is_stationary: bool, True if stationarity was achieved within max_d differences.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> # Random walk needs d=1
>>> rw = pd.Series(np.cumsum(rng.normal(0, 1, 500)))
>>> result = optimal_differencing(rw)
>>> result['optimal_d']
1

References

• Hyndman, R.J. & Khandakar, Y. (2008), “Automatic Time Series Forecasting: The forecast Package for R”, JSS.

variance_ratio_test(data, lags=2, overlap=True, significance=0.05)

Lo-MacKinlay variance ratio test for the random walk hypothesis.

Tests whether the variance of k-period returns scales linearly with k, as implied by a random walk. A variance ratio VR(k) = 1 is consistent with a random walk; VR(k) > 1 suggests positive autocorrelation (momentum); VR(k) < 1 suggests mean reversion.

The test statistic under homoskedasticity is:

VR(k) = Var(r_t(k)) / (k * Var(r_t)) z = (VR(k) - 1) / sqrt(2*(2k-1)*(k-1) / (3*k*T))

A heteroskedasticity-robust version is also computed.

Parameters:

• data (Series) – Time series of prices or levels (NOT returns). NaN values are dropped.

• lags (int, default: 2) – Holding period / aggregation interval k (default 2). Common choices: 2, 4, 8, 16.

• overlap (bool, default: True) – Use overlapping returns for better power (default True). Non-overlapping uses fewer observations.

• significance (float, default: 0.05) – Significance level for the is_random_walk convenience flag (default 0.05).

Returns:

• variance_ratio: float, the VR(k) statistic.

• z_statistic: float, the z-statistic (homoskedastic).

• z_robust: float, heteroskedasticity-robust z-statistic.

• p_value: float, two-sided p-value from the robust z-statistic.

• is_random_walk: bool, True if the random walk null cannot be rejected at the given significance level.

• interpretation: str, human-readable conclusion.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> # Random walk: VR should be close to 1
>>> rw = pd.Series(np.cumsum(rng.normal(0, 1, 1000)))
>>> result = variance_ratio_test(rw, lags=2)
>>> 0.5 < result['variance_ratio'] < 1.5
True

References

• Lo, A.W. & MacKinlay, A.C. (1988), “Stock Market Prices Do Not Follow Random Walks: Evidence from a Simple Specification Test”, Review of Financial Studies.

autocorrelation_features(data, n_lags=40, significance=0.05)

Compute autocorrelation and partial autocorrelation features.

Extracts the ACF and PACF values up to n_lags, identifies statistically significant lags, and runs the Ljung-Box test for overall serial correlation.

Interpretation:

• ACF measures linear dependence between y_t and y_{t-k}. For AR(p) processes, ACF decays exponentially. For MA(q), ACF cuts off after lag q.

• PACF measures the correlation between y_t and y_{t-k} after removing the effect of intermediate lags. For AR(p), PACF cuts off after lag p.

• Ljung-Box: tests H0 that the first m autocorrelations are jointly zero. A significant result indicates the series has exploitable temporal structure.

Parameters:

• data (Series) – Time series. NaN values are dropped.

• n_lags (int, default: 40) – Number of lags to compute (default 40).

• significance (float, default: 0.05) – Significance level for identifying significant lags and Ljung-Box test (default 0.05).

Returns:

• acf: 1-D numpy array of ACF values (lags 0 to n_lags).

• pacf: 1-D numpy array of PACF values (lags 0 to n_lags).

• significant_acf_lags: list of lag indices where ACF exceeds the Bartlett confidence band.

• significant_pacf_lags: list of lag indices where PACF exceeds the confidence band.

• ljung_box: dict with statistic, p_value, and is_significant for the Ljung-Box test at the maximum lag.

• first_significant_lag: int or None, the first lag with significant ACF (excluding lag 0).

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> # AR(1) process
>>> x = np.zeros(500)
>>> for i in range(1, 500):
...     x[i] = 0.7 * x[i-1] + rng.normal()
>>> result = autocorrelation_features(pd.Series(x), n_lags=20)
>>> result['acf'][1] > 0.5  # strong lag-1 autocorrelation
True

spectral_features(data, fs=1.0)

Compute power spectral density features.

Extracts summary statistics from the periodogram including the dominant frequency, spectral entropy (a measure of spectral complexity), and spectral flatness (how noise-like the signal is).

Parameters:

• data (Series) – Time series. NaN values are dropped.

• fs (float, default: 1.0) – Sampling frequency (default 1.0). For daily data with annual cycles, use fs=1.0 and interpret frequencies as cycles per observation.

Returns:

• dominant_frequency: float, the frequency with the highest spectral power (excluding DC component).

• dominant_period: float, the period corresponding to the dominant frequency (1 / dominant_frequency).

• spectral_entropy: float, Shannon entropy of the normalised power spectral density. Higher values indicate a more complex, broadband signal; lower values indicate a few dominant frequencies. Range: [0, log2(N/2)].

• spectral_flatness: float in [0, 1]. Ratio of geometric mean to arithmetic mean of the PSD. A value near 1 indicates white noise; near 0 indicates a tonal/periodic signal. Also known as Wiener entropy.

• spectral_centroid: float, the weighted mean frequency (centre of mass of the spectrum).

• psd: 1-D numpy array of power spectral density values.

• frequencies: 1-D numpy array of corresponding frequencies.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> t = np.arange(500, dtype=float)
>>> pure_tone = pd.Series(np.sin(2 * np.pi * t / 20))  # period=20
>>> result = spectral_features(pure_tone)
>>> abs(result['dominant_period'] - 20) < 2
True
>>> result['spectral_flatness'] < 0.3  # tonal signal
True

complexity_features(data, m=2, r=None, order=3)

Compute entropy-based complexity measures for a time series.

Extracts three entropy measures that quantify the regularity and predictability of the signal:

1. Sample entropy (SampEn): probability that patterns similar for m points remain similar for m+1 points. Lower values indicate more regularity (self-similarity); higher values indicate more complexity. Defined as -ln(A/B) where A is the count of m+1 length template matches and B is the count of m length matches, within tolerance r.

2. Approximate entropy (ApEn): similar to SampEn but includes self-matches, making it slightly biased toward regularity. Historically introduced first (Pincus 1991).

3. Permutation entropy (PeEn): based on the frequency of ordinal patterns of length order. Robust to noise and monotonic transformations. Normalised to [0, 1].

Parameters:

• data (Series) – Time series. NaN values are dropped.

• m (int, default: 2) – Embedding dimension for SampEn and ApEn (default 2). Typical values: 2 or 3.

• r (float | None, default: None) – Tolerance threshold for SampEn and ApEn. If None, defaults to 0.2 * std(data) (standard choice).

• order (int, default: 3) – Ordinal pattern length for permutation entropy (default 3). Typical values: 3-7.

Returns:

• sample_entropy: float, SampEn(m, r). Values typically range from 0 (perfectly regular) to ~2.5 (random).

• approximate_entropy: float, ApEn(m, r).

• permutation_entropy: float, normalised PeEn in [0, 1]. 0 = perfectly predictable, 1 = maximally complex/random.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> regular = pd.Series(np.sin(np.arange(500) * 0.1))
>>> random = pd.Series(rng.normal(0, 1, 500))
>>> reg_feat = complexity_features(regular)
>>> rnd_feat = complexity_features(random)
>>> reg_feat['sample_entropy'] < rnd_feat['sample_entropy']
True

References

• Richman, J.S. & Moorman, J.R. (2000), “Physiological time-series analysis using approximate entropy and sample entropy”, AJP.

• Pincus, S.M. (1991), “Approximate entropy as a measure of system complexity”, PNAS.

• Bandt, C. & Pompe, B. (2002), “Permutation Entropy: A Natural Complexity Measure for Time Series”, PRL.

isolation_forest_ts(data, window=20, contamination=0.05, features=None, random_state=42)

Anomaly detection using Isolation Forest on rolling features.

Rather than applying Isolation Forest directly to the raw values, this function first computes rolling window features (mean, std, skew, min, max) to capture the context of each observation. This approach detects contextual anomalies – values that are unusual given their local context, even if they are not globally extreme.

Parameters:

• data (Series) – Time series. NaN values are dropped.

• window (int, default: 20) – Rolling window size for feature computation (default 20).

• contamination (float, default: 0.05) – Expected proportion of anomalies in the data (default 0.05). Controls the threshold for the Isolation Forest decision function.

• features (list[str] | None, default: None) – List of rolling features to compute. If None, uses ["mean", "std", "skew", "min", "max"]. Supported values: "mean", "std", "skew", "min", "max", "median", "range".

• random_state (int, default: 42) – Random seed for reproducibility (default 42).

Returns:

• anomaly_scores: pd.Series of anomaly scores from the Isolation Forest decision function (lower = more anomalous).

• anomaly_labels: pd.Series of int labels (1 = normal, -1 = anomaly).

• threshold: float, the decision function threshold.

• anomaly_indices: list of index values flagged as anomalies.

• n_anomalies: int, number of detected anomalies.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 500)
>>> x[100] = 20.0  # inject anomaly
>>> x[300] = -15.0  # inject anomaly
>>> data = pd.Series(x)
>>> result = isolation_forest_ts(data, window=20)
>>> result['n_anomalies'] > 0
True

References

• Liu, F.T. et al. (2008), “Isolation Forest”, ICDM.

prophet_anomaly(data, k=3.0, seasonal_period=None, trend=True)

Forecast-based anomaly detection using statsmodels.

Fits an Unobserved Components Model (local level + optional seasonal) and flags observations where the actual value deviates from the smoothed state by more than k standard deviations of the residuals.

This approach does NOT require fbprophet or Prophet. It uses pure statsmodels for a fully open-source, lightweight solution.

Parameters:

• data (Series) – Time series. NaN values are dropped.

• k (float, default: 3.0) – Number of standard deviations for the anomaly threshold (default 3.0). Lower values flag more anomalies.

• seasonal_period (int | None, default: None) – Seasonal period for the model. If None, no seasonal component is included.

• trend (bool, default: True) – Include a trend component (default True).

Returns:

• forecast: pd.Series of smoothed / fitted values.

• residuals: pd.Series of (actual - forecast).

• anomaly_mask: pd.Series of bool, True for anomalies.

• anomaly_indices: list of index values flagged as anomalies.

• upper_bound: pd.Series, forecast + k * sigma.

• lower_bound: pd.Series, forecast - k * sigma.

• sigma: float, standard deviation of residuals.

• n_anomalies: int, number of detected anomalies.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 300)
>>> x[50] = 15.0  # anomaly
>>> x[200] = -12.0  # anomaly
>>> data = pd.Series(x)
>>> result = prophet_anomaly(data, k=3.0)
>>> result['n_anomalies'] >= 2
True

References

• Harvey, A.C. (1989), Forecasting, Structural Time Series Models and the Kalman Filter. Cambridge University Press.

grubbs_test_ts(data, window=50, significance=0.05)

Rolling Grubbs test for outlier detection in time series.

Applies the Grubbs test (maximum normed residual test) within a rolling window to detect observations that are statistically unlikely given the local distribution.

The Grubbs test statistic for a sample of size n is:

G = max|x_i - mean| / std

and the critical value is derived from the t-distribution:

G_crit = ((n-1) / sqrt(n)) * sqrt(t^2 / (n - 2 + t^2))

where t is the critical value of the t-distribution with n-2 degrees of freedom at significance alpha / (2n).

Parameters:

• data (Series) – Time series. NaN values are dropped.

• window (int, default: 50) – Rolling window size (default 50). Should be large enough for the Grubbs test to be meaningful (>= 7).

• significance (float, default: 0.05) – Significance level for the Grubbs test (default 0.05).

Returns:

• outlier_mask: pd.Series of bool, True for detected outliers.

• test_statistics: pd.Series of Grubbs G statistics at each position.

• outlier_indices: list of index values flagged as outliers.

• n_outliers: int, number of detected outliers.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 200)
>>> x[75] = 20.0  # obvious outlier
>>> data = pd.Series(x)
>>> result = grubbs_test_ts(data, window=50)
>>> 75 in result['outlier_indices']
True

References

• Grubbs, F.E. (1950), “Sample Criteria for Testing Outlying Observations”, Annals of Mathematical Statistics.

exponential_smoothing(data, **kwargs)

Fit a Holt-Winters exponential smoothing model.

Holt-Winters decomposes the series into level, trend, and seasonal components, updating each with exponentially decaying weights. It is one of the most widely used forecasting methods for business and macro time series.

When to use:

Use exponential smoothing when: - The series has a clear trend and/or seasonal pattern. - You need fast, interpretable forecasts. - The forecast horizon is short (1-2 seasonal cycles). Prefer ARIMA (auto_arima) when the autocorrelation structure is more complex or when you need differencing to achieve stationarity. Prefer ML models (wraquant.ml) for nonlinear patterns.

Mathematical formulation:

Simple: l_t = alpha * y_t + (1 - alpha) * l_{t-1} Double: adds trend b_t = beta * (l_t - l_{t-1}) + (1 - beta) * b_{t-1} Triple: adds seasonal s_t (additive or multiplicative)

How to interpret:

Call .forecast(h) on the returned result for h-step ahead forecasts. Use .summary() to inspect smoothing parameters (alpha, beta, gamma). Low alpha => model relies on history; high alpha => model follows recent data closely.

Parameters:

• data (Series) – Time series to fit. Should have a DatetimeIndex with a set frequency for seasonal models.

• **kwargs (Any) – Keyword arguments forwarded to statsmodels.tsa.holtwinters.ExponentialSmoothing. Key arguments: trend ("add" or "mul"), seasonal ("add" or "mul"), seasonal_periods (int).

Return type:

Any

Returns:

Fitted HoltWintersResultsWrapper. Call .forecast(h) for predictions.

Example

>>> import pandas as pd, numpy as np
>>> idx = pd.date_range("2020-01-01", periods=120, freq="MS")
>>> data = pd.Series(np.arange(120.0) + np.random.randn(120) * 5, index=idx)
>>> result = exponential_smoothing(data, trend="add")
>>> forecast = result.forecast(12)
>>> len(forecast)
12

See also

auto_arima: Automatic ARIMA model selection. wraquant.ts.decomposition.stl_decompose: Decompose before

forecasting.

References

• Holt (1957), “Forecasting seasonals and trends by exponentially weighted moving averages”

• Winters (1960), “Forecasting Sales by Exponentially Weighted Moving Averages”

auto_arima(data, **kwargs)

Automatically select and fit the best ARIMA model.

Performs a stepwise or grid search over ARIMA(p,d,q) orders, selecting the model that minimises AIC (by default). This automates the Box-Jenkins methodology: differencing to achieve stationarity, then selecting AR and MA orders.

When to use:

Use auto_arima when: - You want ARIMA but are unsure of the correct order. - The series has complex autocorrelation (beyond what ETS

handles).

• You need forecast confidence intervals.

Prefer exponential_smoothing for simple trend/seasonal series where speed matters. Prefer arima_model_selection when you want to compare all candidate orders explicitly.

How to interpret:

The returned model object has: - .predict(n_periods) for point forecasts. - .predict_in_sample() for fitted values. - .summary() for model diagnostics (coefficients,

AIC/BIC, residual tests).

• .order() for the selected (p,d,q) order.

Use arima_diagnostics to validate the residuals.

Parameters:

• data (Series) – Time series to fit. Should have >50 observations for reliable order selection.

• **kwargs (Any) – Keyword arguments forwarded to pmdarima.auto_arima. Key arguments include seasonal (bool), m (seasonal period), stepwise (bool, default True for speed), information_criterion ("aic" or "bic").

Return type:

Any

Returns:

Fitted ARIMA model from pmdarima with .predict(), .summary(), and .order() methods.

Raises:

ImportError – If the timeseries optional dependency group (pmdarima) is not installed.

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(200)))
>>> model = auto_arima(data)
>>> model.order()
(1, 1, 0)

See also

exponential_smoothing: Holt-Winters ETS models. arima_diagnostics: Residual checks for fitted ARIMA models. arima_model_selection: Manual grid search over ARIMA orders.

References

• Hyndman & Khandakar (2008), “Automatic Time Series Forecasting: the forecast Package for R”

arima_diagnostics(model_result, nlags=10, alpha=0.05)

Run comprehensive residual diagnostics on a fitted ARIMA model.

After fitting an ARIMA model, it is critical to check that the residuals behave like white noise. This function runs:

• Ljung-Box test: no remaining autocorrelation in residuals.

• Jarque-Bera test: normality of residuals.

• ARCH-LM test: no remaining heteroskedasticity (ARCH effects).

• Durbin-Watson statistic: first-order autocorrelation check.

• ACF/PACF values: for visual inspection of residual structure.

Parameters:

• model_result (Any) – Fitted model result from statsmodels ARIMA/SARIMAX (must have a .resid attribute) or pmdarima ARIMA (must have a .resid() method).

• nlags (int, default: 10) – Number of lags for autocorrelation tests (default 10).

• alpha (float, default: 0.05) – Significance level for pass/fail decisions (default 0.05).

Returns:

• ljung_box: dict with statistic, p_value, pass.

• jarque_bera: dict with statistic, p_value, pass.

• arch_lm: dict with statistic, p_value, pass.

• durbin_watson: Durbin-Watson statistic (near 2 = no autocorrelation).

• acf_values: autocorrelation function values.

• pacf_values: partial autocorrelation function values.

• model_adequate: True if all diagnostic tests pass.

Return type:

dict

Example

>>> from statsmodels.tsa.arima.model import ARIMA
>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(200)))
>>> fit = ARIMA(data, order=(1, 1, 1)).fit()
>>> arima_diagnostics(fit)

arima_model_selection(data, p_range=range(0, 4), d_range=range(0, 3), q_range=range(0, 4), criterion='aic')

Compare ARIMA(p,d,q) combinations and rank by information criterion.

Performs a grid search over combinations of ARIMA orders and ranks models by AIC or BIC. Use this to systematically find the best ARIMA specification for a time series.

Parameters:

• data (Series) – Time series to model.

• p_range (range | list[int], default: range(0, 4)) – Range of AR orders to test (default 0-3).

• d_range (range | list[int], default: range(0, 3)) – Range of differencing orders to test (default 0-2).

• q_range (range | list[int], default: range(0, 4)) – Range of MA orders to test (default 0-3).

• criterion (str, default: 'aic') – Ranking criterion, "aic" (default) or "bic".

Returns:

order (tuple), aic, bic, converged, sorted by the chosen criterion (ascending).

Return type:

DataFrame

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(200)))
>>> arima_model_selection(data, p_range=range(0, 3), d_range=range(0, 2), q_range=range(0, 3))

auto_forecast(data, h=10, holdout_pct=0.2, seasonal_period=None)

Unified auto-forecasting: try multiple models, pick the best.

Fits several forecasting methods on a training subset, evaluates them on a held-out test set, and returns the forecast from the best model re-fitted on the full data. This is the recommended starting point when you have no prior knowledge about the series dynamics.

Models tried:

• SES: Simple Exponential Smoothing.

• Holt-Winters (additive): trend + additive seasonality.

• Theta: Theta method.

• Seasonal naive: repeat the last seasonal cycle.

• Drift: random walk with drift.

When to use auto_forecast vs manual model selection:

Use auto_forecast for rapid prototyping, exploratory work, or when you need a robust default. Switch to manual model selection (arima_model_selection, auto_arima) when you need fine-grained control, interpretability, or domain-specific constraints.

Parameters:

• data (Series) – Time series (at least 30 observations recommended).

• h (int, default: 10) – Forecast horizon (default 10).

• holdout_pct (float, default: 0.2) – Fraction of data reserved for model evaluation (default 0.20).

• seasonal_period (int | None, default: None) – Seasonal period for Holt-Winters and seasonal naive. Auto-detected if None.

Returns:

• best_model: name of the winning model.

• forecast: pd.Series of point forecasts from the best model (re-fitted on full data).

• confidence_intervals: dict with lower and upper pd.Series (approximate, based on residual std).

• model_comparison: pd.DataFrame with RMSE, MAE, MAPE per model.

• residual_diagnostics: dict with mean, std, ljung_box_pvalue of residuals from the best model.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(200)) + 100)
>>> result = auto_forecast(data, h=10)
>>> result['best_model'] in [
...     'ses', 'holt_winters', 'theta', 'seasonal_naive', 'drift',
... ]
True

theta_forecast(data, h=10, theta=2.0)

Theta method forecast (Assimakopoulos & Nikolopoulos, 2000).

The Theta method decomposes the series into two “theta lines”:

• Theta-0: the linear regression line (long-run trend).

• Theta-2: an amplified version of the original curvature, forecast with Simple Exponential Smoothing (SES).

The final forecast is the average of these two components. Despite its simplicity, Theta consistently ranks among the best methods in forecasting competitions (M3, M4).

When to use:

• Quick, competitive baseline for non-seasonal or deseasonalised series.

• When you need a robust method that does not require order selection (unlike ARIMA).

Math:

Given a series y_t, define the second differences D^2(y_t) = y_t - 2*y_{t-1} + y_{t-2}. The theta line Z_theta(t) satisfies D^2(Z) = theta * D^2(y). For theta=0 the solution is a straight line; for theta=2 the curvature is doubled.

Parameters:

• data (Series) – Time series (at least 6 observations).

• h (int, default: 10) – Forecast horizon (default 10).

• theta (float, default: 2.0) – Theta parameter controlling curvature amplification (default 2.0).

Returns:

• forecast: pd.Series of point forecasts.

• fitted_values: pd.Series of in-sample fitted values.

• theta_params: dict with theta, ses_alpha, drift.

Return type:

dict

References

Assimakopoulos, V., & Nikolopoulos, K. (2000). The Theta model: a decomposition approach to forecasting. International Journal of Forecasting, 16(4), 521-530.

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(100)) + 50)
>>> result = theta_forecast(data, h=12)
>>> result['forecast'].shape
(12,)

ses_forecast(data, h=10, alpha=None)

Simple Exponential Smoothing (SES) with optimal alpha.

SES is the simplest exponential smoothing method. It assigns exponentially decreasing weights to past observations. The single smoothing parameter alpha controls how fast old observations are discounted:

• alpha near 1 => recent data dominates (series is volatile).

• alpha near 0 => old data persists (series is stable).

When to use:

• Stationary series with no clear trend or seasonality.

• Quick baseline forecast.

• As a building block inside ensemble or Theta methods.

Math:

l_t = alpha * y_t + (1 - alpha) * l_{t-1}

Forecast: y_{t+h|t} = l_t (flat forecast for all horizons).

Parameters:

• data (Series) – Time series to forecast.

• h (int, default: 10) – Forecast horizon (default 10).

• alpha (float | None, default: None) – Smoothing parameter in (0, 1). If None (default), the optimal alpha is estimated by minimising SSE.

Returns:

• forecast: pd.Series of point forecasts.

• alpha: optimised smoothing parameter.

• fitted_values: pd.Series of in-sample fitted values.

• residuals: pd.Series of in-sample residuals.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.random.randn(100).cumsum() + 50)
>>> result = ses_forecast(data, h=5)
>>> len(result['forecast'])
5

holt_winters_forecast(data, h=10, seasonal='add', seasonal_periods=None, trend='add')

Holt-Winters exponential smoothing with trend and seasonality.

This is the full triple exponential smoothing method. It extends SES by adding a trend component (Holt) and a seasonal component (Winters), supporting both additive and multiplicative seasonality.

When to use:

• Series with clear trend and seasonality (e.g., monthly sales, quarterly GDP).

• Choose seasonal='mul' when seasonal fluctuations grow with the level of the series.

Math (additive):

l_t = alpha * (y_t - s_{t-m}) + (1 - alpha) * (l_{t-1} + b_{t-1}) b_t = beta * (l_t - l_{t-1}) + (1 - beta) * b_{t-1} s_t = gamma * (y_t - l_t) + (1 - gamma) * s_{t-m}

Parameters:

• data (Series) – Time series (must have at least 2 full seasonal cycles).

• h (int, default: 10) – Forecast horizon (default 10).

• seasonal (str, default: 'add') – "add" (default) or "mul" for additive or multiplicative seasonality.

• seasonal_periods (int | None, default: None) – Number of periods in a season. Auto-detected from the index frequency if None.

• trend (str | None, default: 'add') – Trend component type: "add" (default), "mul", or None for no trend.

Returns:

• forecast: pd.Series of point forecasts.

• params: dict with alpha, beta, gamma.

• fitted_values: pd.Series of in-sample fitted values.

• seasonal_components: pd.Series of estimated seasonal factors.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> idx = pd.date_range('2020-01-01', periods=120, freq='MS')
>>> seasonal = 10 * np.sin(2 * np.pi * np.arange(120) / 12)
>>> data = pd.Series(
...     np.arange(120) * 0.5 + seasonal + np.random.randn(120),
...     index=idx,
... )
>>> result = holt_winters_forecast(data, h=12, seasonal_periods=12)
>>> len(result['forecast'])
12

ensemble_forecast(data, h=10, methods=None, holdout_pct=0.2, seasonal_period=None)

Combine multiple forecasting models via inverse-RMSE weighting.

Ensemble forecasting often outperforms any single model because it diversifies model risk. Weights are proportional to the inverse RMSE on a held-out validation set, so better-performing models contribute more.

When to use:

• When no single model clearly dominates.

• For production pipelines where robustness matters more than interpretability.

Parameters:

• data (Series) – Time series to forecast.

• h (int, default: 10) – Forecast horizon (default 10).

• methods (Optional[Sequence[str]], default: None) – List of method names to include. Supported: "ses", "theta", "holt_winters", "drift". If None, uses all four.

• holdout_pct (float, default: 0.2) – Fraction of data for validation (default 0.20).

• seasonal_period (int | None, default: None) – Seasonal period for Holt-Winters. Auto-detected if None.

Returns:

• forecast: pd.Series of weighted-average forecasts.

• weights: dict mapping method name to its weight.

• individual_forecasts: dict mapping method name to pd.Series forecast.

• rmse_per_model: dict mapping method name to its hold-out RMSE.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(200)) + 100)
>>> result = ensemble_forecast(data, h=10)
>>> abs(sum(result['weights'].values()) - 1.0) < 1e-6
True

forecast_evaluation(actual, forecast, naive_forecast=None, benchmark_forecast=None)

Comprehensive forecast accuracy metrics.

Computes a wide set of point-forecast accuracy measures used in the forecasting literature. When a benchmark_forecast is provided the Diebold-Mariano test is run to determine whether the two forecasts have significantly different accuracy.

Metrics included:

• RMSE (Root Mean Squared Error): penalises large errors.

• MAE (Mean Absolute Error): robust to outliers.

• MAPE (Mean Absolute Percentage Error): scale-independent but undefined when actuals are zero.

• SMAPE (Symmetric MAPE): bounded version of MAPE.

• MdAPE (Median Absolute Percentage Error): robust MAPE.

• MASE (Mean Absolute Scaled Error): scaled by in-sample naive error; needs naive_forecast or falls back to MAE / mean(|diff(actual)|).

• Directional accuracy: fraction of periods where the forecast correctly predicts the sign of the change.

• Diebold-Mariano test: compares this forecast against benchmark_forecast (two-sided, squared-error loss).

Parameters:

• actual (Series | ndarray) – Observed values.

• forecast (Series | ndarray) – Predicted values (same length as actual).

• naive_forecast (Series | ndarray | None, default: None) – Optional naive forecast for MASE scaling. If None, MASE is computed using the mean absolute first difference of actual.

• benchmark_forecast (Series | ndarray | None, default: None) – Optional second forecast for Diebold-Mariano comparison.

Returns:

• rmse, mae, mape, smape, mdape, mase, directional_accuracy.

• diebold_mariano: dict with statistic and p_value (only present when benchmark_forecast is given).

Return type:

dict

References

Diebold, F.X. & Mariano, R.S. (1995). Comparing Predictive Accuracy. Journal of Business & Economic Statistics, 13(3).

Example

>>> import numpy as np
>>> actual = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
>>> pred   = np.array([1.1, 2.2, 2.8, 4.3, 4.7])
>>> metrics = forecast_evaluation(actual, pred)
>>> 'rmse' in metrics
True

rolling_forecast(data, forecast_fn, h=1, initial_window=None, step=1)

Walk-forward (expanding window) out-of-sample forecasting.

At each step the model is re-fitted on all data up to time t, then forecasts h steps ahead. The result is a time series of true out-of-sample predictions that can be compared against actuals.

This is the gold standard for evaluating forecasting models because it mimics how the model would be used in production.

Parameters:

• data (Series) – Full time series.

• forecast_fn (Callable[[Series, int], ndarray | Series]) – Callable (train_series, h) -> array-like that returns h point forecasts.

• h (int, default: 1) – Forecast horizon at each step (default 1).

• initial_window (int | None, default: None) – Minimum training window size. Defaults to max(30, len(data) // 3).

• step (int, default: 1) – Number of observations to advance between re-fits (default 1).

Returns:

• forecasts: pd.Series of out-of-sample predictions.

• actuals: pd.Series of corresponding actual values.

• errors: pd.Series of forecast errors.

• cumulative_metrics: dict with cumulative RMSE, MAE.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(100)) + 50)
>>> def my_model(train, h):
...     return np.full(h, train.iloc[-1])
>>> result = rolling_forecast(
...     data, my_model, h=1, initial_window=50,
... )
>>> len(result['forecasts']) > 0
True

garch_residual_forecast(returns, vol_model='GARCH', forecast_method='ses', horizon=10, garch_p=1, garch_q=1)

Forecast returns using GARCH-filtered residuals.

This function chains the vol and ts modules together:

1. Fit a GARCH model (from wraquant.vol.models) to extract conditional volatility and standardised residuals.

2. Forecast the standardised residuals (which should be approximately iid) using a classical time-series method (SES, Theta, or Holt- Winters from wraquant.ts.forecasting).

3. Forecast conditional volatility forward using the GARCH model.

4. Recombine: return_forecast = residual_forecast * vol_forecast.

The key insight is that raw financial returns are not iid (they exhibit volatility clustering). GARCH filtering removes the heteroskedasticity, producing residuals that classical forecasting methods can handle. The recombination step re-introduces the time-varying volatility into the forecast.

Chains: vol.garch_fit -> ts.ses_forecast (or other) -> vol.garch_forecast -> recombine.

Parameters:

• returns (Series) – Return series (daily or intraday).

• vol_model (str, default: 'GARCH') – GARCH variant to use for filtering. Options: "GARCH", "EGARCH", "GJR".

• forecast_method (str, default: 'ses') – Method for forecasting the standardised residuals. Options: "ses" (default), "theta", "holt_winters".

• horizon (int, default: 10) – Number of periods to forecast ahead (default 10).

• garch_p (int, default: 1) – GARCH lag order (default 1).

• garch_q (int, default: 1) – ARCH lag order (default 1).

Returns:

• 'return_forecast' (np.ndarray) – Forecasted returns (residual_forecast * vol_forecast), shape (horizon,).

• 'vol_forecast' (np.ndarray) – GARCH conditional volatility forecast, shape (horizon,).

• 'residual_forecast' (np.ndarray) – Forecasted standardised residuals, shape (horizon,).

• 'standardized_residuals' (np.ndarray) – In-sample standardised residuals from the GARCH fit.

• 'garch_params' (dict) – Fitted GARCH parameters.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.randn(500) * 0.01)
>>> result = garch_residual_forecast(returns, horizon=5)
>>> len(result['return_forecast'])
5
>>> len(result['vol_forecast'])
5

Notes

This approach is related to the “GARCH-in-mean” literature but separates the volatility and level dynamics for modularity.

See also

wraquant.vol.models.garch_fit: GARCH model fitting. wraquant.vol.models.garch_forecast: GARCH volatility forecasting. ses_forecast: SES for residual forecasting. theta_forecast: Theta method for residual forecasting.

ornstein_uhlenbeck_forecast(data, h=10, dt=None, n_simulations=1000, confidence_level=0.95)

Forecast using an Ornstein-Uhlenbeck (OU) mean-reverting process.

The OU process is the continuous-time analogue of an AR(1) process and is the workhorse model for mean-reverting financial quantities.

Use for mean-reverting assets like spreads, interest rates, volatility (VIX), pairs-trading residuals, or any series that fluctuates around a long-run equilibrium.

SDE:

dX_t = theta * (mu - X_t) * dt + sigma * dW_t

where:

• theta > 0: speed of mean reversion (higher = faster).

• mu: long-run mean level.

• sigma: volatility of the diffusion.

• W_t: standard Brownian motion.

The half-life of mean reversion is ln(2) / theta, giving the expected time for a deviation to halve.

Parameter estimation:

Parameters are estimated via Maximum Likelihood on the discrete-time transition density, which is Gaussian: X_{t+dt} | X_t ~ N(mu + (X_t - mu) * exp(-theta*dt), sigma^2 / (2*theta) * (1 - exp(-2*theta*dt)))

Parameters:

• data (Series) – Time series of observed values.

• h (int, default: 10) – Number of steps to forecast (default 10).

• dt (float | None, default: None) – Time step between observations. If None, inferred from the index (1/252 for business days, 1.0 otherwise).

• n_simulations (int, default: 1000) – Number of Monte Carlo paths for confidence bands (default 1000).

• confidence_level (float, default: 0.95) – Confidence level for intervals (default 0.95).

Returns:

• params: dict with theta, mu, sigma.

• half_life: half-life of mean reversion (in observation units).

• forecast: pd.Series of expected path (conditional mean).

• confidence_intervals: dict with lower and upper pd.Series.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> # Simulate OU process
>>> n, theta, mu, sigma = 500, 5.0, 100.0, 2.0
>>> x = np.zeros(n); x[0] = 100.0
>>> dt = 1 / 252
>>> for i in range(1, n):
...     x[i] = x[i-1] + theta * (mu - x[i-1]) * dt + sigma * np.sqrt(dt) * rng.normal()
>>> data = pd.Series(x)
>>> result = ornstein_uhlenbeck_forecast(data, h=20)
>>> result['half_life'] > 0
True

References

Uhlenbeck, G.E. & Ornstein, L.S. (1930). On the theory of Brownian motion. Physical Review, 36(5), 823.

jump_diffusion_forecast(data, h=10, dt=None, n_paths=1000, confidence_level=0.95, seed=None)

Merton jump-diffusion model forecast via Monte Carlo simulation.

Extends geometric Brownian motion (GBM) by adding Poisson-driven jumps, capturing the fat tails and sudden moves observed in real asset returns.

SDE:

dS/S = (mu - lambda*k) * dt + sigma * dW + J * dN

where:

• mu: drift of the diffusion component.

• sigma: diffusion volatility.

• lambda: jump intensity (expected jumps per unit time).

• J: jump size, ln(1+J) ~ N(mu_j, sigma_j^2).

• N_t: Poisson process with intensity lambda.

• k = exp(mu_j + sigma_j^2/2) - 1.

Parameter estimation:

A moment-matching approach is used on log-returns: - Jump intensity estimated from excess kurtosis. - Jump mean and variance from the non-Gaussian tail structure.

Parameters:

• data (Series) – Price series (positive values).

• h (int, default: 10) – Forecast horizon in steps (default 10).

• dt (float | None, default: None) – Time step between observations (default 1/252 for daily).

• n_paths (int, default: 1000) – Number of Monte Carlo simulation paths (default 1000).

• confidence_level (float, default: 0.95) – Confidence level for intervals (default 0.95).

• seed (int | None, default: None) – Random seed for reproducibility.

Returns:

• params: dict with mu, sigma, lambda_, mu_j, sigma_j.

• forecast_paths: np.ndarray of shape (n_paths, h) with simulated price paths.

• mean_forecast: pd.Series of mean across paths.

• confidence_intervals: dict with lower and upper pd.Series.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> prices = pd.Series(100 * np.exp(np.cumsum(np.random.randn(200) * 0.01)))
>>> result = jump_diffusion_forecast(prices, h=10, n_paths=500, seed=42)
>>> result['forecast_paths'].shape
(500, 10)

References

Merton, R.C. (1976). Option pricing when underlying stock returns are discontinuous. Journal of Financial Economics, 3(1-2), 125-144.

regime_switching_forecast(data, n_regimes=2, h=10)

Forecast using a Markov regime-switching model.

Fits a Markov-switching autoregression where the series parameters (mean, variance) switch between discrete hidden regimes. Forecasts are blended across regimes weighted by the predicted regime probabilities.

When to use:

• Series that alternate between distinct states (e.g., bull/bear markets, high/low volatility regimes).

• When a single-model forecast is inadequate because the data generation process changes over time.

Parameters:

• data (Series) – Time series to forecast.

• n_regimes (int, default: 2) – Number of hidden regimes (default 2).

• h (int, default: 10) – Forecast horizon (default 10).

Returns:

• forecast: pd.Series of blended point forecasts.

• regime_forecasts: dict mapping regime index to its pd.Series forecast.

• regime_probs: np.ndarray of shape (h, n_regimes) with predicted regime probabilities.

• blended: same as forecast (for clarity).

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> # Bull/bear regime data
>>> regimes = np.concatenate([rng.normal(0.1, 0.5, 150),
...                           rng.normal(-0.05, 1.0, 150)])
>>> data = pd.Series(np.cumsum(regimes))
>>> result = regime_switching_forecast(data, n_regimes=2, h=10)
>>> result['regime_probs'].shape[1]
2

References

Hamilton, J.D. (1989). A New Approach to the Economic Analysis of Nonstationary Time Series and the Business Cycle. Econometrica, 57(2), 357-384.

var_forecast(data, h=10, maxlags=None, ic='aic')

Vector Autoregression (VAR) for multi-asset forecasting.

VAR models each variable as a linear function of its own lags and the lags of all other variables. This captures lead-lag relationships (e.g., one asset predicting another) and provides impulse response functions (IRF) and forecast error variance decomposition (FEVD).

When to use:

• Forecasting multiple related time series simultaneously (e.g., a portfolio of assets, yield curve factors).

• Analysing dynamic interactions: which variables Granger-cause which.

• Computing impulse responses: how a shock to one variable propagates to others.

Math:

Y_t = c + A_1 * Y_{t-1} + ... + A_p * Y_{t-p} + e_t

where Y_t is a (k x 1) vector, A_i are (k x k) coefficient matrices, and e_t ~ N(0, Sigma).

Parameters:

• data (DataFrame) – DataFrame with multiple columns (one per variable). Should be stationary (difference if needed).

• h (int, default: 10) – Forecast horizon (default 10).

• maxlags (int | None, default: None) – Maximum lag order to consider. If None, uses 12 * (nobs/100)^(1/4) rule of thumb.

• ic (str, default: 'aic') – Information criterion for lag selection: "aic", "bic", "hqic", "fpe" (default "aic").

Returns:

• forecast: pd.DataFrame of point forecasts.

• irf: impulse response function results (dict of DataFrames keyed by shocked variable).

• fevd: forecast error variance decomposition (dict of DataFrames keyed by target variable).

• granger_causality: dict of p-values for Granger causality tests.

• lag_order: selected lag order.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> n = 200
>>> x1 = np.cumsum(rng.normal(0, 1, n))
>>> x2 = 0.3 * np.roll(x1, 1) + rng.normal(0, 1, n)
>>> data = pd.DataFrame({'x1': np.diff(x1), 'x2': np.diff(x2)})
>>> result = var_forecast(data, h=5)
>>> result['forecast'].shape[0]
5

References

Lutkepohl, H. (2005). New Introduction to Multiple Time Series Analysis. Springer.

tsfresh_features(df, column_id='id', column_sort='time')

Extract time series features using tsfresh.

Wraps tsfresh.extract_features with efficient defaults and returns a cleaned DataFrame with no NaN/infinite columns.

Parameters:

• df (DataFrame) – Long-format DataFrame containing time series data. Must include the columns specified by column_id and column_sort, plus one or more value columns.

• column_id (str, default: 'id') – Column identifying distinct time series.

• column_sort (str, default: 'time') – Column used for temporal ordering.

Returns:

Extracted features indexed by column_id values. Columns with NaN or infinite values are dropped.

Return type:

DataFrame

stumpy_matrix_profile(ts, m=50)

Compute the matrix profile of a time series using STUMPY.

The matrix profile identifies the nearest-neighbour distance for every subsequence of length m, enabling motif discovery and discord (anomaly) detection.

Parameters:

• ts (ndarray | Series) – Univariate time series.

• m (int, default: 50) – Subsequence window length.

Returns:

Dictionary containing:

• matrix_profile – 1-D array of nearest-neighbour distances.

• profile_index – 1-D array of indices of the nearest neighbour.

• motif_idx – index of the top motif (lowest MP value).

• discord_idx – index of the top discord (highest MP value).

Return type:

dict[str, Any]

wavelet_transform(data, wavelet='db4', level=None)

Perform discrete wavelet transform (DWT) using PyWavelets.

Parameters:

• data (ndarray | Series) – Univariate signal to decompose.

• wavelet (str, default: 'db4') – Wavelet name (e.g. 'db4', 'haar', 'sym5').

• level (int | None, default: None) – Decomposition level. When None, the maximum useful level is computed automatically.

Returns:

Dictionary containing:

• coeffs – list of coefficient arrays [cA_n, cD_n, ..., cD_1].

• wavelet – wavelet name used.

• level – decomposition level used.

Return type:

dict[str, Any]

wavelet_denoise(data, wavelet='db4', level=None, threshold='soft')

Denoise a signal using wavelet thresholding.

Applies universal (VisuShrink) thresholding to the detail coefficients after a DWT decomposition.

Parameters:

• data (ndarray | Series) – Noisy univariate signal.

• wavelet (str, default: 'db4') – Wavelet name.

• level (int | None, default: None) – Decomposition level (auto-selected when None).

• threshold (str, default: 'soft') – Thresholding mode: 'soft' or 'hard'.

Returns:

Denoised signal with the same length as the input.

Return type:

ndarray

sktime_forecast(y, model='naive', horizon=10)

Forecast a time series using sktime’s unified interface.

Parameters:

• y (Series) – Historical observations indexed by a pandas PeriodIndex or integer index.

• model (str, default: 'naive') –

  Forecasting model. Supported values:

  • 'naive' – last-value forecast

  • 'theta' – Theta method

  • 'ets' – exponential smoothing (AutoETS)

• horizon (int, default: 10) – Number of steps ahead to forecast.

Returns:

DataFrame with columns forecast and index covering the forecast horizon.

Return type:

DataFrame

statsforecast_auto(y, season_length=1, horizon=10)

Automatic forecasting using Nixtla’s StatsForecast.

Runs AutoARIMA and returns point forecasts.

Parameters:

• y (Series) – Historical observations. Must have a DatetimeIndex or integer index.

• season_length (int, default: 1) – Seasonal period length (e.g. 12 for monthly data with yearly seasonality).

• horizon (int, default: 10) – Number of periods to forecast.

Returns:

DataFrame with column forecast indexed over the forecast horizon.

Return type:

DataFrame

tslearn_dtw(ts1, ts2)

Compute DTW distance between two time series using tslearn.

Parameters:

• ts1 (ndarray | Series) – First time series.

• ts2 (ndarray | Series) – Second time series.

Returns:

Dictionary containing:

• distance – DTW distance.

• path – optimal alignment path as a list of (i, j) tuples.

Return type:

dict[str, Any]

tslearn_kmeans(dataset, n_clusters=3)

Cluster time series using DTW-based K-means from tslearn.

Parameters:

• dataset (ndarray | list[ndarray]) – Collection of time series. If a 2-D array, each row is treated as a separate time series. If a 3-D array, shape should be (n_series, n_timestamps, n_features).

• n_clusters (int, default: 3) – Number of clusters.

Returns:

Dictionary containing:

• labels – cluster assignments (1-D int array).

• inertia – sum of distances to nearest cluster centre.

• cluster_centers – array of cluster centre time series.

• n_clusters – number of clusters used.

Return type:

dict[str, Any]

darts_forecast(data, model='nbeats', horizon=10, **kwargs)

Forecast using Darts deep learning models.

Darts provides production-quality implementations of N-BEATS, N-HiTS, Temporal Fusion Transformer, and other deep forecasting models. Use when you need state-of-the-art accuracy and have sufficient data (>1000 observations).

The model is trained on the provided data and produces a forecast of the specified horizon. Default training parameters are tuned for quick experimentation; pass additional keyword arguments to override them (e.g., n_epochs=200).

Parameters:

• data (ndarray | Series) – Price or return series. Must be univariate.

• model (str, default: 'nbeats') –

  Deep learning model to use:

  • 'nbeats' – N-BEATS (Neural Basis Expansion Analysis).

  • 'nhits' – N-HiTS (Neural Hierarchical Interpolation).

  • 'tcn' – Temporal Convolutional Network.

  • 'rnn' – Simple RNN (LSTM-based).

  • 'transformer' – Vanilla Transformer.

• horizon (int, default: 10) – Number of steps ahead to forecast.

• **kwargs (Any) – Additional keyword arguments passed to the model constructor (e.g., n_epochs, input_chunk_length).

Returns:

Dictionary containing:

• forecast – 1-D numpy array of forecasted values.

• model_name – name of the model used.

• training_loss – final training loss (when available).

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.ts.advanced import darts_forecast
>>> data = np.cumsum(np.random.default_rng(0).normal(0, 1, 500))
>>> result = darts_forecast(data, model="nbeats", horizon=5, n_epochs=5)
>>> len(result["forecast"])
5

Notes

Reference: Oreshkin et al. (2020). “N-BEATS: Neural basis expansion analysis for interpretable time series forecasting.” ICLR 2020.

See also

sktime_forecast

Simpler statistical forecasting models.

statsforecast_auto

AutoARIMA and friends.

auto_forecast

wraquant’s built-in automatic forecasting.

Decomposition

Time series decomposition methods.

Decomposition separates a time series into interpretable components – typically trend, seasonal, and residual – enabling cleaner analysis and forecasting. This module provides three approaches:

1. Classical decomposition (seasonal_decompose) – simple moving- average-based extraction of trend and seasonal components. Fast and easy to understand, but uses fixed seasonal patterns and loses data at endpoints.

2. STL decomposition (stl_decompose) – Seasonal and Trend decomposition using Loess. Robust to outliers, allows the seasonal component to vary over time, and preserves all data points. This is the recommended default for most applications.

3. Hodrick-Prescott filter (trend_filter) – extracts a smooth trend by penalising acceleration. Popular in macroeconomics but controversial due to endpoint instability and spurious cycle extraction.

When to use each:

• STL: best general-purpose choice. Handles outliers, time- varying seasonality, and does not lose endpoints. Use this unless you have a specific reason not to.

• Classical decompose: when you need a simple, fast decomposition and are comfortable with fixed seasonality. Good for EDA and presentation.

• HP filter: for macroeconomic trend extraction where convention dictates its use (GDP, unemployment). Avoid for financial time series where it can create spurious cycles.

References

• Cleveland et al. (1990), “STL: A Seasonal-Trend Decomposition Procedure Based on Loess”

• Hodrick & Prescott (1997), “Postwar U.S. Business Cycles: An Empirical Investigation”

• Hamilton (2018), “Why You Should Never Use the Hodrick-Prescott Filter”

seasonal_decompose(data, period=None, model='additive')

Decompose a time series into trend, seasonal, and residual components.

Classical decomposition extracts the trend via a centred moving average of width equal to the seasonal period, then computes the seasonal component as the average deviation from the trend for each season.

When to use:

Use classical decomposition for quick EDA and when the seasonal pattern is approximately constant over time. For production forecasting or when outliers are present, prefer stl_decompose which is more robust.

Mathematical formulation:

Additive: y_t = T_t + S_t + R_t Multiplicative: y_t = T_t * S_t * R_t

where T is trend, S is seasonal, and R is residual.

How to interpret:

• result.trend: long-term direction of the series. NaN at endpoints (half the period on each side).

• result.seasonal: periodic pattern that repeats every period observations. Constant across years.

• result.resid: what remains after removing trend and seasonal. Should look like white noise if the decomposition is adequate.

Parameters:

• data (Series) – Time series to decompose.

• period (int | None, default: None) – Seasonal period (e.g., 12 for monthly data with annual seasonality, 252 for daily with annual). Inferred from the index frequency when None.

• model (str, default: 'additive') – "additive" (default) – use when seasonal amplitude is roughly constant. "multiplicative" – use when seasonal amplitude scales with the level.

Return type:

Any

Returns:

DecomposeResult with trend, seasonal, and resid attributes, each a pd.Series.

Example

>>> import pandas as pd, numpy as np
>>> idx = pd.date_range("2020-01-01", periods=120, freq="MS")
>>> data = pd.Series(np.arange(120.0) + 5 * np.sin(np.arange(120) * 2 * np.pi / 12), index=idx)
>>> result = seasonal_decompose(data, period=12)
>>> result.seasonal.iloc[0] != 0
True

See also

stl_decompose: Robust decomposition with time-varying seasonality. trend_filter: Hodrick-Prescott trend extraction.

stl_decompose(data, period=None)

STL (Seasonal and Trend decomposition using Loess) decomposition.

STL uses locally weighted regression (Loess) to extract trend and seasonal components. Unlike classical decomposition, STL allows the seasonal pattern to evolve over time and is robust to outliers via iterative re-weighting.

When to use:

STL is the recommended default decomposition method. Use it when: - The seasonal pattern changes strength or shape over time. - The data contains outliers or level shifts. - You need values at the endpoints (no NaN padding). Prefer classical seasonal_decompose only for quick EDA where simplicity matters more than accuracy.

Mathematical background:

STL applies two nested Loess smoothers iteratively: 1. Inner loop: extract seasonal component via Loess on

sub-series (one per season), then extract trend via Loess on the deseasonalised series.

2. Outer loop: compute robustness weights based on residual magnitude, then re-run the inner loop.

How to interpret:

Same as classical decomposition: trend, seasonal, resid. The key difference is that seasonal can vary over time (plot it to see evolution). The resid should be approximately stationary with no remaining seasonal pattern.

Parameters:

• data (Series) – Time series to decompose.

• period (int | None, default: None) – Seasonal period (e.g., 7 for daily data with weekly seasonality). Uses 7 when None and no index frequency is available.

Return type:

Any

Returns:

STLForecast result with trend, seasonal, and resid attributes, each a pd.Series.

Example

>>> import pandas as pd, numpy as np
>>> idx = pd.date_range("2020-01-01", periods=365, freq="D")
>>> data = pd.Series(np.arange(365.0) + 3 * np.sin(np.arange(365) * 2 * np.pi / 7), index=idx)
>>> result = stl_decompose(data, period=7)
>>> hasattr(result, "trend")
True

See also

seasonal_decompose: Simpler classical decomposition. trend_filter: Hodrick-Prescott trend extraction.

References

• Cleveland et al. (1990), “STL: A Seasonal-Trend Decomposition Procedure Based on Loess”

trend_filter(data, method='hp', lamb=1600)

Extract the trend component from a time series.

Currently supports the Hodrick-Prescott (HP) filter, which separates a time series into a smooth trend and a cyclical component by minimising the sum of squared deviations from the trend subject to a penalty on the trend’s second difference (acceleration).

When to use:

The HP filter is standard in macroeconomic analysis (GDP, unemployment, inflation) where convention dictates its use. For financial time series, prefer STL or moving averages – the HP filter can create spurious cycles and has well-documented endpoint instability (Hamilton 2018).

Mathematical formulation:

min_tau sum_{t=1}^T (y_t - tau_t)^2 +

lambda * sum_{t=3}^T ((tau_t - tau_{t-1}) - (tau_{t-1} - tau_{t-2}))^2

where tau is the trend and lambda controls smoothness. Higher lambda => smoother trend.

How to interpret:

The returned series is the trend component. The cyclical component is data - trend. Common lambda values: - Monthly data: lambda = 129600 (Ravn & Uhlig) - Quarterly data: lambda = 1600 (the Hodrick-Prescott default) - Annual data: lambda = 6.25

Parameters:

• data (Series) – Time series.

• method (str, default: 'hp') – Filter method – "hp" (Hodrick-Prescott, default).

• lamb (float, default: 1600) – Smoothing parameter. Default 1600 (appropriate for quarterly data).

Return type:

Series

Returns:

Trend component as a Series.

Raises:

ValueError – If method is not recognized.

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(100)) + np.arange(100) * 0.5)
>>> trend = trend_filter(data, lamb=1600)
>>> len(trend) == len(data)
True

See also

stl_decompose: Preferred for most decomposition tasks. seasonal_decompose: Classical moving-average decomposition.

References

• Hodrick & Prescott (1997), “Postwar U.S. Business Cycles”

• Hamilton (2018), “Why You Should Never Use the HP Filter”

ssa_decompose(data, window=None, n_components=None, groups=None)

Singular Spectrum Analysis (SSA) decomposition.

SSA embeds a time series into a trajectory (Hankel) matrix, applies SVD to extract orthogonal components, and reconstructs interpretable signal components (trend, oscillatory modes, noise).

SSA vs STL:

• STL assumes a fixed seasonal period that you specify up front. It decomposes into exactly three components: trend, seasonal, residual. Works best when the periodicity is known and stable.

• SSA is data-adaptive and makes no assumption about periodicity. It discovers the dominant oscillatory modes directly from the data via singular values. Better for non-stationary signals, signals with multiple or changing periodicities, or when you want to separate trend from multiple oscillatory components.

Algorithm:

1. Embedding: construct the L x K trajectory matrix (Hankel) where L = window length, K = N - L + 1.

2. SVD: decompose the trajectory matrix into singular triplets.

3. Grouping: assign singular triplets to interpretable groups (trend, oscillatory, noise) either automatically or via user specification.

4. Reconstruction: diagonal averaging (Hankelisation) of each group’s partial matrix to recover the time-domain component.

Parameters:

• data (Series) – Time series to decompose. NaN values are dropped.

• window (int | None, default: None) – Embedding window length L. Should be large enough to capture the longest period of interest but not exceed N/2. If None, defaults to N // 2.

• n_components (int | None, default: None) – Number of leading singular triplets to retain. If None, all are retained.

• groups (dict[str, list[int]] | None, default: None) – Optional grouping of singular triplet indices into named components. For example {"trend": [0], "seasonal": [1, 2]}. Each key maps to a list of 0-based component indices. If None, each singular triplet is returned as its own component ("component_0", "component_1", …).

Returns:

• components: dict mapping component name to reconstructed pd.Series.

• singular_values: 1-D numpy array of singular values.

• explained_variance: 1-D numpy array of explained variance ratio per singular value (sums to 1.0 over all values).

• window: embedding window length used.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> t = np.arange(200, dtype=float)
>>> trend = 0.05 * t
>>> osc = 3 * np.sin(2 * np.pi * t / 20)
>>> noise = np.random.default_rng(42).normal(0, 0.3, 200)
>>> data = pd.Series(trend + osc + noise)
>>> result = ssa_decompose(data, window=40, n_components=5)
>>> # Components sum back to approximately the original
>>> recon = sum(result['components'].values())
>>> np.allclose(recon.values, data.values, atol=1e-8)
True

References

• Golyandina, N. & Zhigljavsky, A. (2013), Singular Spectrum Analysis for Time Series. Springer.

• Vautard, R. & Ghil, M. (1989), “Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series”, Physica D.

emd_decompose(data, max_imfs=10, max_siftings=100, sifting_threshold=0.05)

Empirical Mode Decomposition (EMD) via the sifting algorithm.

EMD adaptively decomposes a non-stationary, non-linear signal into a finite set of Intrinsic Mode Functions (IMFs) and a monotonic residual. Each IMF captures a narrow-band oscillatory mode whose frequency and amplitude can vary over time.

Unlike Fourier or wavelet analysis, EMD makes no assumption about basis functions – the decomposition is entirely data-driven.

Algorithm (sifting):

1. Identify local maxima and minima of the signal.

2. Construct upper and lower envelopes by cubic spline interpolation through the extrema.

3. Compute the mean of the envelopes.

4. Subtract the mean from the signal to obtain a candidate IMF.

5. Repeat steps 1-4 until the candidate satisfies the IMF criteria (symmetric envelopes, near-zero mean).

6. Subtract the extracted IMF from the signal and repeat from step 1 to extract the next IMF.

Parameters:

• data (Series) – Time series to decompose. NaN values are dropped.

• max_imfs (int, default: 10) – Maximum number of IMFs to extract (default 10).

• max_siftings (int, default: 100) – Maximum sifting iterations per IMF (default 100).

• sifting_threshold (float, default: 0.05) – Convergence threshold for the normalised squared difference between successive sifting iterations (default 0.05).

Returns:

• imfs: 2-D numpy array of shape (n_imfs, N) where each row is an IMF ordered from highest to lowest frequency.

• residual: 1-D numpy array of the monotonic residual.

• n_imfs: number of IMFs extracted.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> t = np.linspace(0, 1, 500)
>>> sig = np.sin(2 * np.pi * 5 * t) + np.sin(2 * np.pi * 20 * t) + t
>>> data = pd.Series(sig)
>>> result = emd_decompose(data, max_imfs=5)
>>> result['imfs'].shape[0] >= 1
True
>>> np.all(np.isfinite(result['imfs']))
True

References

• Huang, N.E. et al. (1998), “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis”, Proceedings of the Royal Society A.

wavelet_decompose(data, wavelet='db4', level=None)

Multi-level discrete wavelet decomposition.

Decomposes a time series into approximation and detail coefficients at multiple resolution levels using the Discrete Wavelet Transform (DWT). At each level, the signal is split into a low-frequency approximation and a high-frequency detail component.

When to use:

• Analysing time-frequency structure at multiple scales.

• Denoising (threshold the detail coefficients).

• Feature extraction for ML (energy at each scale).

• When the signal has transient or localised frequency content that Fourier analysis would miss.

Parameters:

• data (Series) – Time series to decompose. NaN values are dropped.

• wavelet (str, default: 'db4') – Wavelet family and order. Common choices: "db4" (Daubechies-4, good general purpose), "haar" (simplest, piecewise constant), "sym5" (near-symmetric Daubechies), "coif3" (Coiflet, near-symmetric with vanishing moments).

• level (int | None, default: None) – Decomposition level. If None, the maximum useful level is computed automatically via pywt.dwt_max_level.

Returns:

• approximation: 1-D numpy array of the final-level approximation coefficients (lowest frequency content).

• details: list of 1-D numpy arrays, one per level, ordered from the finest (highest frequency, level 1) to the coarsest (level level). details[0] is level 1 (finest detail).

• wavelet: wavelet name used.

• level: decomposition level used.

• coeffs: raw coefficient list [cA_n, cD_n, ..., cD_1] as returned by pywt.wavedec for advanced use.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> data = pd.Series(np.random.default_rng(42).normal(0, 1, 256))
>>> result = wavelet_decompose(data, wavelet="db4")
>>> len(result['details']) == result['level']
True

References

• Mallat, S. (2009), A Wavelet Tour of Signal Processing. Academic Press.

unobserved_components(data, level=True, trend=False, cycle=False, seasonal=None, stochastic_level=True, stochastic_trend=True, stochastic_cycle=True, stochastic_seasonal=True)

Unobserved Components Model (UCM) decomposition.

Fits a structural time series model (Harvey 1989) that decomposes the series into interpretable latent components estimated via the Kalman filter and smoother:

• Level (local level / random walk): stochastic intercept.

• Trend (local linear trend): level + stochastic slope.

• Cycle: damped stochastic trigonometric cycle.

• Seasonal: time-varying seasonal pattern.

• Irregular: white noise observation error.

The model is:

y_t = mu_t + gamma_t + c_t + epsilon_t

where mu is the trend component (level + slope), gamma is seasonal, c is cycle, and epsilon is irregular.

When to use:

• Macro-economic or business time series where you want interpretable, probabilistic decomposition.

• When you need confidence intervals on components (STL does not provide these).

• Structural break analysis: sudden changes show up as jumps in the stochastic level.

Parameters:

• data (Series) – Time series to decompose. Should have a DatetimeIndex for best results, but integer index also works.

• level (bool, default: True) – Include a stochastic level component (default True).

• trend (bool, default: False) – Include a stochastic trend / slope component (default False). When True, the model is a local linear trend.

• cycle (bool, default: False) – Include a damped stochastic cycle (default False).

• seasonal (int | None, default: None) – Seasonal period. If None, no seasonal component is included.

• stochastic_level (bool, default: True) – Allow the level to vary over time (default True).

• stochastic_trend (bool, default: True) – Allow the slope to vary over time (default True).

• stochastic_cycle (bool, default: True) – Allow the cycle amplitude/phase to vary (default True).

• stochastic_seasonal (bool, default: True) – Allow the seasonal pattern to evolve (default True).

Returns:

• trend: pd.Series of the estimated trend (level + slope).

• trend_ci: pd.DataFrame with lower and upper columns (95% confidence interval for the trend).

• seasonal: pd.Series of the seasonal component (or None if no seasonal was specified).

• seasonal_ci: pd.DataFrame or None.

• cycle: pd.Series of the cycle component (or None).

• cycle_ci: pd.DataFrame or None.

• irregular: pd.Series of the irregular / residual component.

• model: the fitted statsmodels UnobservedComponentsResults object for further inspection.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> idx = pd.date_range("2015-01-01", periods=120, freq="MS")
>>> t = np.arange(120, dtype=float)
>>> y = 0.5 * t + 3 * np.sin(2 * np.pi * t / 12) + np.random.default_rng(42).normal(0, 1, 120)
>>> data = pd.Series(y, index=idx)
>>> result = unobserved_components(data, level=True, trend=True, seasonal=12)
>>> result['trend'] is not None
True

References

• Harvey, A.C. (1989), Forecasting, Structural Time Series Models and the Kalman Filter. Cambridge University Press.

• Durbin, J. & Koopman, S.J. (2012), Time Series Analysis by State Space Methods. Oxford University Press.

Forecasting

Time series forecasting wrappers.

Provides classical univariate forecasting models (ARIMA, exponential smoothing), automated model selection, and residual diagnostics.

The two primary approaches are:

1. Exponential Smoothing (ETS / Holt-Winters) – weighted averages of past observations with exponentially decaying weights. Best for series with clear trend and/or seasonal patterns and relatively stable structure. Fast to fit, interpretable, and often competitive with more complex models for short-horizon forecasts.

2. ARIMA / auto_arima – Box-Jenkins methodology that models the series as a linear combination of its own lags (AR), differences (I), and past forecast errors (MA). More flexible than ETS for capturing complex autocorrelation patterns, but requires careful order selection (handled automatically by auto_arima).

When to use each:

• ETS: strong seasonality, smooth trend, short horizon (1-12 steps). Typical for retail sales, energy demand, macro indicators.

• ARIMA: irregular patterns, no obvious seasonality, or when you need confidence intervals that account for parameter uncertainty. auto_arima handles order selection via AIC minimisation.

• Neither: for high-frequency financial returns (near-random- walk), consider GARCH (wraquant.vol) for volatility forecasting or ML models (wraquant.ml).

References

• Hyndman & Athanasopoulos (2021), “Forecasting: Principles and Practice” (3rd ed.)

• Box, Jenkins & Reinsel (2015), “Time Series Analysis”

exponential_smoothing(data, **kwargs)

Fit a Holt-Winters exponential smoothing model.

Holt-Winters decomposes the series into level, trend, and seasonal components, updating each with exponentially decaying weights. It is one of the most widely used forecasting methods for business and macro time series.

When to use:

Use exponential smoothing when: - The series has a clear trend and/or seasonal pattern. - You need fast, interpretable forecasts. - The forecast horizon is short (1-2 seasonal cycles). Prefer ARIMA (auto_arima) when the autocorrelation structure is more complex or when you need differencing to achieve stationarity. Prefer ML models (wraquant.ml) for nonlinear patterns.

Mathematical formulation:

Simple: l_t = alpha * y_t + (1 - alpha) * l_{t-1} Double: adds trend b_t = beta * (l_t - l_{t-1}) + (1 - beta) * b_{t-1} Triple: adds seasonal s_t (additive or multiplicative)

How to interpret:

Call .forecast(h) on the returned result for h-step ahead forecasts. Use .summary() to inspect smoothing parameters (alpha, beta, gamma). Low alpha => model relies on history; high alpha => model follows recent data closely.

Parameters:

• data (Series) – Time series to fit. Should have a DatetimeIndex with a set frequency for seasonal models.

• **kwargs (Any) – Keyword arguments forwarded to statsmodels.tsa.holtwinters.ExponentialSmoothing. Key arguments: trend ("add" or "mul"), seasonal ("add" or "mul"), seasonal_periods (int).

Return type:

Any

Returns:

Fitted HoltWintersResultsWrapper. Call .forecast(h) for predictions.

Example

>>> import pandas as pd, numpy as np
>>> idx = pd.date_range("2020-01-01", periods=120, freq="MS")
>>> data = pd.Series(np.arange(120.0) + np.random.randn(120) * 5, index=idx)
>>> result = exponential_smoothing(data, trend="add")
>>> forecast = result.forecast(12)
>>> len(forecast)
12

See also

auto_arima: Automatic ARIMA model selection. wraquant.ts.decomposition.stl_decompose: Decompose before

forecasting.

References

• Holt (1957), “Forecasting seasonals and trends by exponentially weighted moving averages”

• Winters (1960), “Forecasting Sales by Exponentially Weighted Moving Averages”

auto_arima(data, **kwargs)

Automatically select and fit the best ARIMA model.

Performs a stepwise or grid search over ARIMA(p,d,q) orders, selecting the model that minimises AIC (by default). This automates the Box-Jenkins methodology: differencing to achieve stationarity, then selecting AR and MA orders.

When to use:

Use auto_arima when: - You want ARIMA but are unsure of the correct order. - The series has complex autocorrelation (beyond what ETS

handles).

• You need forecast confidence intervals.

Prefer exponential_smoothing for simple trend/seasonal series where speed matters. Prefer arima_model_selection when you want to compare all candidate orders explicitly.

How to interpret:

The returned model object has: - .predict(n_periods) for point forecasts. - .predict_in_sample() for fitted values. - .summary() for model diagnostics (coefficients,

AIC/BIC, residual tests).

• .order() for the selected (p,d,q) order.

Use arima_diagnostics to validate the residuals.

Parameters:

• data (Series) – Time series to fit. Should have >50 observations for reliable order selection.

• **kwargs (Any) – Keyword arguments forwarded to pmdarima.auto_arima. Key arguments include seasonal (bool), m (seasonal period), stepwise (bool, default True for speed), information_criterion ("aic" or "bic").

Return type:

Any

Returns:

Fitted ARIMA model from pmdarima with .predict(), .summary(), and .order() methods.

Raises:

ImportError – If the timeseries optional dependency group (pmdarima) is not installed.

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(200)))
>>> model = auto_arima(data)
>>> model.order()
(1, 1, 0)

See also

exponential_smoothing: Holt-Winters ETS models. arima_diagnostics: Residual checks for fitted ARIMA models. arima_model_selection: Manual grid search over ARIMA orders.

References

• Hyndman & Khandakar (2008), “Automatic Time Series Forecasting: the forecast Package for R”

arima_diagnostics(model_result, nlags=10, alpha=0.05)

Run comprehensive residual diagnostics on a fitted ARIMA model.

After fitting an ARIMA model, it is critical to check that the residuals behave like white noise. This function runs:

• Ljung-Box test: no remaining autocorrelation in residuals.

• Jarque-Bera test: normality of residuals.

• ARCH-LM test: no remaining heteroskedasticity (ARCH effects).

• Durbin-Watson statistic: first-order autocorrelation check.

• ACF/PACF values: for visual inspection of residual structure.

Parameters:

• model_result (Any) – Fitted model result from statsmodels ARIMA/SARIMAX (must have a .resid attribute) or pmdarima ARIMA (must have a .resid() method).

• nlags (int, default: 10) – Number of lags for autocorrelation tests (default 10).

• alpha (float, default: 0.05) – Significance level for pass/fail decisions (default 0.05).

Returns:

• ljung_box: dict with statistic, p_value, pass.

• jarque_bera: dict with statistic, p_value, pass.

• arch_lm: dict with statistic, p_value, pass.

• durbin_watson: Durbin-Watson statistic (near 2 = no autocorrelation).

• acf_values: autocorrelation function values.

• pacf_values: partial autocorrelation function values.

• model_adequate: True if all diagnostic tests pass.

Return type:

dict

Example

>>> from statsmodels.tsa.arima.model import ARIMA
>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(200)))
>>> fit = ARIMA(data, order=(1, 1, 1)).fit()
>>> arima_diagnostics(fit)

arima_model_selection(data, p_range=range(0, 4), d_range=range(0, 3), q_range=range(0, 4), criterion='aic')

Compare ARIMA(p,d,q) combinations and rank by information criterion.

Performs a grid search over combinations of ARIMA orders and ranks models by AIC or BIC. Use this to systematically find the best ARIMA specification for a time series.

Parameters:

• data (Series) – Time series to model.

• p_range (range | list[int], default: range(0, 4)) – Range of AR orders to test (default 0-3).

• d_range (range | list[int], default: range(0, 3)) – Range of differencing orders to test (default 0-2).

• q_range (range | list[int], default: range(0, 4)) – Range of MA orders to test (default 0-3).

• criterion (str, default: 'aic') – Ranking criterion, "aic" (default) or "bic".

Returns:

order (tuple), aic, bic, converged, sorted by the chosen criterion (ascending).

Return type:

DataFrame

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(200)))
>>> arima_model_selection(data, p_range=range(0, 3), d_range=range(0, 2), q_range=range(0, 3))

theta_forecast(data, h=10, theta=2.0)

Theta method forecast (Assimakopoulos & Nikolopoulos, 2000).

The Theta method decomposes the series into two “theta lines”:

• Theta-0: the linear regression line (long-run trend).

• Theta-2: an amplified version of the original curvature, forecast with Simple Exponential Smoothing (SES).

The final forecast is the average of these two components. Despite its simplicity, Theta consistently ranks among the best methods in forecasting competitions (M3, M4).

When to use:

• Quick, competitive baseline for non-seasonal or deseasonalised series.

• When you need a robust method that does not require order selection (unlike ARIMA).

Math:

Given a series y_t, define the second differences D^2(y_t) = y_t - 2*y_{t-1} + y_{t-2}. The theta line Z_theta(t) satisfies D^2(Z) = theta * D^2(y). For theta=0 the solution is a straight line; for theta=2 the curvature is doubled.

Parameters:

• data (Series) – Time series (at least 6 observations).

• h (int, default: 10) – Forecast horizon (default 10).

• theta (float, default: 2.0) – Theta parameter controlling curvature amplification (default 2.0).

Returns:

• forecast: pd.Series of point forecasts.

• fitted_values: pd.Series of in-sample fitted values.

• theta_params: dict with theta, ses_alpha, drift.

Return type:

dict

References

Assimakopoulos, V., & Nikolopoulos, K. (2000). The Theta model: a decomposition approach to forecasting. International Journal of Forecasting, 16(4), 521-530.

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(100)) + 50)
>>> result = theta_forecast(data, h=12)
>>> result['forecast'].shape
(12,)

ses_forecast(data, h=10, alpha=None)

Simple Exponential Smoothing (SES) with optimal alpha.

SES is the simplest exponential smoothing method. It assigns exponentially decreasing weights to past observations. The single smoothing parameter alpha controls how fast old observations are discounted:

• alpha near 1 => recent data dominates (series is volatile).

• alpha near 0 => old data persists (series is stable).

When to use:

• Stationary series with no clear trend or seasonality.

• Quick baseline forecast.

• As a building block inside ensemble or Theta methods.

Math:

l_t = alpha * y_t + (1 - alpha) * l_{t-1}

Forecast: y_{t+h|t} = l_t (flat forecast for all horizons).

Parameters:

• data (Series) – Time series to forecast.

• h (int, default: 10) – Forecast horizon (default 10).

• alpha (float | None, default: None) – Smoothing parameter in (0, 1). If None (default), the optimal alpha is estimated by minimising SSE.

Returns:

• forecast: pd.Series of point forecasts.

• alpha: optimised smoothing parameter.

• fitted_values: pd.Series of in-sample fitted values.

• residuals: pd.Series of in-sample residuals.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.random.randn(100).cumsum() + 50)
>>> result = ses_forecast(data, h=5)
>>> len(result['forecast'])
5

holt_winters_forecast(data, h=10, seasonal='add', seasonal_periods=None, trend='add')

Holt-Winters exponential smoothing with trend and seasonality.

This is the full triple exponential smoothing method. It extends SES by adding a trend component (Holt) and a seasonal component (Winters), supporting both additive and multiplicative seasonality.

When to use:

• Series with clear trend and seasonality (e.g., monthly sales, quarterly GDP).

• Choose seasonal='mul' when seasonal fluctuations grow with the level of the series.

Math (additive):

l_t = alpha * (y_t - s_{t-m}) + (1 - alpha) * (l_{t-1} + b_{t-1}) b_t = beta * (l_t - l_{t-1}) + (1 - beta) * b_{t-1} s_t = gamma * (y_t - l_t) + (1 - gamma) * s_{t-m}

Parameters:

• data (Series) – Time series (must have at least 2 full seasonal cycles).

• h (int, default: 10) – Forecast horizon (default 10).

• seasonal (str, default: 'add') – "add" (default) or "mul" for additive or multiplicative seasonality.

• seasonal_periods (int | None, default: None) – Number of periods in a season. Auto-detected from the index frequency if None.

• trend (str | None, default: 'add') – Trend component type: "add" (default), "mul", or None for no trend.

Returns:

• forecast: pd.Series of point forecasts.

• params: dict with alpha, beta, gamma.

• fitted_values: pd.Series of in-sample fitted values.

• seasonal_components: pd.Series of estimated seasonal factors.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> idx = pd.date_range('2020-01-01', periods=120, freq='MS')
>>> seasonal = 10 * np.sin(2 * np.pi * np.arange(120) / 12)
>>> data = pd.Series(
...     np.arange(120) * 0.5 + seasonal + np.random.randn(120),
...     index=idx,
... )
>>> result = holt_winters_forecast(data, h=12, seasonal_periods=12)
>>> len(result['forecast'])
12

forecast_evaluation(actual, forecast, naive_forecast=None, benchmark_forecast=None)

Comprehensive forecast accuracy metrics.

Computes a wide set of point-forecast accuracy measures used in the forecasting literature. When a benchmark_forecast is provided the Diebold-Mariano test is run to determine whether the two forecasts have significantly different accuracy.

Metrics included:

• RMSE (Root Mean Squared Error): penalises large errors.

• MAE (Mean Absolute Error): robust to outliers.

• MAPE (Mean Absolute Percentage Error): scale-independent but undefined when actuals are zero.

• SMAPE (Symmetric MAPE): bounded version of MAPE.

• MdAPE (Median Absolute Percentage Error): robust MAPE.

• MASE (Mean Absolute Scaled Error): scaled by in-sample naive error; needs naive_forecast or falls back to MAE / mean(|diff(actual)|).

• Directional accuracy: fraction of periods where the forecast correctly predicts the sign of the change.

• Diebold-Mariano test: compares this forecast against benchmark_forecast (two-sided, squared-error loss).

Parameters:

• actual (Series | ndarray) – Observed values.

• forecast (Series | ndarray) – Predicted values (same length as actual).

• naive_forecast (Series | ndarray | None, default: None) – Optional naive forecast for MASE scaling. If None, MASE is computed using the mean absolute first difference of actual.

• benchmark_forecast (Series | ndarray | None, default: None) – Optional second forecast for Diebold-Mariano comparison.

Returns:

• rmse, mae, mape, smape, mdape, mase, directional_accuracy.

• diebold_mariano: dict with statistic and p_value (only present when benchmark_forecast is given).

Return type:

dict

References

Diebold, F.X. & Mariano, R.S. (1995). Comparing Predictive Accuracy. Journal of Business & Economic Statistics, 13(3).

Example

>>> import numpy as np
>>> actual = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
>>> pred   = np.array([1.1, 2.2, 2.8, 4.3, 4.7])
>>> metrics = forecast_evaluation(actual, pred)
>>> 'rmse' in metrics
True

auto_forecast(data, h=10, holdout_pct=0.2, seasonal_period=None)

Unified auto-forecasting: try multiple models, pick the best.

Fits several forecasting methods on a training subset, evaluates them on a held-out test set, and returns the forecast from the best model re-fitted on the full data. This is the recommended starting point when you have no prior knowledge about the series dynamics.

Models tried:

• SES: Simple Exponential Smoothing.

• Holt-Winters (additive): trend + additive seasonality.

• Theta: Theta method.

• Seasonal naive: repeat the last seasonal cycle.

• Drift: random walk with drift.

When to use auto_forecast vs manual model selection:

Use auto_forecast for rapid prototyping, exploratory work, or when you need a robust default. Switch to manual model selection (arima_model_selection, auto_arima) when you need fine-grained control, interpretability, or domain-specific constraints.

Parameters:

• data (Series) – Time series (at least 30 observations recommended).

• h (int, default: 10) – Forecast horizon (default 10).

• holdout_pct (float, default: 0.2) – Fraction of data reserved for model evaluation (default 0.20).

• seasonal_period (int | None, default: None) – Seasonal period for Holt-Winters and seasonal naive. Auto-detected if None.

Returns:

• best_model: name of the winning model.

• forecast: pd.Series of point forecasts from the best model (re-fitted on full data).

• confidence_intervals: dict with lower and upper pd.Series (approximate, based on residual std).

• model_comparison: pd.DataFrame with RMSE, MAE, MAPE per model.

• residual_diagnostics: dict with mean, std, ljung_box_pvalue of residuals from the best model.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(200)) + 100)
>>> result = auto_forecast(data, h=10)
>>> result['best_model'] in [
...     'ses', 'holt_winters', 'theta', 'seasonal_naive', 'drift',
... ]
True

ensemble_forecast(data, h=10, methods=None, holdout_pct=0.2, seasonal_period=None)

Combine multiple forecasting models via inverse-RMSE weighting.

Ensemble forecasting often outperforms any single model because it diversifies model risk. Weights are proportional to the inverse RMSE on a held-out validation set, so better-performing models contribute more.

When to use:

• When no single model clearly dominates.

• For production pipelines where robustness matters more than interpretability.

Parameters:

• data (Series) – Time series to forecast.

• h (int, default: 10) – Forecast horizon (default 10).

• methods (Optional[Sequence[str]], default: None) – List of method names to include. Supported: "ses", "theta", "holt_winters", "drift". If None, uses all four.

• holdout_pct (float, default: 0.2) – Fraction of data for validation (default 0.20).

• seasonal_period (int | None, default: None) – Seasonal period for Holt-Winters. Auto-detected if None.

Returns:

• forecast: pd.Series of weighted-average forecasts.

• weights: dict mapping method name to its weight.

• individual_forecasts: dict mapping method name to pd.Series forecast.

• rmse_per_model: dict mapping method name to its hold-out RMSE.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(200)) + 100)
>>> result = ensemble_forecast(data, h=10)
>>> abs(sum(result['weights'].values()) - 1.0) < 1e-6
True

rolling_forecast(data, forecast_fn, h=1, initial_window=None, step=1)

Walk-forward (expanding window) out-of-sample forecasting.

At each step the model is re-fitted on all data up to time t, then forecasts h steps ahead. The result is a time series of true out-of-sample predictions that can be compared against actuals.

This is the gold standard for evaluating forecasting models because it mimics how the model would be used in production.

Parameters:

• data (Series) – Full time series.

• forecast_fn (Callable[[Series, int], ndarray | Series]) – Callable (train_series, h) -> array-like that returns h point forecasts.

• h (int, default: 1) – Forecast horizon at each step (default 1).

• initial_window (int | None, default: None) – Minimum training window size. Defaults to max(30, len(data) // 3).

• step (int, default: 1) – Number of observations to advance between re-fits (default 1).

Returns:

• forecasts: pd.Series of out-of-sample predictions.

• actuals: pd.Series of corresponding actual values.

• errors: pd.Series of forecast errors.

• cumulative_metrics: dict with cumulative RMSE, MAE.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> data = pd.Series(np.cumsum(np.random.randn(100)) + 50)
>>> def my_model(train, h):
...     return np.full(h, train.iloc[-1])
>>> result = rolling_forecast(
...     data, my_model, h=1, initial_window=50,
... )
>>> len(result['forecasts']) > 0
True

garch_residual_forecast(returns, vol_model='GARCH', forecast_method='ses', horizon=10, garch_p=1, garch_q=1)

Forecast returns using GARCH-filtered residuals.

This function chains the vol and ts modules together:

1. Fit a GARCH model (from wraquant.vol.models) to extract conditional volatility and standardised residuals.

2. Forecast the standardised residuals (which should be approximately iid) using a classical time-series method (SES, Theta, or Holt- Winters from wraquant.ts.forecasting).

3. Forecast conditional volatility forward using the GARCH model.

4. Recombine: return_forecast = residual_forecast * vol_forecast.

The key insight is that raw financial returns are not iid (they exhibit volatility clustering). GARCH filtering removes the heteroskedasticity, producing residuals that classical forecasting methods can handle. The recombination step re-introduces the time-varying volatility into the forecast.

Chains: vol.garch_fit -> ts.ses_forecast (or other) -> vol.garch_forecast -> recombine.

Parameters:

• returns (Series) – Return series (daily or intraday).

• vol_model (str, default: 'GARCH') – GARCH variant to use for filtering. Options: "GARCH", "EGARCH", "GJR".

• forecast_method (str, default: 'ses') – Method for forecasting the standardised residuals. Options: "ses" (default), "theta", "holt_winters".

• horizon (int, default: 10) – Number of periods to forecast ahead (default 10).

• garch_p (int, default: 1) – GARCH lag order (default 1).

• garch_q (int, default: 1) – ARCH lag order (default 1).

Returns:

• 'return_forecast' (np.ndarray) – Forecasted returns (residual_forecast * vol_forecast), shape (horizon,).

• 'vol_forecast' (np.ndarray) – GARCH conditional volatility forecast, shape (horizon,).

• 'residual_forecast' (np.ndarray) – Forecasted standardised residuals, shape (horizon,).

• 'standardized_residuals' (np.ndarray) – In-sample standardised residuals from the GARCH fit.

• 'garch_params' (dict) – Fitted GARCH parameters.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.randn(500) * 0.01)
>>> result = garch_residual_forecast(returns, horizon=5)
>>> len(result['return_forecast'])
5
>>> len(result['vol_forecast'])
5

Notes

This approach is related to the “GARCH-in-mean” literature but separates the volatility and level dynamics for modularity.

See also

wraquant.vol.models.garch_fit: GARCH model fitting. wraquant.vol.models.garch_forecast: GARCH volatility forecasting. ses_forecast: SES for residual forecasting. theta_forecast: Theta method for residual forecasting.

Changepoint Detection

Change-point detection for time series.

cusum(data, threshold=1.0)

Cumulative sum (CUSUM) control chart for change detection.

Returns a CUSUM statistic series. Values exceeding threshold standard deviations indicate a potential shift.

Parameters:

• data (Series) – Time series.

• threshold (float, default: 1.0) – Detection threshold in standard deviation units.

Return type:

Series

Returns:

CUSUM statistic series.

detect_changepoints(data, method='pelt', penalty=None)

Detect change-points using the ruptures library.

Parameters:

• data (Series) – Time series.

• method (str, default: 'pelt') – Algorithm — "pelt" (default), "binseg", or "window".

• penalty (float | None, default: None) – Penalty value for the algorithm. When None, a sensible default is chosen.

Return type:

list[int]

Returns:

List of change-point indices (positions in the array, excluding the final point).

Raises:

ValueError – If method is not recognized.

Stationarity Tests

Stationarity tests and transformations for time series.

Provides both transformations (differencing, detrending) and formal statistical tests (ADF, KPSS, Phillips-Perron) for assessing and achieving stationarity in time series data.

difference(data, order=1)

Apply integer differencing to a time series.

Parameters:

• data (Series) – Time series.

• order (int, default: 1) – Differencing order (1 = first difference).

Return type:

Series

Returns:

Differenced series with NaN values dropped.

fractional_difference(data, d=0.5, threshold=1e-05)

Apply fractional differencing to preserve long-memory information.

Implements the fixed-width window fracdiff method from Advances in Financial Machine Learning (Lopez de Prado).

Parameters:

• data (Series) – Time series.

• d (float, default: 0.5) – Fractional differencing parameter (0 < d < 1).

• threshold (float, default: 1e-05) – Weight cutoff threshold for the window.

Return type:

Series

Returns:

Fractionally differenced series.

detrend(data, method='linear')

Remove trend from a time series.

Parameters:

• data (Series) – Time series.

• method (str, default: 'linear') – Detrending method — "linear" (default) or "constant" (demean).

Return type:

Series

Returns:

Detrended series.

Raises:

ValueError – If method is not recognized.

adf_test(data, max_lags=None, regression='c', significance=0.05)

Augmented Dickey-Fuller (ADF) unit root test.

Tests the null hypothesis that a unit root is present (i.e., the series is non-stationary). A low p-value leads to rejection of the null, concluding stationarity.

The test regression is:

Delta y_t = alpha + beta*t + gamma*y_{t-1} + sum_i delta_i * Delta y_{t-i} + e_t

where alpha is a constant (if regression='c'), beta*t is a time trend (if regression='ct'), and the number of lagged differences is chosen to remove serial correlation in the residuals.

ADF vs KPSS:

• ADF: null = unit root (non-stationary). Rejecting H0 supports stationarity.

• KPSS: null = stationary. Rejecting H0 supports non-stationarity.

• Best practice: run both. If ADF rejects and KPSS does not reject, strong evidence of stationarity. If both reject or both fail to reject, the series may be trend-stationary or difference-stationary.

Parameters:

• data (Series) – Time series to test. NaN values are dropped.

• max_lags (int | None, default: None) – Maximum number of lags for the lagged differences. If None, the optimal lag is selected automatically using the AIC criterion.

• regression (str, default: 'c') – Deterministic terms to include: "c" – constant only (default, most common). "ct" – constant and linear trend. "n" – no constant, no trend.

• significance (float, default: 0.05) – Significance level for the is_stationary convenience flag (default 0.05).

Returns:

• test_statistic: float, the ADF t-statistic.

• p_value: float, MacKinnon approximate p-value.

• optimal_lag: int, number of lags used.

• n_obs: int, number of observations used in the regression.

• critical_values: dict mapping significance levels ("1%", "5%", "10%") to critical values.

• is_stationary: bool, True if p-value < significance.

• interpretation: str, human-readable conclusion.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> # White noise is stationary
>>> stationary = pd.Series(rng.normal(0, 1, 500))
>>> result = adf_test(stationary)
>>> result['is_stationary']
True

References

• Dickey, D.A. & Fuller, W.A. (1979), “Distribution of the Estimators for Autoregressive Time Series With a Unit Root”, JASA.

• Said, S.E. & Dickey, D.A. (1984), “Testing for Unit Roots in Autoregressive-Moving Average Models of Unknown Order”, Biometrika.

kpss_test(data, regression='c', n_lags='auto', significance=0.05)

Kwiatkowski-Phillips-Schmidt-Shin (KPSS) stationarity test.

Tests the null hypothesis that the series is stationary around a deterministic trend. This is the opposite of the ADF test: here, rejecting H0 implies non-stationarity.

When to use ADF vs KPSS vs both:

• ADF alone: quick check; but has low power against near-unit-root alternatives.

• KPSS alone: confirms stationarity; but may over-reject for long-memory processes.

• Both (recommended): a confirmatory strategy. Four cases: 1. ADF rejects, KPSS does not reject -> stationary. 2. ADF does not reject, KPSS rejects -> non-stationary. 3. Both reject -> trend-stationary (difference to remove trend). 4. Neither rejects -> inconclusive, may need more data.

Parameters:

• data (Series) – Time series to test. NaN values are dropped.

• regression (str, default: 'c') – Deterministic component: "c" – test for level stationarity (default). "ct" – test for trend stationarity.

• n_lags (int | str, default: 'auto') – Number of lags for the Newey-West estimator. "auto" (default) uses the data-dependent rule. An integer value fixes the lag truncation.

• significance (float, default: 0.05) – Significance level for the is_stationary convenience flag (default 0.05).

Returns:

• test_statistic: float, the KPSS statistic.

• p_value: float, interpolated p-value (may be bounded at 0.01 or 0.10 depending on the table).

• n_lags: int, number of lags used.

• critical_values: dict mapping significance levels to critical values.

• is_stationary: bool, True if p-value >= significance (i.e., cannot reject the null of stationarity).

• interpretation: str, human-readable conclusion.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> stationary = pd.Series(rng.normal(0, 1, 500))
>>> result = kpss_test(stationary)
>>> result['is_stationary']
True

References

• Kwiatkowski, D. et al. (1992), “Testing the null hypothesis of stationarity against the alternative of a unit root”, Journal of Econometrics.

phillips_perron(data, regression='c', significance=0.05)

Phillips-Perron (PP) unit root test.

Like the ADF test, the PP test has the null hypothesis of a unit root (non-stationarity). However, the PP test uses a non-parametric correction to the t-statistic that is robust to heteroskedasticity and serial correlation in the error term, without needing to specify a lag order.

When to prefer PP over ADF:

• When the error process exhibits heteroskedasticity (e.g., financial returns with volatility clustering).

• When you are uncertain about the appropriate lag length for ADF.

• As a robustness check alongside ADF.

The PP test statistic modifies the ADF statistic with a correction factor based on the Newey-West long-run variance estimate.

Parameters:

• data (Series) – Time series to test. NaN values are dropped.

• regression (str, default: 'c') – Deterministic terms: "c" – constant only (default). "ct" – constant and trend. "n" – no constant, no trend.

• significance (float, default: 0.05) – Significance level for the is_stationary convenience flag (default 0.05).

Returns:

• test_statistic: float, the PP t-statistic.

• p_value: float, MacKinnon approximate p-value.

• n_lags: int, truncation lag for Newey-West.

• critical_values: dict mapping significance levels to critical values.

• is_stationary: bool, True if p-value < significance.

• interpretation: str, human-readable conclusion.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> stationary = pd.Series(rng.normal(0, 1, 500))
>>> result = phillips_perron(stationary)
>>> result['is_stationary']
True

References

• Phillips, P.C.B. & Perron, P. (1988), “Testing for a Unit Root in Time Series Regression”, Biometrika.

optimal_differencing(data, max_d=2, significance=0.05)

Automatically determine the optimal differencing order for stationarity.

Sequentially applies integer differencing (d = 0, 1, 2, …) and runs the ADF test at each order. Returns the smallest d for which the ADF test rejects the unit root null.

Parameters:

• data (Series) – Time series. NaN values are dropped.

• max_d (int, default: 2) – Maximum differencing order to try (default 2). Higher orders are rarely needed in practice.

• significance (float, default: 0.05) – Significance level for the ADF test (default 0.05).

Returns:

• optimal_d: int, the smallest differencing order that achieves stationarity (or max_d if stationarity is not achieved).

• test_results: dict mapping each d to its adf_test result dictionary.

• is_stationary: bool, True if stationarity was achieved within max_d differences.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> # Random walk needs d=1
>>> rw = pd.Series(np.cumsum(rng.normal(0, 1, 500)))
>>> result = optimal_differencing(rw)
>>> result['optimal_d']
1

References

• Hyndman, R.J. & Khandakar, Y. (2008), “Automatic Time Series Forecasting: The forecast Package for R”, JSS.

variance_ratio_test(data, lags=2, overlap=True, significance=0.05)

Lo-MacKinlay variance ratio test for the random walk hypothesis.

Tests whether the variance of k-period returns scales linearly with k, as implied by a random walk. A variance ratio VR(k) = 1 is consistent with a random walk; VR(k) > 1 suggests positive autocorrelation (momentum); VR(k) < 1 suggests mean reversion.

The test statistic under homoskedasticity is:

VR(k) = Var(r_t(k)) / (k * Var(r_t)) z = (VR(k) - 1) / sqrt(2*(2k-1)*(k-1) / (3*k*T))

A heteroskedasticity-robust version is also computed.

Parameters:

• data (Series) – Time series of prices or levels (NOT returns). NaN values are dropped.

• lags (int, default: 2) – Holding period / aggregation interval k (default 2). Common choices: 2, 4, 8, 16.

• overlap (bool, default: True) – Use overlapping returns for better power (default True). Non-overlapping uses fewer observations.

• significance (float, default: 0.05) – Significance level for the is_random_walk convenience flag (default 0.05).

Returns:

• variance_ratio: float, the VR(k) statistic.

• z_statistic: float, the z-statistic (homoskedastic).

• z_robust: float, heteroskedasticity-robust z-statistic.

• p_value: float, two-sided p-value from the robust z-statistic.

• is_random_walk: bool, True if the random walk null cannot be rejected at the given significance level.

• interpretation: str, human-readable conclusion.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> # Random walk: VR should be close to 1
>>> rw = pd.Series(np.cumsum(rng.normal(0, 1, 1000)))
>>> result = variance_ratio_test(rw, lags=2)
>>> 0.5 < result['variance_ratio'] < 1.5
True

References

• Lo, A.W. & MacKinlay, A.C. (1988), “Stock Market Prices Do Not Follow Random Walks: Evidence from a Simple Specification Test”, Review of Financial Studies.

Seasonality

Seasonality detection and feature engineering.

detect_seasonality(data, max_period=365)

Detect dominant seasonal periods via spectral analysis.

Uses Welch’s periodogram to identify significant frequency peaks.

Parameters:

• data (Series) – Time series to analyse.

• max_period (int, default: 365) – Maximum period to consider.

Return type:

list[int]

Returns:

List of detected seasonal periods (in number of observations), sorted by spectral power (strongest first).

fourier_features(index, period, n_harmonics)

Generate Fourier sine/cosine features for a datetime index.

Useful for encoding seasonality as regression features.

Parameters:

• index (DatetimeIndex) – Datetime index.

• period (int) – Seasonal period (in the same time unit as the index).

• n_harmonics (int) – Number of Fourier harmonics to generate.

Return type:

DataFrame

Returns:

DataFrame with sin_k and cos_k columns for each harmonic.

multi_fourier_features(index, periods, n_harmonics=3)

Generate Fourier terms for multiple seasonal periods.

Creates sin/cos pairs for each period-harmonic combination, producing a feature matrix suitable for regression-based seasonal modelling (e.g., linear regression, Prophet-style models, or as exogenous regressors for ARIMA).

For a period P and harmonic k, the features are:

sin(2 * pi * k * t / P) and cos(2 * pi * k * t / P)

This captures up to k-th order seasonal variation within each period.

Parameters:

• index (DatetimeIndex) – Datetime index of the time series.

• periods (list[int]) – List of seasonal periods (in the same units as the index). For example, [7, 365] for daily data with weekly and yearly seasonality.

• n_harmonics (int | list[int], default: 3) – Number of harmonics per period. Can be a single int (applied to all periods) or a list of ints (one per period). Default 3.

Return type:

DataFrame

Returns:

DataFrame with columns named sin_P{period}_H{harmonic} and cos_P{period}_H{harmonic} for each period-harmonic pair. Shape is (len(index), 2 * sum(harmonics)).

Example

>>> import pandas as pd
>>> idx = pd.date_range("2020-01-01", periods=365, freq="D")
>>> df = multi_fourier_features(idx, periods=[7, 365], n_harmonics=3)
>>> df.shape[1]
12

seasonal_strength(data, period=None)

Quantify the strength of seasonality in a time series.

Uses the Wang, Smith, and Hyndman (2006) measure:

F_s = max(0, 1 - Var(R_t) / Var(S_t + R_t))

where S_t is the seasonal component and R_t is the remainder from an STL decomposition.

Returns a float in [0, 1]:

• 0.0: no detectable seasonality (remainder dominates).

• 1.0: perfectly seasonal (no noise).

• > 0.9 is typically considered “strong” seasonality.

• < 0.4 is typically “weak”.

Parameters:

• data (Series) – Time series to evaluate. NaN values are dropped.

• period (int | None, default: None) – Seasonal period for the STL decomposition. If None, defaults to 7.

Return type:

float

Returns:

Strength of seasonality as a float in [0, 1].

Example

>>> import numpy as np, pandas as pd
>>> t = np.arange(200, dtype=float)
>>> pure_seasonal = pd.Series(10 * np.sin(2 * np.pi * t / 20))
>>> strength = seasonal_strength(pure_seasonal, period=20)
>>> strength > 0.9
True

References

• Wang, X., Smith, K. & Hyndman, R. (2006), “Characteristic- Based Clustering for Time Series Data”, Data Mining and Knowledge Discovery.

multi_seasonal_decompose(data, periods)

Decompose a time series with multiple seasonal periods.

Handles complex seasonality (e.g., daily data with both weekly and yearly patterns) by using MSTL (Multiple Seasonal-Trend decomposition using Loess) from statsmodels when available, or an iterative STL approach as fallback.

Parameters:

• data (Series) – Time series to decompose. NaN values are dropped.

• periods (list[int]) – List of seasonal periods, ordered from shortest to longest. For example, [7, 365] for weekly and yearly seasonality in daily data.

Returns:

• trend: pd.Series of the trend component.

• seasonal: dict mapping each period to its pd.Series seasonal component.

• residual: pd.Series of the residual after removing trend and all seasonal components.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> t = np.arange(730, dtype=float)
>>> weekly = 3 * np.sin(2 * np.pi * t / 7)
>>> yearly = 5 * np.sin(2 * np.pi * t / 365)
>>> trend = 0.01 * t
>>> data = pd.Series(trend + weekly + yearly)
>>> result = multi_seasonal_decompose(data, periods=[7, 365])
>>> sorted(result['seasonal'].keys())
[7, 365]

References

• Bandara, K. et al. (2021), “MSTL: A Seasonal-Trend Decomposition Algorithm for Time Series with Multiple Seasonal Patterns”, arXiv:2107.13462.

Feature Extraction

Time series feature extraction.

Provides functions to compute summary features from time series data for use in classification, clustering, or exploratory analysis. Features fall into three categories:

1. Autocorrelation features: ACF, PACF, Ljung-Box test for serial dependence structure.

2. Spectral features: power spectral density, dominant frequency, spectral entropy, spectral flatness.

3. Complexity features: sample entropy, approximate entropy, permutation entropy for measuring signal regularity/predictability.

autocorrelation_features(data, n_lags=40, significance=0.05)

Compute autocorrelation and partial autocorrelation features.

Extracts the ACF and PACF values up to n_lags, identifies statistically significant lags, and runs the Ljung-Box test for overall serial correlation.

Interpretation:

• ACF measures linear dependence between y_t and y_{t-k}. For AR(p) processes, ACF decays exponentially. For MA(q), ACF cuts off after lag q.

• PACF measures the correlation between y_t and y_{t-k} after removing the effect of intermediate lags. For AR(p), PACF cuts off after lag p.

• Ljung-Box: tests H0 that the first m autocorrelations are jointly zero. A significant result indicates the series has exploitable temporal structure.

Parameters:

• data (Series) – Time series. NaN values are dropped.

• n_lags (int, default: 40) – Number of lags to compute (default 40).

• significance (float, default: 0.05) – Significance level for identifying significant lags and Ljung-Box test (default 0.05).

Returns:

• acf: 1-D numpy array of ACF values (lags 0 to n_lags).

• pacf: 1-D numpy array of PACF values (lags 0 to n_lags).

• significant_acf_lags: list of lag indices where ACF exceeds the Bartlett confidence band.

• significant_pacf_lags: list of lag indices where PACF exceeds the confidence band.

• ljung_box: dict with statistic, p_value, and is_significant for the Ljung-Box test at the maximum lag.

• first_significant_lag: int or None, the first lag with significant ACF (excluding lag 0).

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> # AR(1) process
>>> x = np.zeros(500)
>>> for i in range(1, 500):
...     x[i] = 0.7 * x[i-1] + rng.normal()
>>> result = autocorrelation_features(pd.Series(x), n_lags=20)
>>> result['acf'][1] > 0.5  # strong lag-1 autocorrelation
True

spectral_features(data, fs=1.0)

Compute power spectral density features.

Extracts summary statistics from the periodogram including the dominant frequency, spectral entropy (a measure of spectral complexity), and spectral flatness (how noise-like the signal is).

Parameters:

• data (Series) – Time series. NaN values are dropped.

• fs (float, default: 1.0) – Sampling frequency (default 1.0). For daily data with annual cycles, use fs=1.0 and interpret frequencies as cycles per observation.

Returns:

• dominant_frequency: float, the frequency with the highest spectral power (excluding DC component).

• dominant_period: float, the period corresponding to the dominant frequency (1 / dominant_frequency).

• spectral_entropy: float, Shannon entropy of the normalised power spectral density. Higher values indicate a more complex, broadband signal; lower values indicate a few dominant frequencies. Range: [0, log2(N/2)].

• spectral_flatness: float in [0, 1]. Ratio of geometric mean to arithmetic mean of the PSD. A value near 1 indicates white noise; near 0 indicates a tonal/periodic signal. Also known as Wiener entropy.

• spectral_centroid: float, the weighted mean frequency (centre of mass of the spectrum).

• psd: 1-D numpy array of power spectral density values.

• frequencies: 1-D numpy array of corresponding frequencies.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> t = np.arange(500, dtype=float)
>>> pure_tone = pd.Series(np.sin(2 * np.pi * t / 20))  # period=20
>>> result = spectral_features(pure_tone)
>>> abs(result['dominant_period'] - 20) < 2
True
>>> result['spectral_flatness'] < 0.3  # tonal signal
True

complexity_features(data, m=2, r=None, order=3)

Compute entropy-based complexity measures for a time series.

Extracts three entropy measures that quantify the regularity and predictability of the signal:

1. Sample entropy (SampEn): probability that patterns similar for m points remain similar for m+1 points. Lower values indicate more regularity (self-similarity); higher values indicate more complexity. Defined as -ln(A/B) where A is the count of m+1 length template matches and B is the count of m length matches, within tolerance r.

2. Approximate entropy (ApEn): similar to SampEn but includes self-matches, making it slightly biased toward regularity. Historically introduced first (Pincus 1991).

3. Permutation entropy (PeEn): based on the frequency of ordinal patterns of length order. Robust to noise and monotonic transformations. Normalised to [0, 1].

Parameters:

• data (Series) – Time series. NaN values are dropped.

• m (int, default: 2) – Embedding dimension for SampEn and ApEn (default 2). Typical values: 2 or 3.

• r (float | None, default: None) – Tolerance threshold for SampEn and ApEn. If None, defaults to 0.2 * std(data) (standard choice).

• order (int, default: 3) – Ordinal pattern length for permutation entropy (default 3). Typical values: 3-7.

Returns:

• sample_entropy: float, SampEn(m, r). Values typically range from 0 (perfectly regular) to ~2.5 (random).

• approximate_entropy: float, ApEn(m, r).

• permutation_entropy: float, normalised PeEn in [0, 1]. 0 = perfectly predictable, 1 = maximally complex/random.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> regular = pd.Series(np.sin(np.arange(500) * 0.1))
>>> random = pd.Series(rng.normal(0, 1, 500))
>>> reg_feat = complexity_features(regular)
>>> rnd_feat = complexity_features(random)
>>> reg_feat['sample_entropy'] < rnd_feat['sample_entropy']
True

References

• Richman, J.S. & Moorman, J.R. (2000), “Physiological time-series analysis using approximate entropy and sample entropy”, AJP.

• Pincus, S.M. (1991), “Approximate entropy as a measure of system complexity”, PNAS.

• Bandt, C. & Pompe, B. (2002), “Permutation Entropy: A Natural Complexity Measure for Time Series”, PRL.

Anomaly Detection

Time series anomaly detection.

Provides methods for detecting anomalous observations in time series data using statistical and machine learning approaches:

1. Isolation Forest on rolling features – detects contextual anomalies by featurising the local window (mean, std, skew) and running an isolation forest on the feature space.

2. Forecast-based anomaly detection – fits a time series model and flags observations that deviate from the forecast by more than k standard deviations. Does NOT require fbprophet.

3. Rolling Grubbs test – applies the classical Grubbs outlier test in a rolling window to detect local outliers.

isolation_forest_ts(data, window=20, contamination=0.05, features=None, random_state=42)

Anomaly detection using Isolation Forest on rolling features.

Rather than applying Isolation Forest directly to the raw values, this function first computes rolling window features (mean, std, skew, min, max) to capture the context of each observation. This approach detects contextual anomalies – values that are unusual given their local context, even if they are not globally extreme.

Parameters:

• data (Series) – Time series. NaN values are dropped.

• window (int, default: 20) – Rolling window size for feature computation (default 20).

• contamination (float, default: 0.05) – Expected proportion of anomalies in the data (default 0.05). Controls the threshold for the Isolation Forest decision function.

• features (list[str] | None, default: None) – List of rolling features to compute. If None, uses ["mean", "std", "skew", "min", "max"]. Supported values: "mean", "std", "skew", "min", "max", "median", "range".

• random_state (int, default: 42) – Random seed for reproducibility (default 42).

Returns:

• anomaly_scores: pd.Series of anomaly scores from the Isolation Forest decision function (lower = more anomalous).

• anomaly_labels: pd.Series of int labels (1 = normal, -1 = anomaly).

• threshold: float, the decision function threshold.

• anomaly_indices: list of index values flagged as anomalies.

• n_anomalies: int, number of detected anomalies.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 500)
>>> x[100] = 20.0  # inject anomaly
>>> x[300] = -15.0  # inject anomaly
>>> data = pd.Series(x)
>>> result = isolation_forest_ts(data, window=20)
>>> result['n_anomalies'] > 0
True

References

• Liu, F.T. et al. (2008), “Isolation Forest”, ICDM.

prophet_anomaly(data, k=3.0, seasonal_period=None, trend=True)

Forecast-based anomaly detection using statsmodels.

Fits an Unobserved Components Model (local level + optional seasonal) and flags observations where the actual value deviates from the smoothed state by more than k standard deviations of the residuals.

This approach does NOT require fbprophet or Prophet. It uses pure statsmodels for a fully open-source, lightweight solution.

Parameters:

• data (Series) – Time series. NaN values are dropped.

• k (float, default: 3.0) – Number of standard deviations for the anomaly threshold (default 3.0). Lower values flag more anomalies.

• seasonal_period (int | None, default: None) – Seasonal period for the model. If None, no seasonal component is included.

• trend (bool, default: True) – Include a trend component (default True).

Returns:

• forecast: pd.Series of smoothed / fitted values.

• residuals: pd.Series of (actual - forecast).

• anomaly_mask: pd.Series of bool, True for anomalies.

• anomaly_indices: list of index values flagged as anomalies.

• upper_bound: pd.Series, forecast + k * sigma.

• lower_bound: pd.Series, forecast - k * sigma.

• sigma: float, standard deviation of residuals.

• n_anomalies: int, number of detected anomalies.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 300)
>>> x[50] = 15.0  # anomaly
>>> x[200] = -12.0  # anomaly
>>> data = pd.Series(x)
>>> result = prophet_anomaly(data, k=3.0)
>>> result['n_anomalies'] >= 2
True

References

• Harvey, A.C. (1989), Forecasting, Structural Time Series Models and the Kalman Filter. Cambridge University Press.

grubbs_test_ts(data, window=50, significance=0.05)

Rolling Grubbs test for outlier detection in time series.

Applies the Grubbs test (maximum normed residual test) within a rolling window to detect observations that are statistically unlikely given the local distribution.

The Grubbs test statistic for a sample of size n is:

G = max|x_i - mean| / std

and the critical value is derived from the t-distribution:

G_crit = ((n-1) / sqrt(n)) * sqrt(t^2 / (n - 2 + t^2))

where t is the critical value of the t-distribution with n-2 degrees of freedom at significance alpha / (2n).

Parameters:

• data (Series) – Time series. NaN values are dropped.

• window (int, default: 50) – Rolling window size (default 50). Should be large enough for the Grubbs test to be meaningful (>= 7).

• significance (float, default: 0.05) – Significance level for the Grubbs test (default 0.05).

Returns:

• outlier_mask: pd.Series of bool, True for detected outliers.

• test_statistics: pd.Series of Grubbs G statistics at each position.

• outlier_indices: list of index values flagged as outliers.

• n_outliers: int, number of detected outliers.

Return type:

dict

Example

>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 200)
>>> x[75] = 20.0  # obvious outlier
>>> data = pd.Series(x)
>>> result = grubbs_test_ts(data, window=50)
>>> 75 in result['outlier_indices']
True

References

• Grubbs, F.E. (1950), “Sample Criteria for Testing Outlying Observations”, Annals of Mathematical Statistics.

Stochastic Processes

Stochastic process forecasting for financial time series.

Provides parameter estimation and forecasting for continuous-time stochastic processes commonly used in quantitative finance:

• Ornstein-Uhlenbeck (OU): mean-reverting diffusion for spreads, rates, and volatility.

• Merton Jump-Diffusion: GBM with Poisson jumps for fat-tailed asset returns.

• Regime-Switching: separate models per hidden regime, blended by regime probabilities.

• Vector Autoregression (VAR): multi-asset linear forecasting with impulse response and variance decomposition.

ornstein_uhlenbeck_forecast(data, h=10, dt=None, n_simulations=1000, confidence_level=0.95)

Forecast using an Ornstein-Uhlenbeck (OU) mean-reverting process.

The OU process is the continuous-time analogue of an AR(1) process and is the workhorse model for mean-reverting financial quantities.

Use for mean-reverting assets like spreads, interest rates, volatility (VIX), pairs-trading residuals, or any series that fluctuates around a long-run equilibrium.

SDE:

dX_t = theta * (mu - X_t) * dt + sigma * dW_t

where:

• theta > 0: speed of mean reversion (higher = faster).

• mu: long-run mean level.

• sigma: volatility of the diffusion.

• W_t: standard Brownian motion.

The half-life of mean reversion is ln(2) / theta, giving the expected time for a deviation to halve.

Parameter estimation:

Parameters are estimated via Maximum Likelihood on the discrete-time transition density, which is Gaussian: X_{t+dt} | X_t ~ N(mu + (X_t - mu) * exp(-theta*dt), sigma^2 / (2*theta) * (1 - exp(-2*theta*dt)))

Parameters:

• data (Series) – Time series of observed values.

• h (int, default: 10) – Number of steps to forecast (default 10).

• dt (float | None, default: None) – Time step between observations. If None, inferred from the index (1/252 for business days, 1.0 otherwise).

• n_simulations (int, default: 1000) – Number of Monte Carlo paths for confidence bands (default 1000).

• confidence_level (float, default: 0.95) – Confidence level for intervals (default 0.95).

Returns:

• params: dict with theta, mu, sigma.

• half_life: half-life of mean reversion (in observation units).

• forecast: pd.Series of expected path (conditional mean).

• confidence_intervals: dict with lower and upper pd.Series.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> # Simulate OU process
>>> n, theta, mu, sigma = 500, 5.0, 100.0, 2.0
>>> x = np.zeros(n); x[0] = 100.0
>>> dt = 1 / 252
>>> for i in range(1, n):
...     x[i] = x[i-1] + theta * (mu - x[i-1]) * dt + sigma * np.sqrt(dt) * rng.normal()
>>> data = pd.Series(x)
>>> result = ornstein_uhlenbeck_forecast(data, h=20)
>>> result['half_life'] > 0
True

References

Uhlenbeck, G.E. & Ornstein, L.S. (1930). On the theory of Brownian motion. Physical Review, 36(5), 823.

jump_diffusion_forecast(data, h=10, dt=None, n_paths=1000, confidence_level=0.95, seed=None)

Merton jump-diffusion model forecast via Monte Carlo simulation.

Extends geometric Brownian motion (GBM) by adding Poisson-driven jumps, capturing the fat tails and sudden moves observed in real asset returns.

SDE:

dS/S = (mu - lambda*k) * dt + sigma * dW + J * dN

where:

• mu: drift of the diffusion component.

• sigma: diffusion volatility.

• lambda: jump intensity (expected jumps per unit time).

• J: jump size, ln(1+J) ~ N(mu_j, sigma_j^2).

• N_t: Poisson process with intensity lambda.

• k = exp(mu_j + sigma_j^2/2) - 1.

Parameter estimation:

A moment-matching approach is used on log-returns: - Jump intensity estimated from excess kurtosis. - Jump mean and variance from the non-Gaussian tail structure.

Parameters:

• data (Series) – Price series (positive values).

• h (int, default: 10) – Forecast horizon in steps (default 10).

• dt (float | None, default: None) – Time step between observations (default 1/252 for daily).

• n_paths (int, default: 1000) – Number of Monte Carlo simulation paths (default 1000).

• confidence_level (float, default: 0.95) – Confidence level for intervals (default 0.95).

• seed (int | None, default: None) – Random seed for reproducibility.

Returns:

• params: dict with mu, sigma, lambda_, mu_j, sigma_j.

• forecast_paths: np.ndarray of shape (n_paths, h) with simulated price paths.

• mean_forecast: pd.Series of mean across paths.

• confidence_intervals: dict with lower and upper pd.Series.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> prices = pd.Series(100 * np.exp(np.cumsum(np.random.randn(200) * 0.01)))
>>> result = jump_diffusion_forecast(prices, h=10, n_paths=500, seed=42)
>>> result['forecast_paths'].shape
(500, 10)

References

Merton, R.C. (1976). Option pricing when underlying stock returns are discontinuous. Journal of Financial Economics, 3(1-2), 125-144.

regime_switching_forecast(data, n_regimes=2, h=10)

Forecast using a Markov regime-switching model.

Fits a Markov-switching autoregression where the series parameters (mean, variance) switch between discrete hidden regimes. Forecasts are blended across regimes weighted by the predicted regime probabilities.

When to use:

• Series that alternate between distinct states (e.g., bull/bear markets, high/low volatility regimes).

• When a single-model forecast is inadequate because the data generation process changes over time.

Parameters:

• data (Series) – Time series to forecast.

• n_regimes (int, default: 2) – Number of hidden regimes (default 2).

• h (int, default: 10) – Forecast horizon (default 10).

Returns:

• forecast: pd.Series of blended point forecasts.

• regime_forecasts: dict mapping regime index to its pd.Series forecast.

• regime_probs: np.ndarray of shape (h, n_regimes) with predicted regime probabilities.

• blended: same as forecast (for clarity).

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> # Bull/bear regime data
>>> regimes = np.concatenate([rng.normal(0.1, 0.5, 150),
...                           rng.normal(-0.05, 1.0, 150)])
>>> data = pd.Series(np.cumsum(regimes))
>>> result = regime_switching_forecast(data, n_regimes=2, h=10)
>>> result['regime_probs'].shape[1]
2

References

Hamilton, J.D. (1989). A New Approach to the Economic Analysis of Nonstationary Time Series and the Business Cycle. Econometrica, 57(2), 357-384.

var_forecast(data, h=10, maxlags=None, ic='aic')

Vector Autoregression (VAR) for multi-asset forecasting.

VAR models each variable as a linear function of its own lags and the lags of all other variables. This captures lead-lag relationships (e.g., one asset predicting another) and provides impulse response functions (IRF) and forecast error variance decomposition (FEVD).

When to use:

• Forecasting multiple related time series simultaneously (e.g., a portfolio of assets, yield curve factors).

• Analysing dynamic interactions: which variables Granger-cause which.

• Computing impulse responses: how a shock to one variable propagates to others.

Math:

Y_t = c + A_1 * Y_{t-1} + ... + A_p * Y_{t-p} + e_t

where Y_t is a (k x 1) vector, A_i are (k x k) coefficient matrices, and e_t ~ N(0, Sigma).

Parameters:

• data (DataFrame) – DataFrame with multiple columns (one per variable). Should be stationary (difference if needed).

• h (int, default: 10) – Forecast horizon (default 10).

• maxlags (int | None, default: None) – Maximum lag order to consider. If None, uses 12 * (nobs/100)^(1/4) rule of thumb.

• ic (str, default: 'aic') – Information criterion for lag selection: "aic", "bic", "hqic", "fpe" (default "aic").

Returns:

• forecast: pd.DataFrame of point forecasts.

• irf: impulse response function results (dict of DataFrames keyed by shocked variable).

• fevd: forecast error variance decomposition (dict of DataFrames keyed by target variable).

• granger_causality: dict of p-values for Granger causality tests.

• lag_order: selected lag order.

Return type:

dict

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> n = 200
>>> x1 = np.cumsum(rng.normal(0, 1, n))
>>> x2 = 0.3 * np.roll(x1, 1) + rng.normal(0, 1, n)
>>> data = pd.DataFrame({'x1': np.diff(x1), 'x2': np.diff(x2)})
>>> result = var_forecast(data, h=5)
>>> result['forecast'].shape[0]
5

References

Lutkepohl, H. (2005). New Introduction to Multiple Time Series Analysis. Springer.

Advanced Methods

Advanced time series integrations using optional packages.

Provides wrappers around tsfresh, stumpy, pywavelets, sktime, statsforecast, and tslearn for feature extraction, matrix profiles, wavelet analysis, forecasting, and time series clustering.

tsfresh_features(df, column_id='id', column_sort='time')

Extract time series features using tsfresh.

Wraps tsfresh.extract_features with efficient defaults and returns a cleaned DataFrame with no NaN/infinite columns.

Parameters:

• df (DataFrame) – Long-format DataFrame containing time series data. Must include the columns specified by column_id and column_sort, plus one or more value columns.

• column_id (str, default: 'id') – Column identifying distinct time series.

• column_sort (str, default: 'time') – Column used for temporal ordering.

Returns:

Extracted features indexed by column_id values. Columns with NaN or infinite values are dropped.

Return type:

DataFrame

stumpy_matrix_profile(ts, m=50)

Compute the matrix profile of a time series using STUMPY.

The matrix profile identifies the nearest-neighbour distance for every subsequence of length m, enabling motif discovery and discord (anomaly) detection.

Parameters:

• ts (ndarray | Series) – Univariate time series.

• m (int, default: 50) – Subsequence window length.

Returns:

Dictionary containing:

• matrix_profile – 1-D array of nearest-neighbour distances.

• profile_index – 1-D array of indices of the nearest neighbour.

• motif_idx – index of the top motif (lowest MP value).

• discord_idx – index of the top discord (highest MP value).

Return type:

dict[str, Any]

wavelet_transform(data, wavelet='db4', level=None)

Perform discrete wavelet transform (DWT) using PyWavelets.

Parameters:

• data (ndarray | Series) – Univariate signal to decompose.

• wavelet (str, default: 'db4') – Wavelet name (e.g. 'db4', 'haar', 'sym5').

• level (int | None, default: None) – Decomposition level. When None, the maximum useful level is computed automatically.

Returns:

Dictionary containing:

• coeffs – list of coefficient arrays [cA_n, cD_n, ..., cD_1].

• wavelet – wavelet name used.

• level – decomposition level used.

Return type:

dict[str, Any]

wavelet_denoise(data, wavelet='db4', level=None, threshold='soft')

Denoise a signal using wavelet thresholding.

Applies universal (VisuShrink) thresholding to the detail coefficients after a DWT decomposition.

Parameters:

• data (ndarray | Series) – Noisy univariate signal.

• wavelet (str, default: 'db4') – Wavelet name.

• level (int | None, default: None) – Decomposition level (auto-selected when None).

• threshold (str, default: 'soft') – Thresholding mode: 'soft' or 'hard'.

Returns:

Denoised signal with the same length as the input.

Return type:

ndarray

sktime_forecast(y, model='naive', horizon=10)

Forecast a time series using sktime’s unified interface.

Parameters:

• y (Series) – Historical observations indexed by a pandas PeriodIndex or integer index.

• model (str, default: 'naive') –

  Forecasting model. Supported values:

  • 'naive' – last-value forecast

  • 'theta' – Theta method

  • 'ets' – exponential smoothing (AutoETS)

• horizon (int, default: 10) – Number of steps ahead to forecast.

Returns:

DataFrame with columns forecast and index covering the forecast horizon.

Return type:

DataFrame

statsforecast_auto(y, season_length=1, horizon=10)

Automatic forecasting using Nixtla’s StatsForecast.

Runs AutoARIMA and returns point forecasts.

Parameters:

• y (Series) – Historical observations. Must have a DatetimeIndex or integer index.

• season_length (int, default: 1) – Seasonal period length (e.g. 12 for monthly data with yearly seasonality).

• horizon (int, default: 10) – Number of periods to forecast.

Returns:

DataFrame with column forecast indexed over the forecast horizon.

Return type:

DataFrame

tslearn_dtw(ts1, ts2)

Compute DTW distance between two time series using tslearn.

Parameters:

• ts1 (ndarray | Series) – First time series.

• ts2 (ndarray | Series) – Second time series.

Returns:

Dictionary containing:

• distance – DTW distance.

• path – optimal alignment path as a list of (i, j) tuples.

Return type:

dict[str, Any]

tslearn_kmeans(dataset, n_clusters=3)

Cluster time series using DTW-based K-means from tslearn.

Parameters:

• dataset (ndarray | list[ndarray]) – Collection of time series. If a 2-D array, each row is treated as a separate time series. If a 3-D array, shape should be (n_series, n_timestamps, n_features).

• n_clusters (int, default: 3) – Number of clusters.

Returns:

Dictionary containing:

• labels – cluster assignments (1-D int array).

• inertia – sum of distances to nearest cluster centre.

• cluster_centers – array of cluster centre time series.

• n_clusters – number of clusters used.

Return type:

dict[str, Any]

darts_forecast(data, model='nbeats', horizon=10, **kwargs)

Forecast using Darts deep learning models.

Darts provides production-quality implementations of N-BEATS, N-HiTS, Temporal Fusion Transformer, and other deep forecasting models. Use when you need state-of-the-art accuracy and have sufficient data (>1000 observations).

The model is trained on the provided data and produces a forecast of the specified horizon. Default training parameters are tuned for quick experimentation; pass additional keyword arguments to override them (e.g., n_epochs=200).

Parameters:

• data (ndarray | Series) – Price or return series. Must be univariate.

• model (str, default: 'nbeats') –

  Deep learning model to use:

  • 'nbeats' – N-BEATS (Neural Basis Expansion Analysis).

  • 'nhits' – N-HiTS (Neural Hierarchical Interpolation).

  • 'tcn' – Temporal Convolutional Network.

  • 'rnn' – Simple RNN (LSTM-based).

  • 'transformer' – Vanilla Transformer.

• horizon (int, default: 10) – Number of steps ahead to forecast.

• **kwargs (Any) – Additional keyword arguments passed to the model constructor (e.g., n_epochs, input_chunk_length).

Returns:

Dictionary containing:

• forecast – 1-D numpy array of forecasted values.

• model_name – name of the model used.

• training_loss – final training loss (when available).

Return type:

dict[str, Any]

Example

>>> import numpy as np
>>> from wraquant.ts.advanced import darts_forecast
>>> data = np.cumsum(np.random.default_rng(0).normal(0, 1, 500))
>>> result = darts_forecast(data, model="nbeats", horizon=5, n_epochs=5)
>>> len(result["forecast"])
5

Notes

Reference: Oreshkin et al. (2020). “N-BEATS: Neural basis expansion analysis for interpretable time series forecasting.” ICLR 2020.

See also

sktime_forecast

Simpler statistical forecasting models.

statsforecast_auto

AutoARIMA and friends.

auto_forecast

wraquant’s built-in automatic forecasting.