Statistics ( wraquant.stats )<a class="headerlink" href="#statistics-wraquant-stats" title="Link to this heading">¶

Returns:

Annualized return as a float.

annualized_volatility(returns, periods_per_year=252)[source]¶

Compute annualized volatility from a simple return series.

Parameters:

returns (Series) – Simple return series.
periods_per_year (int, default: 252) – Number of periods per year (252 for daily).

Return type:

Returns:

Annualized volatility as a float.

max_drawdown(prices)[source]¶

Compute maximum drawdown from a price series.

Parameters:: prices (Series) – Price series (not returns).
Return type:: float
Returns:: Maximum drawdown as a negative float (e.g., -0.25 for 25% drawdown).

calmar_ratio(returns, periods_per_year=252)[source]¶

Compute the Calmar ratio (annualized return / max drawdown).

Parameters:

returns (Series) – Simple return series.
periods_per_year (int, default: 252) – Number of periods per year.

Return type:

Returns:

Calmar ratio as a float.

omega_ratio(returns, threshold=0.0)[source]¶

Compute the Omega ratio.

The Omega ratio is the probability-weighted ratio of gains versus losses relative to a threshold.

Parameters:

returns (Series) – Simple return series.
threshold (float, default: 0.0) – Return threshold (default 0).

Return type:

Returns:

Omega ratio as a float.

rolling_sharpe(returns, window=60, risk_free_rate=0.0, periods_per_year=252)[source]¶

Compute the rolling Sharpe ratio over a moving window.

The Sharpe ratio is the most widely used risk-adjusted performance measure. The rolling variant shows how risk-adjusted performance evolves over time, revealing periods of strong and weak risk-adjusted returns.

When to use:

To monitor strategy performance stability over time.
To detect regime changes in risk-adjusted returns (e.g., a strategy that worked pre-2020 but degraded post-2020).
To compare two strategies’ time-varying risk-adjusted performance.

Mathematical formulation:

For each window of length w:

\[\text{Sharpe}_t = \frac{\bar{r}_t - r_f}{\sigma_t} \cdot \sqrt{P}\]

where \bar{r}_t and \sigma_t are the rolling mean and standard deviation of returns, r_f is the per-period risk-free rate, and P is the annualisation factor (e.g., 252 for daily).

How to interpret:

Sharpe > 1.0: good risk-adjusted performance (annualised).
Sharpe > 2.0: very strong.
Sharpe < 0.0: losing money on a risk-adjusted basis.
Large swings in rolling Sharpe indicate unstable performance.

Parameters:

returns (Series) – Simple return series.
window (int, default: 60) – Rolling window size in periods (default 60, roughly 3 months of daily data).
risk_free_rate (float, default: 0.0) – Per-period risk-free rate (default 0.0).
periods_per_year (int, default: 252) – Annualisation factor (252 for daily data).

Return type:

Returns:

Rolling Sharpe ratio as a pd.Series. First window - 1 values are NaN.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> ret = pd.Series(np.random.normal(0.001, 0.02, 252))
>>> rs = rolling_sharpe(ret, window=60)
>>> rs.dropna().shape[0]
193

See also

annualized_volatility: Annualised standard deviation. calmar_ratio: Drawdown-based risk-adjusted return.

rolling_drawdown(returns, window=60)[source]¶

Compute the rolling maximum drawdown over a moving window.

For each time step, the maximum drawdown is computed using only the most recent window observations. This provides a time-varying measure of downside risk.

When to use:

To monitor the worst-case loss over a recent period.
To detect periods of elevated tail risk that may not show up in rolling volatility.
As an input to risk overlays that tighten exposure when recent drawdowns are deep.
To compare the downside risk profile of different strategies over time.

Mathematical formulation:

For each window [t - w + 1, t], compute the cumulative return series from the window’s returns, find the peak, and measure the maximum drop from peak to trough:

\[\text{MDD}_t = \min_{s \in [t-w+1, t]} \frac{P_s - \max_{u \le s} P_u}{\max_{u \le s} P_u}\]

How to interpret:

Values are negative (or zero). More negative = deeper drawdown.
A rolling drawdown of -0.10 means the portfolio lost 10% from its peak within the window.
Compare to static max drawdown to see if the worst period is concentrated or distributed.

Parameters:

returns (Series) – Simple return series.
window (int, default: 60) – Rolling window size in periods (default 60).

Return type:

Returns:

Rolling maximum drawdown as a pd.Series. First window - 1 values are NaN.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> ret = pd.Series(np.random.normal(0, 0.02, 252))
>>> rd = rolling_drawdown(ret, window=60)
>>> (rd.dropna() <= 0).all()
True

See also

max_drawdown: Full-sample maximum drawdown. calmar_ratio: Return / drawdown ratio.

return_attribution(portfolio_weights, benchmark_weights, portfolio_returns, benchmark_returns)[source]¶

Decompose portfolio excess return using the Brinson-Fachler model.

The Brinson model is the industry standard for performance attribution, decomposing the active return (portfolio minus benchmark) into three components: asset allocation, security selection, and interaction effects.

When to use:

To explain why a portfolio outperformed or underperformed its benchmark.
To separate the contribution of top-down asset allocation decisions from bottom-up security selection.
For reporting to investors or risk committees.

Mathematical formulation:

For each asset/sector i:

Allocation effect: (w^P_i - w^B_i) * (r^B_i - r^B_total)
Selection effect: w^B_i * (r^P_i - r^B_i)
Interaction effect: (w^P_i - w^B_i) * (r^P_i - r^B_i)
Total active return: sum(allocation + selection + interaction) = r^P_total - r^B_total

How to interpret:

Positive allocation: the portfolio overweighted sectors that outperformed the benchmark.
Positive selection: within each sector, the portfolio held better-performing securities.
The interaction term captures the joint effect.
The three components sum to the total excess return.

Parameters:

portfolio_weights (Series) – Portfolio weights per asset/sector.
benchmark_weights (Series) – Benchmark weights per asset/sector (same index).
portfolio_returns (Series) – Portfolio returns per asset/sector.
benchmark_returns (Series) – Benchmark returns per asset/sector (same index).

Returns:

allocation: total allocation effect (float).
selection: total selection effect (float).
interaction: total interaction effect (float).
total_excess: total excess return (float).
detail: DataFrame with per-asset breakdown.

Return type:

Example

>>> import pandas as pd
>>> pw = pd.Series({"Tech": 0.4, "Fin": 0.3, "Health": 0.3})
>>> bw = pd.Series({"Tech": 0.3, "Fin": 0.4, "Health": 0.3})
>>> pr = pd.Series({"Tech": 0.05, "Fin": 0.02, "Health": 0.03})
>>> br = pd.Series({"Tech": 0.04, "Fin": 0.03, "Health": 0.03})
>>> result = return_attribution(pw, bw, pr, br)
>>> abs(result["total_excess"] - (pw @ pr - bw @ br)) < 1e-10
True

See also

risk_contribution: Per-asset risk decomposition.

risk_contribution(weights, cov_matrix)[source]¶

Compute marginal risk contributions per asset.

Risk contribution measures how much each asset contributes to the total portfolio risk (standard deviation). This is the foundation of risk-parity and risk-budgeting portfolio construction.

When to use:

To understand where portfolio risk comes from.
To build risk-parity portfolios where each asset contributes equally to total risk.
For risk monitoring: detect when a single position dominates portfolio risk.
To compare intended risk budgets with realised risk allocation.

Mathematical formulation:

The marginal contribution to risk (MCR) for asset i is:

\[\text{MCR}_i = w_i \cdot \frac{(\Sigma w)_i}{\sigma_p}\]

where \Sigma is the covariance matrix, w is the weight vector, and \sigma_p = \sqrt{w' \Sigma w} is the portfolio standard deviation.

The marginal contributions sum to the total portfolio risk:

\[\sum_i \text{MCR}_i = \sigma_p\]

How to interpret:

Values are in the same units as portfolio standard deviation.
Each value represents the portion of total portfolio risk attributable to that asset.
Negative contributions are possible for assets that hedge overall portfolio risk.

Parameters:

weights (Series | ndarray) – Portfolio weights (1-D array or Series).
cov_matrix (DataFrame | ndarray) – Covariance matrix (2-D array or DataFrame).

Return type:

Returns:

pd.Series of marginal risk contributions per asset.

Example

>>> import pandas as pd, numpy as np
>>> w = pd.Series({"A": 0.5, "B": 0.3, "C": 0.2})
>>> cov = pd.DataFrame(
...     np.diag([0.04, 0.09, 0.01]),
...     index=["A", "B", "C"], columns=["A", "B", "C"],
... )
>>> rc = risk_contribution(w, cov)
>>> abs(rc.sum() - np.sqrt(w.values @ cov.values @ w.values)) < 1e-10
True

See also

return_attribution: Return decomposition (Brinson model). shrunk_covariance: Better covariance input for risk contributions.

test_normality(data, method='jarque_bera')[source]¶

Test whether a series is normally distributed.

Parameters:

data (Series) – Data series to test.
method (str, default: 'jarque_bera') – Test method — "jarque_bera" (default), "shapiro", or "dagostino".

Return type:

Returns:

Dictionary with statistic, p_value, and is_normal (at 5% significance level).

Raises:

ValueError – If method is not recognized.

test_stationarity(data, method='adf')[source]¶

Test whether a time series is stationary.

Parameters:

data (Series) – Time series to test.
method (str, default: 'adf') – Test method — "adf" (Augmented Dickey-Fuller, default) or "kpss".

Return type:

Returns:

Dictionary with statistic, p_value, and is_stationary (at 5% significance level).

Raises:

ValueError – If method is not recognized.

test_autocorrelation(data, nlags=10)[source]¶

Ljung-Box test for autocorrelation.

Parameters:

data (Series) – Time series to test.
nlags (int, default: 10) – Number of lags to test.

Return type:

Returns:

Dictionary with statistic (at max lag), p_value, is_autocorrelated (at 5% significance), and the full results DataFrame.

shapiro_wilk(data)[source]¶

Shapiro-Wilk test for normality.

The Shapiro-Wilk test is widely regarded as the most powerful normality test for small to moderate sample sizes (n < 5000). It is more sensitive than the Jarque-Bera test, which relies only on skewness and kurtosis, because it considers the full empirical distribution.

When to use:

When you have fewer than 2000 observations and need a reliable normality assessment (e.g., validating assumptions before parametric VaR, calibrating option pricing models).
As a complement to Jarque-Bera: Shapiro-Wilk catches departures in the center of the distribution that JB (which focuses on moments 3 and 4) may miss.
For validating regression residuals before using t-based confidence intervals.

Mathematical formulation:

\[W = \frac{\left(\sum_{i=1}^n a_i x_{(i)}\right)^2}{\sum_{i=1}^n (x_i - \bar{x})^2}\]

where x_{(i)} are the order statistics and a_i are tabulated constants derived from the expected values and covariance matrix of order statistics from a normal distribution.

How to interpret:

W is in (0, 1]. Values near 1 indicate normality.
Reject normality if p_value < 0.05.
For financial returns, rejection is typical (fat tails), confirming that Gaussian-based risk measures are unreliable.

Parameters:

data (Series | ndarray) – Data series or 1-D array. Sample size should be between 3 and 5000 (scipy limitation).

Returns:

statistic: Shapiro-Wilk W statistic.
p_value: p-value for H0: data is normally distributed.
is_normal: bool, True if p > 0.05.

Return type:

Example

>>> import numpy as np
>>> data = np.random.default_rng(42).normal(0, 1, 200)
>>> result = shapiro_wilk(data)
>>> result["is_normal"]
True

durbin_watson(residuals)[source]¶

Durbin-Watson test for first-order autocorrelation in residuals.

The Durbin-Watson statistic tests whether the residuals of a regression model exhibit first-order serial correlation. This is critical in financial econometrics where autocorrelated residuals invalidate standard OLS inference.

When to use:

After fitting any OLS regression to time-series data (e.g., CAPM beta estimation, factor models). Autocorrelated residuals mean standard errors are biased and t-statistics are unreliable.
As a diagnostic before deciding whether to use Newey-West (HAC) standard errors.
For model validation: significant autocorrelation suggests a missing variable or incorrect functional form.

Mathematical formulation:

\[DW = \frac{\sum_{t=2}^T (e_t - e_{t-1})^2}{\sum_{t=1}^T e_t^2}\]

How to interpret:

DW ≈ 2.0: no autocorrelation.
DW < 2.0: positive autocorrelation (residuals tend to have the same sign as their predecessor).
DW > 2.0: negative autocorrelation.
Rule of thumb: DW < 1.5 or DW > 2.5 indicates significant autocorrelation. For precise inference, compare to the Durbin-Watson tables for dL and dU critical values.

Parameters:

residuals (Series | ndarray) – Regression residuals (1-D array or Series).

Returns:

statistic: Durbin-Watson statistic (range [0, 4]).
interpretation: string describing the result.

Return type:

Example

>>> import numpy as np
>>> residuals = np.random.default_rng(42).normal(0, 1, 100)
>>> result = durbin_watson(residuals)
>>> 1.5 < result["statistic"] < 2.5  # no autocorrelation
True

See also

test_autocorrelation: Ljung-Box test for higher-order autocorrelation.

breusch_pagan(residuals, exog)[source]¶

Breusch-Pagan Lagrange Multiplier test for heteroskedasticity.

Tests whether the variance of regression residuals depends on the values of the independent variables. If heteroskedasticity is present, OLS standard errors are biased and inference is invalid.

When to use:

After OLS regression on financial data where the volatility of returns (and hence residuals) may depend on market conditions, firm size, or other regressors.
To decide between OLS and WLS, or whether to use White/HC robust standard errors.
For validating GARCH model residuals: after fitting a GARCH model, the standardized residuals should be homoskedastic.

Mathematical formulation:

Regress squared residuals e^2 on the original regressors.
The LM statistic is n * R^2 from this auxiliary regression.
Under H0 (homoskedasticity), LM ~ chi^2(k) where k is the number of regressors.

How to interpret:

Low p-value (< 0.05): reject H0, heteroskedasticity is present. Use robust standard errors or WLS.
High p-value: no evidence of heteroskedasticity. OLS inference is valid.

Parameters:

residuals (Series | ndarray) – OLS regression residuals (1-D array or Series).
exog (DataFrame | ndarray) – Design matrix of independent variables used in the original regression (should include constant if one was used).

Returns:

lm_stat: Lagrange Multiplier statistic.
p_value: p-value from chi-squared distribution.
f_stat: F-statistic variant.
f_p_value: p-value from F-distribution.
is_heteroskedastic: bool, True if p_value < 0.05.

Return type:

Example

>>> import numpy as np, statsmodels.api as sm
>>> rng = np.random.default_rng(42)
>>> X = rng.normal(0, 1, (200, 2))
>>> X = sm.add_constant(X)
>>> y = X @ [1, 0.5, -0.3] + rng.normal(0, 1, 200)
>>> from wraquant.stats.regression import ols
>>> resid = ols(y, X, add_constant=False)["residuals"]
>>> result = breusch_pagan(resid, X)
>>> "lm_stat" in result
True

See also

white_test: More general heteroskedasticity test. durbin_watson: Test for autocorrelation instead.

white_test(residuals, exog)[source]¶

White’s test for heteroskedasticity.

White’s test is a more general heteroskedasticity test than Breusch-Pagan. It does not assume a specific functional form for the heteroskedasticity — it includes squares and cross-products of all regressors in the auxiliary regression, so it can detect non-linear forms of heteroskedasticity.

When to use:

When you want a comprehensive heteroskedasticity diagnostic that does not assume the variance is a linear function of regressors (which Breusch-Pagan assumes).
When the Breusch-Pagan test fails to reject but you still suspect non-linear heteroskedasticity.
Note: White’s test has lower power than BP when BP’s assumptions are correct, and requires more observations because it estimates more parameters.

Mathematical formulation:

Regress squared residuals e^2 on the original regressors, their squares, and all pairwise cross-products. The test statistic is n * R^2 from this auxiliary regression, which follows a chi-squared distribution under H0.

Parameters:

residuals (Series | ndarray) – OLS regression residuals (1-D array or Series).
exog (DataFrame | ndarray) – Design matrix of independent variables (should include constant if used in the original regression).

Returns:

lm_stat: White LM statistic.
p_value: p-value from chi-squared distribution.
f_stat: F-statistic variant.
f_p_value: p-value from F-distribution.
is_heteroskedastic: bool, True if p_value < 0.05.

Return type:

Example

>>> import numpy as np, statsmodels.api as sm
>>> rng = np.random.default_rng(42)
>>> X = rng.normal(0, 1, (200, 2))
>>> X = sm.add_constant(X)
>>> # Heteroskedastic errors: variance depends on X
>>> y = X @ [1, 0.5, -0.3] + rng.normal(0, 1, 200) * (1 + np.abs(X[:, 1]))
>>> from wraquant.stats.regression import ols
>>> resid = ols(y, X, add_constant=False)["residuals"]
>>> result = white_test(resid, X)
>>> "lm_stat" in result
True

See also

breusch_pagan: Simpler but less general heteroskedasticity test.

chow_test(y, X, break_point, add_constant=True)[source]¶

Chow test for structural break at a known break point.

The Chow test examines whether the regression coefficients differ between two sub-periods, i.e., whether a structural break occurred at the specified point. This is fundamental in finance for detecting regime changes, policy shifts, or market structure changes.

When to use:

To test whether a known event (e.g., a policy announcement, market crash, regulatory change) caused a structural change in the relationship between variables.
As a diagnostic for rolling regression: if the Chow test rejects stability, rolling or regime-switching models are warranted.
To validate that a backtested model’s parameters are stable across in-sample and out-of-sample periods.

Mathematical formulation:

Fit the regression on the full sample, sub-sample 1 (before break), and sub-sample 2 (after break). The F-statistic is:

\[F = \frac{(\text{RSS}_{\text{full}} - \text{RSS}_1 - \text{RSS}_2) / k}{(\text{RSS}_1 + \text{RSS}_2) / (n - 2k)}\]

where k is the number of parameters and n is the total sample size.

How to interpret:

Large F-stat (small p-value < 0.05): reject the null of stable coefficients. A structural break is detected.
Small F-stat: no evidence of a break. The relationship appears stable across the two sub-periods.

Parameters:

y (Series | ndarray) – Dependent variable (1-D array or Series).
X (DataFrame | ndarray) – Independent variables.
break_point (int) – Index (0-based row number) at which to split the sample. Must be at least k + 1 from either end.
add_constant (bool, default: True) – Whether to add an intercept to X.

Returns:

f_stat: Chow F-statistic.
p_value: p-value from the F-distribution.
break_detected: bool, True if p_value < 0.05.

Return type:

Raises:

ValueError – If break_point is too close to the endpoints.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> X = rng.normal(0, 1, (200, 1))
>>> y = np.concatenate([
...     X[:100] @ [1.0] + rng.normal(0, 0.5, 100),
...     X[100:] @ [3.0] + rng.normal(0, 0.5, 100),
... ])
>>> result = chow_test(y, X, break_point=100)
>>> result["break_detected"]
True

variance_inflation_factor(X)[source]¶

Compute the Variance Inflation Factor (VIF) for each feature.

VIF measures how much the variance of a regression coefficient is inflated due to multicollinearity with the other features. High VIF indicates that a feature is nearly a linear combination of other features, making its coefficient estimate unstable.

When to use:

Before running any multiple regression (OLS, factor model, Fama-MacBeth) to check for multicollinearity.
When regression coefficients have unexpected signs or large standard errors despite significant F-statistics.
As a feature selection diagnostic in ML pipelines.

Mathematical formulation:

For each feature X_j, regress it on all other features and compute:

\[\text{VIF}_j = \frac{1}{1 - R_j^2}\]

where R_j^2 is the R-squared from regressing X_j on the remaining features.

How to interpret:

VIF = 1: no collinearity.
VIF < 5: low collinearity, generally acceptable.
5 <= VIF < 10: moderate collinearity, warrants attention.
VIF >= 10: severe collinearity. The coefficient is poorly estimated. Consider removing the feature, combining features, or using regularization (ridge regression).

Parameters:: X (DataFrame) – DataFrame of independent variables (each column is a feature). Do not include an intercept/constant column.
Return type:: Series
Returns:: pd.Series of VIF values indexed by feature name.

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> x1 = rng.normal(0, 1, 200)
>>> x2 = x1 + rng.normal(0, 0.1, 200)  # nearly collinear with x1
>>> x3 = rng.normal(0, 1, 200)
>>> X = pd.DataFrame({"x1": x1, "x2": x2, "x3": x3})
>>> vif = variance_inflation_factor(X)
>>> vif["x1"] > 10  # collinear pair
True

See also

ols: OLS regression where VIF diagnostics are needed.

correlation_matrix(returns, method='pearson')[source]¶

Compute a correlation matrix from asset returns.

Pearson correlation measures linear dependence, Spearman measures monotonic rank dependence, and Kendall measures concordance of pairs. For financial returns, Pearson is standard but understates co-movement in the tails; Spearman and Kendall are more robust to outliers and non-linearity.

When to use:

"pearson" (default): standard linear correlation. Use for most portfolio and factor analyses.
"spearman": rank correlation. Use when you suspect non-linear but monotonic relationships, or when returns have heavy tails / outliers.
"kendall": concordance-based. More robust than Spearman for small samples. Also connects naturally to copula models (Kendall’s tau has a direct relationship to copula parameters).

Parameters:

returns (DataFrame) – DataFrame of asset returns (columns = assets).
method (str, default: 'pearson') – Correlation method – "pearson", "spearman", or "kendall".

Return type:

Returns:

Correlation matrix as a DataFrame (p x p, symmetric, diagonal entries = 1.0).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame(np.random.randn(100, 3), columns=["A", "B", "C"])
>>> corr = correlation_matrix(returns)
>>> corr.shape
(3, 3)

See also

shrunk_covariance: Regularised covariance estimation. rolling_correlation: Time-varying pairwise correlation. wraquant.risk.copulas.rank_correlation: Kendall/Spearman for

copula analysis.

shrunk_covariance(returns, method='ledoit_wolf')[source]¶

Compute a shrinkage-estimated covariance matrix.

Shrinkage estimators blend the sample covariance with a structured target to reduce estimation error, producing a better-conditioned matrix that is especially valuable when the number of assets (p) is large relative to the number of observations (T).

When to use:

Always prefer shrinkage over the raw sample covariance for portfolio optimisation. The improvement is largest when p/T is close to or exceeds 1 (e.g., 500 stocks with 252 daily observations).

"ledoit_wolf" (default): analytically optimal shrinkage toward a structured target. Best general-purpose choice. Automatically determines the optimal shrinkage intensity.
"oas" (Oracle Approximating Shrinkage): assumes the underlying distribution is Gaussian and computes the oracle- approximating shrinkage intensity. Slightly better than Ledoit-Wolf when normality holds.
"basic": simple shrinkage toward the diagonal with a fixed (non-optimal) shrinkage coefficient. Use only as a baseline.

Mathematical formulation:

Sigma_shrunk = (1 - alpha) * S + alpha * F

where S is the sample covariance, F is the shrinkage target (e.g., identity or diagonal), and alpha is the shrinkage intensity (0 = no shrinkage, 1 = full shrinkage to target).

How to interpret:

The returned matrix is guaranteed positive semi-definite. Its eigenvalues are more dispersed than the sample covariance (less extreme), leading to more stable portfolio weights. Compare the condition number (ratio of max to min eigenvalue) before and after shrinkage to see the regularisation effect.

Parameters:

returns (DataFrame) – DataFrame of asset returns (columns = assets).
method (str, default: 'ledoit_wolf') – Shrinkage method – "ledoit_wolf" (default), "oas", or "basic".

Return type:

Returns:

Shrunk covariance matrix as a DataFrame (p x p).

Raises:

ValueError – If method is not recognized.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame(np.random.randn(100, 5), columns=list("ABCDE"))
>>> cov = shrunk_covariance(returns, method="ledoit_wolf")
>>> cov.shape
(5, 5)

See also

correlation_matrix: Correlation (standardised covariance). wraquant.ml.preprocessing.denoised_correlation: Random matrix

theory-based denoising.

References

Ledoit & Wolf (2004), “A well-conditioned estimator for large-dimensional covariance matrices”
Chen et al. (2010), “Shrinkage Algorithms for MMSE Covariance Estimation”

rolling_correlation(x, y, window)[source]¶

Compute rolling Pearson correlation between two series.

Rolling correlation reveals how the linear relationship between two assets evolves over time. Stable correlation is a key assumption in portfolio construction; large swings in rolling correlation indicate that static portfolio weights may be suboptimal.

When to use:

Use rolling correlation to: - Monitor diversification benefit over time (correlation

rising toward 1.0 means diversification is eroding).

Detect correlation regime changes for pairs trading or hedging ratio adjustment.
Validate the stationarity assumption of portfolio optimisation inputs.

How to interpret:

Values near +1.0: strong positive co-movement (little diversification benefit).
Values near 0.0: approximately uncorrelated.
Values near -1.0: strong negative co-movement (excellent diversification or natural hedge).
Spikes toward +1.0 during sell-offs are typical (“correlation goes to 1 in a crisis”).

Parameters:

x (Series) – First return series.
y (Series) – Second return series (same index).
window (int) – Rolling window size (e.g., 60 for ~3 months of daily data).

Return type:

Returns:

Rolling Pearson correlation series. First window - 1 values are NaN.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> x = pd.Series(np.random.randn(200))
>>> y = pd.Series(0.5 * x + np.random.randn(200) * 0.5)
>>> rc = rolling_correlation(x, y, window=60)
>>> rc.dropna().iloc[0] > 0
True

See also

correlation_matrix: Full cross-asset correlation matrix. wraquant.risk.dcc.rolling_correlation_dcc: DCC-GARCH-based

dynamic correlation.

partial_correlation(data)[source]¶

Compute the partial correlation matrix, controlling for all other variables.

Partial correlation measures the linear association between two variables after removing the linear effect of all other variables in the dataset. This is essential in finance for understanding direct relationships between assets, factors, or macro variables — as opposed to associations that are mediated through a common driver.

When to use:

Use partial correlation when you suspect that the observed correlation between two assets (or factors) is driven by a shared exposure to a third variable. For example, two energy stocks may appear highly correlated, but partial correlation can reveal that after controlling for oil prices, the direct relationship is weak.

Mathematical formulation:

For each pair (i, j), regress both X_i and X_j on all remaining variables, then compute the Pearson correlation of the residuals:

\[\rho_{ij \cdot \text{rest}} = \text{corr}(e_i, e_j)\]

where e_i = X_i - \hat{X}_i is the residual from regressing X_i on all other columns.

Equivalently, partial correlations can be obtained from the inverse of the correlation matrix (the precision matrix):

\[\rho_{ij \cdot \text{rest}} = -\frac{P_{ij}}{\sqrt{P_{ii} P_{jj}}}\]

where P = R^{-1} is the precision matrix.

How to interpret:

Values near 0 indicate no direct linear relationship once shared drivers are removed.
A large drop from raw correlation to partial correlation signals that the association is mostly indirect (mediated).
The diagonal is always 1.0.

Parameters:: data (DataFrame) – DataFrame with columns as variables (assets, factors, etc.) and rows as observations. Must have at least 3 columns.
Return type:: DataFrame
Returns:: Partial correlation matrix as a DataFrame (p x p, symmetric, diagonal = 1.0).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> z = np.random.randn(200)
>>> data = pd.DataFrame({
...     "A": z + np.random.randn(200) * 0.3,
...     "B": z + np.random.randn(200) * 0.3,
...     "C": np.random.randn(200),
... })
>>> pcorr = partial_correlation(data)
>>> pcorr.shape
(3, 3)

See also

correlation_matrix: Standard (marginal) correlation matrix. mutual_information: Non-linear dependence measure.

distance_correlation(x, y)[source]¶

Compute the Brownian distance correlation between two variables.

Distance correlation (Szekely, Rizzo & Bakirov, 2007) is a measure of dependence that equals zero if and only if the two variables are independent — unlike Pearson correlation, which only captures linear dependence. This makes it invaluable for detecting nonlinear relationships in financial data (e.g., option-like payoffs, regime- dependent correlations, or tail dependence).

When to use:

You suspect a nonlinear relationship that Pearson/Spearman will miss (e.g., a U-shaped or threshold relationship).
You want a single-number summary of any type of dependence.
You need a test statistic for independence that is consistent against all alternatives with finite first moments.

Mathematical formulation:

Compute the pairwise Euclidean distance matrices a_{kl} = |X_k - X_l| and b_{kl} = |Y_k - Y_l|.
Double-center each matrix: A_{kl} = a_{kl} - \bar{a}_{k\cdot} - \bar{a}_{\cdot l} + \bar{a}_{\cdot\cdot}
Distance covariance squared: \text{dCov}^2(X, Y) = \frac{1}{n^2} \sum_{k,l} A_{kl} B_{kl}
Distance correlation: \text{dCor}(X, Y) = \frac{\text{dCov}(X, Y)}{\sqrt{\text{dVar}(X) \cdot \text{dVar}(Y)}}

How to interpret:

0.0: independence (no dependence of any kind).
1.0: perfect dependence (deterministic relationship).
Values between 0 and 1 indicate partial dependence.
Distance correlation >= |Pearson correlation|, so it always detects at least as much dependence.

Parameters:

x (Series | ndarray) – First variable (1-D array or Series).
y (Series | ndarray) – Second variable (1-D array or Series, same length).

Return type:

Returns:

Distance correlation as a float in [0, 1].

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 200)
>>> y = x ** 2 + rng.normal(0, 0.3, 200)  # nonlinear
>>> dcor = distance_correlation(x, y)
>>> dcor > 0.3  # detects nonlinear dependence
True

References

Szekely, G. J., Rizzo, M. L. & Bakirov, N. K. (2007). “Measuring and testing dependence by correlation of distances.” Annals of Statistics, 35(6), 2769-2794.

See also

correlation_matrix: Linear (Pearson) correlation. mutual_information: Information-theoretic dependence measure.

kendall_tau(x, y)[source]¶

Compute Kendall’s tau-b rank correlation coefficient with p-value.

Kendall’s tau measures the ordinal association between two variables. It counts the number of concordant and discordant pairs: a pair (x_i, y_i), (x_j, y_j) is concordant if the ranks agree and discordant if they disagree. The tau-b variant adjusts for ties.

When to use:

When data has heavy tails or outliers that distort Pearson correlation.
For small samples where Spearman is less reliable.
When you need a rank-based measure that connects naturally to copula parameters (Kendall’s tau has a one-to-one mapping to the parameter of many copula families).

Mathematical formulation:

\[\tau_b = \frac{C - D}{\sqrt{(C + D + T_x)(C + D + T_y)}}\]

where C = concordant pairs, D = discordant pairs, T_x and T_y = pairs tied only on x or y.

How to interpret:

+1: perfect concordance (monotonically increasing relationship).
-1: perfect discordance (monotonically decreasing).
0: no ordinal association.
|tau| > 0.3 is generally considered a moderate association in financial data.
The p-value tests H0: tau = 0 (independence).

Parameters:

x (Series | ndarray) – First variable (1-D array or Series).
y (Series | ndarray) – Second variable (1-D array or Series, same length).

Returns:

tau: Kendall’s tau-b statistic.
p_value: two-sided p-value for H0: tau = 0.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 100)
>>> y = 0.8 * x + rng.normal(0, 0.5, 100)
>>> result = kendall_tau(x, y)
>>> result["tau"] > 0
True
>>> result["p_value"] < 0.05
True

See also

correlation_matrix: Pearson/Spearman/Kendall full matrix. distance_correlation: Nonlinear dependence measure.

mutual_information(x, y, n_bins=20, method='binning')[source]¶

Estimate mutual information between two continuous variables.

Mutual information (MI) quantifies the amount of information obtained about one variable by observing the other. Unlike correlation, MI captures any type of statistical dependence — linear, nonlinear, or even purely distributional.

When to use:

Feature selection for ML-based trading models: MI identifies features with any predictive signal, not just linear ones.
Comparing the information content of different alpha signals.
Measuring the ``true’’ dependence between returns and macro indicators that may have complex, non-monotonic relationships.

Mathematical formulation:

\[I(X; Y) = \sum_{x} \sum_{y} p(x, y) \log \frac{p(x, y)}{p(x) p(y)}\]

For continuous variables, the sums become integrals. The "binning" method discretises both variables into n_bins bins and computes MI on the resulting contingency table. The "kde" method uses kernel density estimation for the joint and marginal densities.

How to interpret:

MI = 0: independence (knowing X tells you nothing about Y).
MI > 0: some dependence exists.
MI is measured in nats (when using natural log) and is non-negative.
There is no upper bound in general, but normalised MI (MI / sqrt(H(X)*H(Y))) can be used for comparisons.

Parameters:

x (Series | ndarray) – First continuous variable.
y (Series | ndarray) – Second continuous variable (same length).
n_bins (int, default: 20) – Number of bins for discretisation ("binning" method). More bins capture finer structure but need more data.
method (str, default: 'binning') – Estimation method – "binning" (default) or "kde".

Return type:

Returns:

Estimated mutual information in nats (>= 0).

Raises:

ValueError – If method is not recognized.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 500)
>>> y = x + rng.normal(0, 0.5, 500)
>>> mi = mutual_information(x, y)
>>> mi > 0
True

See also

distance_correlation: Another non-linear dependence measure. correlation_matrix: Linear dependence only.

correlation_significance(x, y, method='pearson', confidence=0.95)[source]¶

Test whether the correlation between two variables is significantly non-zero.

Computes the sample correlation, performs a t-test for the null hypothesis H0: rho = 0, and constructs a confidence interval using Fisher’s z-transformation.

When to use:

To confirm that an observed correlation is statistically significant and not just sampling noise.
To obtain confidence intervals for reporting correlation estimates with uncertainty.
Before using a correlation estimate in portfolio construction or risk models — insignificant correlations may be unreliable.

Mathematical formulation:

Test statistic:

\[t = r \sqrt{\frac{n - 2}{1 - r^2}}\]

which follows a t-distribution with n - 2 degrees of freedom under H0.

Confidence interval via Fisher z-transformation:

\[z = \text{arctanh}(r), \quad SE = \frac{1}{\sqrt{n - 3}}\]

The interval [z - z_{\alpha/2} \cdot SE, z + z_{\alpha/2} \cdot SE] is back-transformed via tanh() to the correlation scale.

Parameters:

x (Series | ndarray) – First variable.
y (Series | ndarray) – Second variable (same length).
method (str, default: 'pearson') – Correlation method – "pearson" (default) or "spearman".
confidence (float, default: 0.95) – Confidence level for the interval (default 0.95).

Returns:

r: sample correlation coefficient.
t_stat: t-test statistic.
p_value: two-sided p-value for H0: rho = 0.
ci_lower: lower bound of the confidence interval.
ci_upper: upper bound of the confidence interval.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 100)
>>> y = 0.5 * x + rng.normal(0, 1, 100)
>>> result = correlation_significance(x, y)
>>> result["p_value"] < 0.05
True
>>> result["ci_lower"] < result["r"] < result["ci_upper"]
True

See also

correlation_matrix: Compute correlations without significance test. kendall_tau: Rank correlation with p-value.

minimum_spanning_tree_correlation(corr_matrix)[source]¶

Compute the minimum spanning tree (MST) of a correlation matrix.

The MST is a connected, acyclic subgraph that connects all assets with the minimum total distance, where distance is derived from correlation. It reveals the hierarchical structure of the market: which assets are the most “central” and how clusters of correlated assets are organized.

When to use:

To visualise market structure and identify clusters of related assets (sectors, factor groups, regimes).
As input to hierarchical risk parity (HRP) portfolio construction (Lopez de Prado, 2016).
To detect changes in market structure over time by comparing MSTs across different periods.
To reduce dimensionality of the correlation matrix for network-based analysis.

Mathematical formulation:

The correlation matrix is converted to a distance matrix:

\[d_{ij} = \sqrt{2(1 - \rho_{ij})}\]

This metric satisfies the triangle inequality and maps perfect correlation (rho = 1) to zero distance and zero correlation (rho = 0) to distance sqrt(2).

Prim’s or Kruskal’s algorithm is then applied to find the MST of the complete weighted graph.

How to interpret:

The returned adjacency matrix has non-zero entries only for edges in the MST. The values are the correlation-derived distances. Assets connected by short edges are highly correlated; the “hub” asset with the most edges is the most central.

Parameters:: corr_matrix (DataFrame) – Correlation matrix as a DataFrame (p x p, symmetric).
Return type:: DataFrame
Returns:: Adjacency matrix of the MST as a DataFrame (p x p). Non-zero entries indicate edges in the tree, with values equal to the distance sqrt(2 * (1 - rho)).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame(np.random.randn(100, 4), columns=list("ABCD"))
>>> corr = returns.corr()
>>> mst = minimum_spanning_tree_correlation(corr)
>>> mst.shape
(4, 4)
>>> (mst.values > 0).sum()  # MST has p-1 edges, each counted twice
6

References

Mantegna, R. N. (1999). “Hierarchical structure in financial markets.”
Lopez de Prado, M. (2016). “Building diversified portfolios that outperform out-of-sample.”

See also

correlation_matrix: Compute the input correlation matrix. shrunk_covariance: Regularised covariance for more stable MSTs.

tail_dependence_coefficient(x, y, threshold=0.05)[source]¶

Estimate upper and lower tail dependence coefficients from empirical data.

Tail dependence measures the probability that one variable is extremely large (small) given that the other is also extremely large (small). This is crucial for understanding joint tail risk in portfolios — standard correlation says nothing about co-movement in the tails.

When to use:

To quantify the risk of joint extreme losses in a portfolio.
To assess whether diversification benefits disappear during market crashes (asymmetric tail dependence).
To select the appropriate copula family: Gaussian copulas have zero tail dependence, while Clayton (lower) and Gumbel (upper) copulas can model it.
To compare the tail behaviour of different asset pairs.

Mathematical formulation:

The upper tail dependence coefficient is:

\[\lambda_U = \lim_{q \to 1} P(Y > F_Y^{-1}(q) \mid X > F_X^{-1}(q))\]

The lower tail dependence coefficient is:

\[\lambda_L = \lim_{q \to 0} P(Y \le F_Y^{-1}(q) \mid X \le F_X^{-1}(q))\]

We estimate these empirically using the rank-transformed data (pseudo-observations) and counting joint exceedances.

How to interpret:

lambda = 0: no tail dependence (e.g., Gaussian copula). Diversification holds in the tails.
lambda > 0: positive tail dependence. Extreme events tend to happen together.
lambda_L > lambda_U: lower tail dependence is stronger than upper (common in equity markets — crashes are more contagious than rallies).
Values are bounded in [0, 1].

Parameters:

x (Series | ndarray) – First variable (1-D array or Series).
y (Series | ndarray) – Second variable (1-D array or Series, same length).
threshold (float, default: 0.05) – Quantile threshold for defining “extreme” (default 0.05, meaning the top/bottom 5%).

Returns:

upper_lambda: estimated upper tail dependence coefficient.
lower_lambda: estimated lower tail dependence coefficient.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> # Gaussian data has zero tail dependence
>>> x = rng.normal(0, 1, 5000)
>>> y = 0.7 * x + rng.normal(0, 0.71, 5000)
>>> result = tail_dependence_coefficient(x, y)
>>> 0 <= result["upper_lambda"] <= 1
True
>>> 0 <= result["lower_lambda"] <= 1
True

References

Joe, H. (2014). Dependence Modeling with Copulas, Ch. 2.
McNeil et al. (2015). Quantitative Risk Management, Ch. 7.

See also

copula_selection: Fit copulas that model tail dependence. rank_correlation_matrix: Rank-based dependence matrix.

copula_selection(x, y)[source]¶

Fit multiple copula families and select the best by AIC.

Copulas separate the marginal distributions from the dependence structure, allowing flexible modelling of how two variables move together. This function fits several parametric copula families and ranks them by AIC to identify the best model for the data.

When to use:

To model joint distributions for portfolio risk (e.g., joint simulation of asset returns for VaR/CVaR).
To capture tail dependence or asymmetric dependence that Gaussian models miss.
To select the appropriate copula for bivariate analysis before using it in a larger risk framework.

Copula families fitted:

Gaussian: symmetric, zero tail dependence.
Student-t (approximated): symmetric, positive tail dependence in both tails.
Clayton: lower tail dependence (joint crashes).
Gumbel: upper tail dependence (joint rallies).

Mathematical formulation:

A copula C(u, v) is a joint CDF on [0,1]^2 whose marginals are uniform. By Sklar’s theorem, any joint distribution can be written as:

\[F(x, y) = C(F_X(x), F_Y(y))\]

Each copula family has a parameter theta that controls the strength and shape of dependence. We estimate theta by maximum likelihood on the pseudo-observations.

Parameters:

x (Series | ndarray) – First variable (1-D array or Series).
y (Series | ndarray) – Second variable (1-D array or Series, same length).

Returns:

best_copula: name of the best-fitting copula.
all_fits: DataFrame with columns copula, parameter, log_likelihood, aic, sorted by AIC ascending.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 500)
>>> y = 0.7 * x + rng.normal(0, 0.71, 500)
>>> result = copula_selection(x, y)
>>> isinstance(result["all_fits"], pd.DataFrame)
True
>>> result["best_copula"] in ["gaussian", "student_t", "clayton", "gumbel"]
True

References

Joe, H. (2014). Dependence Modeling with Copulas.
Nelsen, R. B. (2006). An Introduction to Copulas.

See also

tail_dependence_coefficient: Empirical tail dependence. rank_correlation_matrix: Rank-based dependence matrix.

rank_correlation_matrix(data, method='spearman')[source]¶

Compute the Spearman rank correlation matrix.

Spearman correlation assesses monotonic relationships by computing Pearson correlation on the rank-transformed data. It is more robust to outliers and non-linearity than Pearson correlation, making it well-suited for financial data with heavy tails.

When to use:

When the relationship between variables is monotonic but not necessarily linear.
When data contains outliers that would distort Pearson correlation.
For copula parameter estimation (Spearman’s rho has a direct relationship to many copula parameters).
As a robustness check on Pearson correlation: if the two differ substantially, the relationship may be nonlinear.

Mathematical formulation:

\[\rho_S(X, Y) = \text{Pearson}(\text{rank}(X), \text{rank}(Y))\]

Equivalently:

\[\rho_S = 1 - \frac{6 \sum d_i^2}{n(n^2 - 1)}\]

where d_i is the difference in ranks.

How to interpret:

Same scale as Pearson: [-1, 1].
rho_S > rho_P: the relationship is stronger in the ranks than in the raw values (concave/convex relationship).
rho_S ≈ rho_P: the relationship is approximately linear.

Parameters:

data (DataFrame) – DataFrame with columns as variables and rows as observations.
method (str, default: 'spearman') – Rank correlation method – "spearman" (default) or "kendall".

Return type:

Returns:

Rank correlation matrix as a DataFrame (p x p).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> data = pd.DataFrame(np.random.randn(100, 4), columns=list("ABCD"))
>>> rcm = rank_correlation_matrix(data)
>>> rcm.shape
(4, 4)

See also

correlation_matrix: Pearson correlation matrix. partial_correlation: Partial correlation controlling for others.

concordance_index(predicted, observed)[source]¶

Compute Harrell’s concordance index (C-index).

The C-index measures the probability that for a random pair of observations, the one with the higher predicted value also has the higher observed value. It is a generalised rank correlation measure widely used in survival analysis, credit risk, and model validation.

When to use:

To evaluate the discriminatory power of a predictive model (e.g., a credit scoring model, a default probability model, or an alpha signal).
When you care about ordinal accuracy (ranking) rather than calibration (magnitude).
As an alternative to AUC for continuous outcomes (the C-index generalises AUC to continuous data).

Mathematical formulation:

\[C = \frac{\text{concordant pairs}}{\text{concordant + discordant pairs}}\]

A pair (i, j) is concordant if the predicted and observed orderings agree: (\hat{y}_i > \hat{y}_j \text{ and } y_i > y_j) or (\hat{y}_i < \hat{y}_j \text{ and } y_i < y_j).

How to interpret:

C = 0.5: random (no predictive power).
C = 1.0: perfect concordance (perfect ranking).
C < 0.5: worse than random (model has the sign wrong).
C > 0.7: generally considered acceptable in finance.
C > 0.8: strong discriminatory power.

Parameters:

predicted (Series | ndarray) – Predicted values or scores.
observed (Series | ndarray) – Observed/actual values.

Return type:

Returns:

Concordance index as a float in [0, 1].

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> predicted = rng.normal(0, 1, 100)
>>> observed = predicted + rng.normal(0, 0.5, 100)
>>> c = concordance_index(predicted, observed)
>>> c > 0.7  # good concordance
True

References

Harrell, F. E., Lee, K. L. & Mark, D. B. (1996). “Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors.” Statistics in Medicine, 15(4), 361-387.

See also

kendall_tau: Rank correlation (related to C-index).

fit_distribution(data, dist='norm')[source]¶

Fit a parametric distribution to data.

Parameters:

data (Series) – Data series to fit.
dist (str, default: 'norm') – Name of a scipy.stats distribution (e.g., "norm", "t", "lognorm").

Return type:

Returns:

Dictionary with params (tuple of fitted parameters), ks_statistic, and ks_pvalue from a Kolmogorov-Smirnov goodness-of-fit test.

Raises:

AttributeError – If dist is not a valid scipy distribution.

fit_stable_distribution(data)[source]¶

Fit a stable (Levy) distribution to data.

Uses scipy’s levy_stable distribution to estimate the four parameters: alpha (stability), beta (skewness), loc, and scale.

Parameters:

data (Series | ndarray) – Data array or series.

Returns:

alpha: stability parameter (0, 2].
beta: skewness parameter [-1, 1].
loc: location parameter.
scale: scale parameter.
ks_statistic: Kolmogorov-Smirnov statistic.
ks_pvalue: KS test p-value.

Return type:

tail_ratio(returns, quantile=0.05)[source]¶

Compute the tail ratio (right tail / left tail).

A tail ratio > 1 indicates a fatter right tail (more extreme gains) relative to the left tail.

Parameters:

returns (Series) – Return series.
quantile (float, default: 0.05) – Quantile for tail measurement (default 5%).

Return type:

Returns:

Tail ratio as a float.

tail_index(data, method='hill', threshold_quantile=0.9)[source]¶

Estimate the tail index of a distribution.

The tail index characterises the heaviness of distribution tails. A finite tail index indicates power-law tails (Pareto-like).

Parameters:

data (Series | ndarray) – Data array or series.
method (str, default: 'hill') – Estimation method. One of "hill" (Hill estimator), "pickands" (Pickands estimator), or "moment" (moment estimator of Dekkers-Einmahl-de Haan).
threshold_quantile (float, default: 0.9) – Quantile above which tail observations are used (default 0.9, i.e. top 10%).

Returns:

tail_index: estimated tail index (xi).
method: method used.
n_tail: number of observations in the tail.

Return type:

Raises:

ValueError – If method is not one of the supported estimators.

hurst_exponent(data)[source]¶

Estimate the Hurst exponent via rescaled range (R/S) analysis.

The Hurst exponent characterises the long-term memory of a series:

H < 0.5: mean-reverting
H = 0.5: random walk
H > 0.5: trending / persistent

Parameters:: data (Series) – Time series (prices or returns).
Return type:: float
Returns:: Estimated Hurst exponent as a float.

qqplot_data(data, dist='norm')[source]¶

Generate quantile-quantile plot data.

Computes theoretical and sample quantiles for constructing a Q-Q plot against a reference distribution.

Parameters:

data (Series | ndarray) – Data array or series.
dist (str, default: 'norm') – Name of a scipy.stats distribution to use as the theoretical reference (default "norm").

Returns:

theoretical_quantiles: array of theoretical quantiles.
sample_quantiles: array of ordered sample values.
slope: slope of the best-fit line through the Q-Q plot.
intercept: intercept of the best-fit line.

Return type:

jarque_bera(data)[source]¶

Perform the Jarque-Bera test for normality.

Tests the null hypothesis that the data is normally distributed, based on sample skewness and kurtosis.

Parameters:

data (Series | ndarray) – Data array or series.

Returns:

statistic: Jarque-Bera test statistic.
p_value: p-value of the test.
skewness: sample skewness.
kurtosis: sample excess kurtosis.

Return type:

kolmogorov_smirnov(data, dist='norm')[source]¶

Perform the Kolmogorov-Smirnov goodness-of-fit test.

Tests the null hypothesis that data was drawn from the specified distribution. The distribution parameters are first estimated via MLE.

Parameters:

data (Series | ndarray) – Data array or series.
dist (str, default: 'norm') – Name of a scipy.stats distribution (default "norm").

Returns:

statistic: KS test statistic.
p_value: p-value of the test.
dist: distribution name tested.
params: fitted distribution parameters.

Return type:

anderson_darling(data, dist='norm')[source]¶

Perform the Anderson-Darling goodness-of-fit test.

The Anderson-Darling test is more sensitive to deviations in the tails than the Kolmogorov-Smirnov test, making it more appropriate for financial data where tail behaviour matters most (e.g., VaR and CVaR estimation).

Parameters:

data (Series | ndarray) – Data array or series.
dist (str, default: 'norm') – Distribution to test against. Supported values depend on scipy.stats.anderson and include "norm", "expon", "logistic", "gumbel", "gumbel_l", "gumbel_r".

Returns:

statistic: Anderson-Darling test statistic.
critical_values: array of critical values for each significance level.
significance_levels: corresponding significance levels (%).

Return type:

Example

>>> import numpy as np
>>> data = np.random.default_rng(42).normal(0, 1, 1000)
>>> anderson_darling(data)

best_fit_distribution(data, candidates=None)[source]¶

Rank candidate distributions by goodness of fit.

Fits multiple parametric distributions to the data and ranks them by AIC and KS/AD statistics. Use this to choose the best model for return distributions when the assumption of normality fails (which it usually does in finance).

Parameters:

data (Series | ndarray) – Data array or series.
candidates (list[str] | None, default: None) – List of scipy.stats distribution names to test. Defaults to ["norm", "t", "skewnorm", "gennorm", "nct", "johnsonsu"] – a set well-suited for financial returns.

Returns:

distribution, params, ks_statistic, ad_statistic, aic, sorted by AIC (ascending).

Return type:

Example

>>> import numpy as np
>>> data = np.random.default_rng(42).standard_t(df=5, size=1000)
>>> best_fit_distribution(data)

kernel_density_estimate(data, n_points=200, bandwidth='scott')[source]¶

Kernel density estimation using Gaussian kernels.

Non-parametric density estimation that makes no assumptions about the underlying distribution shape. Useful for visualizing return distributions, computing non-parametric VaR, and comparing regime-specific densities.

Parameters:

data (ndarray | Series) – Sample data (1D array of returns or prices).
n_points (int, default: 200) – Number of evaluation points for the density curve.
bandwidth (str | float, default: 'scott') – Bandwidth method (“scott”, “silverman”) or float.

Returns:

x – Evaluation points (n_points,).
density – Estimated density values (n_points,).
bandwidth – Bandwidth used.
mode – Location of peak density.
cdf – Cumulative distribution values (n_points,).

Return type:

dict[str, ndarray]

Example

>>> from wraquant.stats.distributions import kernel_density_estimate
>>> kde = kernel_density_estimate(returns, n_points=500)
>>> var_95 = kde['x'][np.searchsorted(kde['cdf'], 0.05)]

mad(data, scale='normal')[source]¶

Compute the Median Absolute Deviation (MAD).

MAD is a robust measure of dispersion. Unlike standard deviation, it is not influenced by a few extreme values, making it ideal for financial return distributions with fat tails.

Parameters:

data (Series | ndarray) – Data series or array.
scale (str, default: 'normal') – Scaling factor. Use "normal" (default) so that the result is consistent with standard deviation for normally distributed data. Use None for the raw MAD.

Return type:

Returns:

MAD as a float.

Example

>>> import pandas as pd
>>> returns = pd.Series([0.01, 0.02, -0.01, -0.05, 0.10])
>>> mad(returns)

winsorize(data, lower=0.05, upper=0.05)[source]¶

Cap extreme values at given percentiles (Winsorization).

Winsorization limits extreme values to reduce the influence of outliers without removing observations. This is preferable to trimming when you want to keep the same sample size.

Parameters:

data (Series | ndarray) – Data series or array.
lower (float, default: 0.05) – Fraction to clip on the lower tail (default 5%).
upper (float, default: 0.05) – Fraction to clip on the upper tail (default 5%).

Return type:

Series | ndarray

Returns:

Winsorized data, same type as input.

Example

>>> import pandas as pd
>>> returns = pd.Series([0.01, 0.02, -0.50, 0.03, 0.80])
>>> winsorize(returns, lower=0.1, upper=0.1)

trimmed_mean(data, proportiontocut=0.05)[source]¶

Compute the trimmed mean, excluding extreme observations.

The trimmed mean removes a fraction of the highest and lowest values before computing the average. Use it when the mean is distorted by outliers (e.g., flash-crash returns).

Parameters:

data (Series | ndarray) – Data series or array.
proportiontocut (float, default: 0.05) – Fraction to cut from each tail (default 5%).

Return type:

Returns:

Trimmed mean as a float.

Example

>>> import numpy as np
>>> data = np.array([1, 2, 3, 4, 100])
>>> trimmed_mean(data, proportiontocut=0.2)

trimmed_std(data, proportiontocut=0.05)[source]¶

Compute the standard deviation after trimming extreme values.

Combines trimming with standard deviation computation for a measure of dispersion that is less sensitive to outliers than standard std but retains more information than MAD.

Parameters:

data (Series | ndarray) – Data series or array.
proportiontocut (float, default: 0.05) – Fraction to cut from each tail (default 5%).

Return type:

Returns:

Trimmed standard deviation as a float.

Example

>>> import numpy as np
>>> data = np.array([1, 2, 3, 4, 100])
>>> trimmed_std(data, proportiontocut=0.2)

robust_zscore(data)[source]¶

Compute robust z-scores using median and MAD.

Standard z-scores (x - mean) / std are heavily influenced by outliers. Robust z-scores replace mean with median and std with MAD, providing a more reliable outlier detection metric for financial data.

Parameters:: data (Series | ndarray) – Data series or array.
Return type:: Series
Returns:: Robust z-scores as a pd.Series.

Example

>>> import pandas as pd
>>> returns = pd.Series([0.01, 0.02, -0.01, -0.05, 0.50])
>>> robust_zscore(returns)

robust_covariance(data, support_fraction=None)[source]¶

Estimate a robust covariance matrix via Minimum Covariance Determinant.

The MCD estimator finds the subset of observations (of a given fraction) whose classical covariance has the smallest determinant. This makes it highly resistant to outliers – essential when computing portfolio covariance from return data that may contain erroneous prints or fat-tailed events.

Parameters:

data (DataFrame) – DataFrame of asset returns (columns = assets).
support_fraction (float | None, default: None) – Fraction of data to use in support (default None lets sklearn choose).

Returns:

covariance: robust covariance matrix (np.ndarray).
location: robust location estimate (np.ndarray).
support_fraction: fraction of data used.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> returns = pd.DataFrame(np.random.randn(100, 3), columns=['A', 'B', 'C'])
>>> robust_covariance(returns)

huber_mean(data, delta=1.5, max_iter=50, tol=1e-08)[source]¶

Compute the Huber M-estimator of location.

The Huber estimator behaves like the mean for observations within delta MAD-scaled deviations of the center, but limits the influence of observations beyond that threshold via iteratively reweighted least squares. It provides a smooth trade-off between efficiency (mean) and robustness (median).

Parameters:

data (Series | ndarray) – Data series or array.
delta (float, default: 1.5) – Threshold parameter controlling robustness. Smaller values give more robustness (closer to median). Default 1.5 is a standard choice.
max_iter (int, default: 50) – Maximum number of IRLS iterations.
tol (float, default: 1e-08) – Convergence tolerance for the location estimate.

Return type:

Returns:

Huber location estimate as a float.

Example

>>> import numpy as np
>>> data = np.array([1, 2, 3, 4, 100])
>>> huber_mean(data, delta=1.5)

outlier_detection(data, method='mad', threshold=3.0)[source]¶

Flag outliers using a robust detection method.

Outlier detection is critical in finance for identifying data errors (bad ticks), extreme events, or contaminated observations before computing risk metrics.

Parameters:

data (Series | ndarray) – Data series or array.
method (Literal['mad', 'iqr', 'grubbs'], default: 'mad') –
Detection method: - "mad": Median Absolute Deviation (default). Flag

points whose robust z-score exceeds threshold. Best general-purpose choice for financial data.
- "iqr": Interquartile Range. Flag points outside [Q1 - threshold*IQR, Q3 + threshold*IQR]. Classic box-plot method.
- "grubbs": Grubbs’ test for a single outlier. Tests whether the most extreme value is an outlier assuming approximate normality.
threshold (float, default: 3.0) – Sensitivity parameter (default 3.0). For MAD this is the z-score cutoff; for IQR it is the multiplier.

Returns:

outliers: boolean array (True = outlier).
n_outliers: count of flagged outliers.
method: method used.

Return type:

Raises:

ValueError – If method is not recognized.

Example

>>> import pandas as pd
>>> returns = pd.Series([0.01, 0.02, -0.01, -0.50, 0.03])
>>> result = outlier_detection(returns, method="mad")
>>> result["n_outliers"]

engle_granger(y1, y2, max_lag=None)[source]¶

Engle-Granger two-step cointegration test.

Regresses y1 on y2 via OLS, then tests the residuals for a unit root using the Augmented Dickey-Fuller test.

Parameters:

y1 (Series) – First price series.
y2 (Series) – Second price series.
max_lag (int | None, default: None) – Maximum number of lags for the ADF test. When None, adfuller selects lags automatically via AIC.

Return type:

Returns:

Dictionary with statistic (ADF test statistic), p_value, is_cointegrated (at 5 % significance), hedge_ratio (OLS slope coefficient), and residuals (pd.Series).

johansen(data, det_order=0, k_ar_diff=1)[source]¶

Johansen cointegration test for multiple time series.

Requires the timeseries optional dependency group (provides statsmodels.tsa.vector_ar.vecm).

Parameters:

data (DataFrame) – DataFrame of price series (columns = assets).
det_order (int, default: 0) – Deterministic term order. -1 for no deterministic term, 0 for constant, 1 for linear trend.
k_ar_diff (int, default: 1) – Number of lagged differences in the model.

Return type:

Returns:

Dictionary with trace_stats (array of trace statistics), eigenvalues, critical_values_95 (95 % critical values), and n_cointegrating (number of cointegrating relationships at the 5 % level).

half_life(spread)[source]¶

Estimate the half-life of mean reversion for a spread series.

Fits an OLS regression of the change in spread on the lagged spread level: delta_spread = phi * spread_lag + eps. The half-life is -log(2) / log(1 + phi).

Parameters:: spread (Series) – Spread (or residual) series.
Return type:: float
Returns:: Half-life in the same time units as the spread index. Returns float('inf') when the spread is not mean-reverting.

spread(y1, y2, hedge_ratio=None)[source]¶

Compute the spread between two price series.

When hedge_ratio is None it is estimated via OLS.

Parameters:

y1 (Series) – First price series (the dependent variable).
y2 (Series) – Second price series (the independent variable).
hedge_ratio (float | None, default: None) – Explicit hedge ratio. If None, the ratio is estimated from the data using OLS.

Return type:

Returns:

Spread series (y1 - hedge_ratio * y2).

zscore_signal(spread, window=20)[source]¶

Compute a rolling z-score of the spread for trading signals.

Parameters:

spread (Series) – Spread series.
window (int, default: 20) – Rolling window size for mean and standard deviation.

Return type:

Returns:

Rolling z-score series.

hedge_ratio(y1, y2, method='ols')[source]¶

Estimate the hedge ratio between two price series.

Parameters:

y1 (Series) – First (dependent) price series.
y2 (Series) – Second (independent) price series.
method (str, default: 'ols') – Estimation method — "ols" for ordinary least squares or "tls" for total least squares (orthogonal regression).

Return type:

Returns:

Hedge ratio as a float.

Raises:

ValueError – If method is not recognized.

pairs_backtest_signals(spread, entry_z=2.0, exit_z=0.5)[source]¶

Generate pairs trading signals based on z-score thresholds.

The strategy goes short the spread when z-score > entry_z, goes long when z-score < -entry_z, and exits when the z-score crosses back inside [-exit_z, exit_z].

Parameters:

spread (Series) – Spread series (raw, not z-scored).
entry_z (float, default: 2.0) – Z-score threshold for entry (absolute value).
exit_z (float, default: 0.5) – Z-score threshold for exit (absolute value).

Return type:

Returns:

Signal series with values in {-1, 0, 1}. 1 means long the spread, -1 means short, 0 means flat.

find_cointegrated_pairs(prices_df, significance=0.05)[source]¶

Scan a DataFrame of price series and find all cointegrated pairs.

For each pair of columns, the Engle-Granger cointegration test is applied. Pairs with a p-value below significance are returned.

Parameters:

prices_df (DataFrame) – DataFrame of price series (columns = asset names).
significance (float, default: 0.05) – Significance level for the cointegration test.

Return type:

list[tuple]

Returns:

List of tuples (asset1, asset2, p_value, hedge_ratio) for each cointegrated pair, sorted by p-value ascending.

ols(y, X, add_constant=True)[source]¶

Ordinary least squares regression.

OLS finds the linear coefficients that minimise the sum of squared residuals. It is the foundation of empirical finance – used for CAPM beta estimation, factor model fitting, and return attribution.

When to use:

Use OLS as the default regression when you have a cross-section or time-series of returns to explain. Switch to WLS if errors are heteroskedastic, or to Newey-West if errors are also autocorrelated (common in overlapping return regressions).

Mathematical formulation:

y = X * beta + epsilon beta_hat = (X’X)^{-1} X’y

How to interpret:

coefficients[0] is the intercept (alpha in CAPM).
coefficients[1:] are the factor loadings (betas).
t_stats and p_values: test H0: beta_i = 0. Reject if |t| > 2 (roughly, p < 0.05).
r_squared: fraction of variance explained. For CAPM on individual stocks, R^2 of 0.10-0.30 is typical.
residuals: unexplained returns. Check for autocorrelation, heteroskedasticity, and normality.

Parameters:

y (Series | ndarray) – Dependent variable (e.g., asset returns).
X (DataFrame | ndarray) – Independent variables (e.g., factor returns). Can be 1-D (single factor) or 2-D (multiple factors).
add_constant (bool, default: True) – Whether to add an intercept column to X.

Return type:

Returns:

Dictionary with coefficients (array), t_stats (array), p_values (array), r_squared, adj_r_squared, and residuals (array).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = np.random.normal(0.0005, 0.01, 252)
>>> stock = 0.001 + 1.2 * market + np.random.normal(0, 0.005, 252)
>>> result = ols(stock, market)
>>> abs(result["coefficients"][1] - 1.2) < 0.3  # beta near 1.2
True

See also

wls: Weighted least squares. newey_west_ols: OLS with HAC-robust standard errors. rolling_ols: Time-varying coefficient estimation.

rolling_ols(y, X, window=60, add_constant=True)[source]¶

Rolling window OLS regression.

Fits OLS independently at each time step using only the most recent window observations, producing time-varying coefficient estimates. This is essential for detecting parameter instability – a common feature of financial data where betas, hedge ratios, and factor loadings evolve over time.

When to use:

Use rolling OLS when: - You suspect the relationship between variables changes over

time (e.g., a stock’s beta increasing during crises).

You need time-varying hedge ratios for pairs trading.
You want to validate that a full-sample OLS result is stable.

For a more adaptive approach, consider Kalman regression (wraquant.regimes.kalman_regression), which estimates time-varying coefficients with a state-space model rather than fixed windows.

How to interpret:

coefficients: DataFrame of rolling betas. Plot to see how each coefficient evolves. Large swings indicate parameter instability.
r_squared: rolling R-squared. A declining R-squared suggests the model’s explanatory power is deteriorating.

Parameters:

y (Series) – Dependent variable series.
X (DataFrame | Series) – Independent variable(s). A Series is treated as a single regressor; a DataFrame may contain multiple regressors.
window (int, default: 60) – Rolling window size (e.g., 60 for ~3 months of daily data). Shorter windows are more responsive but noisier.
add_constant (bool, default: True) – Whether to add an intercept column.

Return type:

Returns:

Dictionary with coefficients (DataFrame of rolling betas, NaN before the first full window) and r_squared (Series of rolling R-squared values).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = pd.Series(np.random.normal(0, 0.01, 252))
>>> stock = 1.2 * market + np.random.normal(0, 0.005, 252)
>>> result = rolling_ols(stock, market, window=60)
>>> result["coefficients"].dropna().iloc[-1, 1]  # beta estimate
1.2...

See also

ols: Full-sample OLS regression. wraquant.regimes.kalman_regression: Kalman-filter-based

time-varying regression.

wls(y, X, weights, add_constant=True)[source]¶

Weighted least squares regression.

WLS accounts for heteroskedasticity by assigning different weights to observations. Observations with higher weight have more influence on the coefficient estimates. This is the appropriate estimator when the variance of the error term differs across observations.

When to use:

Use WLS when: - You know or suspect heteroskedasticity (non-constant error

variance). For example, high-cap stocks have less noisy returns than micro-caps.

You want to give more weight to recent observations (exponentially decaying weights for time-series regression).
You have grouped data where some groups are more precisely measured than others.

If you do not know the weights, use OLS with Newey-West (HAC) standard errors (newey_west_ols) instead.

Mathematical formulation:

Minimises sum_i w_i * (y_i - X_i * beta)^2

Equivalent to OLS on the transformed system: sqrt(w_i) * y_i = sqrt(w_i) * X_i * beta + epsilon_i

How to interpret:

Same output structure as ols. The difference is that coefficients are efficient under heteroskedasticity (lower standard errors than OLS when the weights are correctly specified). If the weights are misspecified, WLS can be worse than OLS.

Parameters:

y (Series | ndarray) – Dependent variable.
X (DataFrame | ndarray) – Independent variables.
weights (Series | ndarray) – Observation weights (higher weight = more influence). Common choices: inverse variance, exponential decay, or sample size per group.
add_constant (bool, default: True) – Whether to add an intercept column.

Return type:

Returns:

Dictionary with coefficients, t_stats, p_values, r_squared, adj_r_squared, and residuals.

Example

>>> import numpy as np
>>> np.random.seed(42)
>>> X = np.random.randn(100, 2)
>>> y = X @ [1.5, -0.5] + np.random.randn(100) * (1 + np.abs(X[:, 0]))
>>> w = 1.0 / (1 + np.abs(X[:, 0]))  # inverse heteroskedasticity
>>> result = wls(y, X, weights=w)
>>> len(result["coefficients"])
3

See also

ols: Unweighted OLS (assumes homoskedasticity). newey_west_ols: OLS with robust standard errors.

fama_macbeth(panel_y, panel_X)[source]¶

Fama-MacBeth two-pass cross-sectional regression.

The Fama-MacBeth procedure is the standard methodology for testing whether a factor commands a risk premium in cross-sectional asset pricing. It handles the errors-in-variables problem that arises when estimated betas are used as regressors.

When to use:

Use Fama-MacBeth when: - You want to estimate risk premia for factors (e.g., “does

size, value, or momentum have a positive risk premium?”).

You have panel data: many assets observed over many time periods.
You want t-statistics that properly account for cross- sectional correlation (unlike pooled OLS, which overstates significance).

For time-series factor models (e.g., CAPM alpha of a single fund), use plain OLS instead.

Mathematical formulation:

Pass 1: For each time period t, run cross-sectional OLS:: r_{i,t} = gamma_{0,t} + gamma_{1,t} * beta_{i,t} + e_{i,t}
Pass 2: Average the cross-sectional slopes over time:: gamma_k = mean(gamma_{k,t}) t_k = gamma_k / (std(gamma_{k,t}) / sqrt(T))

How to interpret:

risk_premia: average slope coefficients. A positive, statistically significant risk premium means the factor is priced.
t_stats: Fama-MacBeth t-statistics. |t| > 2 suggests the premium is statistically significant.
gamma_series: period-by-period slopes. Plot to see time variation in risk premia.
r_squared: average explanatory power of the factors in the cross-section.

Parameters:

panel_y (DataFrame) – DataFrame of returns with shape (T, N) where T is the number of time periods and N is the number of assets.
panel_X (DataFrame | dict[str, DataFrame]) – Factor exposures. Either a single DataFrame with the same shape as panel_y (single factor), or a dictionary mapping factor names to DataFrames of exposures.

Return type:

Returns:

Dictionary with risk_premia (array), t_stats (array), r_squared (float), and gamma_series (DataFrame of period-by-period slope coefficients).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> T, N = 120, 30
>>> betas = np.random.randn(T, N)
>>> returns = pd.DataFrame(0.005 * betas + np.random.randn(T, N) * 0.01)
>>> exposures = pd.DataFrame(betas)
>>> result = fama_macbeth(returns, exposures)
>>> len(result["risk_premia"])
2

See also

ols: Single time-series OLS regression. newey_west_ols: OLS with HAC-robust standard errors.

References

Fama & MacBeth (1973), “Risk, Return, and Equilibrium: Empirical Tests”

newey_west_ols(y, X, max_lags=None, add_constant=True)[source]¶

OLS regression with Newey-West HAC standard errors.

Produces the same coefficient estimates as OLS, but computes standard errors that are robust to both heteroskedasticity and autocorrelation (HAC). This is the standard approach in financial econometrics where both issues are almost always present.

When to use:

Use Newey-West whenever: - The regression involves overlapping returns (e.g., monthly

returns sampled from daily data) which introduce mechanical autocorrelation.

You suspect GARCH-type volatility clustering in the residuals (heteroskedasticity + autocorrelation).
You want valid inference without specifying the exact form of heteroskedasticity (unlike WLS which requires known weights).

Newey-West is strictly better than plain OLS for inference in financial data. The point estimates are the same; only the standard errors, t-statistics, and p-values differ.

Mathematical formulation:

The HAC variance estimator is:: V_HAC = (X’X)^{-1} S (X’X)^{-1}

where S = sum_{j=-L}^{L} w(j) * Gamma_j is the kernel- weighted sum of autocovariance matrices of X * epsilon, using the Bartlett kernel w(j) = 1 - |j|/(L+1).

How to interpret:

Same output as ols, plus hac_se (the HAC standard errors). Compare hac_se to the OLS standard errors: if HAC SE > OLS SE, the OLS was underestimating uncertainty (common in finance). Use the HAC-based t_stats and p_values for inference.

Parameters:

y (Series | ndarray) – Dependent variable.
X (DataFrame | ndarray) – Independent variables.
max_lags (int | None, default: None) – Maximum number of lags for the Newey-West kernel. When None, uses floor(4 * (T/100)^(2/9)) (Andrews 1991 rule of thumb).
add_constant (bool, default: True) – Whether to add an intercept column.

Return type:

Returns:

Dictionary with coefficients, t_stats, p_values, r_squared, adj_r_squared, residuals, and hac_se (HAC standard errors).

Example

>>> import numpy as np
>>> np.random.seed(42)
>>> X = np.random.randn(200, 2)
>>> y = X @ [1.0, 0.5] + np.random.randn(200) * 0.5
>>> result = newey_west_ols(y, X)
>>> len(result["hac_se"])
3

See also

ols: OLS with classical (non-robust) standard errors. wls: Weighted least squares for known heteroskedasticity.

References

Newey & West (1987), “A Simple, Positive Semi-Definite, Heteroskedasticity and Autocorrelation Consistent Covariance Matrix”

fama_french_regression(returns, factors_df)[source]¶

Regress asset returns on Fama-French factors.

The regression is R_i - R_f = alpha + beta_1 * F_1 + ... + eps. The factors DataFrame should contain the factor returns (e.g., Mkt-RF, SMB, HML, and optionally RMW, CMA). If a column named RF is present it is used to compute excess returns; otherwise returns are assumed to already be excess returns.

Parameters:

returns (Series) – Asset return series.
factors_df (DataFrame) – DataFrame of factor returns. Columns are factor names. An optional RF column is the risk-free rate.

Return type:

Returns:

Dictionary with alpha (intercept), betas (dict mapping factor name to coefficient), t_stats (dict mapping name to t-statistic), p_values (dict), and r_squared.

factor_attribution(returns, factor_returns)[source]¶

Decompose returns into factor contributions and specific return.

Runs a regression of returns on factor_returns and attributes the mean return to each factor.

Parameters:

returns (Series) – Asset return series.
factor_returns (DataFrame) – DataFrame of factor return series.

Return type:

Returns:

Dictionary with factor_contributions (dict mapping factor name to its average contribution), specific_return (mean residual return), total_return (mean of returns), and r_squared.

information_coefficient(predictions, returns)[source]¶

Compute the information coefficient (Spearman rank correlation).

The IC measures the predictive power of a signal: the rank correlation between the cross-sectional predictions and subsequent realised returns.

Parameters:

predictions (Series | ndarray) – Predicted values or signals.
returns (Series | ndarray) – Subsequent realised returns.

Return type:

Returns:

Spearman rank correlation coefficient (between -1 and 1).

quantile_analysis(predictions, returns, n_quantiles=5)[source]¶

Analyse returns by prediction quantile.

Sorts observations into quantiles based on predictions and computes summary statistics for each quantile bucket.

Parameters:

predictions (Series) – Predicted values or signals.
returns (Series) – Subsequent realised returns.
n_quantiles (int, default: 5) – Number of quantile buckets (default 5 = quintiles).

Return type:

Returns:

DataFrame indexed by quantile (1 = lowest, n_quantiles = highest) with columns mean_return, std_return, hit_rate (fraction of positive returns), and count.

pca_factors(returns, n_components=3, method='svd')[source]¶

Extract statistical factors from a returns matrix using PCA.

Parameters:

returns (DataFrame) – DataFrame of asset returns with shape (T, N) where T is the number of observations and N is the number of assets.
n_components (int, default: 3) – Number of principal components to retain.
method (str, default: 'svd') – Decomposition method. "svd" (default) uses numpy.linalg.svd on the demeaned returns; "eig" uses eigendecomposition of the covariance matrix.

Returns:

factors: DataFrame of extracted factors (T, n_components).
loadings: DataFrame of factor loadings (N, n_components).
explained_variance: array of variance explained by each component.
explained_variance_ratio: array of fraction of total variance explained.

Return type:

factor_loadings(returns, factors)[source]¶

Compute factor loadings by regressing each asset’s returns on factors.

For each asset column in returns, an OLS regression is run against factors (with intercept) to obtain betas (loadings), alphas, and R-squared values.

Parameters:

returns (DataFrame | ndarray) – Asset returns matrix (T, N).
factors (DataFrame | ndarray) – Factor returns matrix (T, K).

Returns:

loadings: ndarray (N, K) of factor betas.
alphas: ndarray (N,) of intercepts.
r_squared: ndarray (N,) of regression R-squared values.
residuals: ndarray (T, N) of regression residuals.

Return type:

scree_plot_data(returns)[source]¶

Return eigenvalues and explained variance ratios for a scree plot.

Parameters:

returns (DataFrame | ndarray) – Asset returns matrix (T, N).

Returns:

eigenvalues: array of eigenvalues in descending order.
explained_variance_ratio: array of fraction of total variance explained by each component.
cumulative_variance_ratio: cumulative sum of explained variance ratios.

Return type:

varimax_rotation(loadings, max_iter=500, tol=1e-06)[source]¶

Apply varimax rotation to a factor loadings matrix.

Varimax maximises the sum of variances of squared loadings within each factor, producing a simpler (more interpretable) structure.

Parameters:

loadings (ndarray) – Factor loadings matrix (N, K) where N is the number of variables and K is the number of factors.
max_iter (int, default: 500) – Maximum number of iterations.
tol (float, default: 1e-06) – Convergence tolerance on the rotation criterion change.

Returns:

rotated_loadings: ndarray (N, K) of rotated loadings.
rotation_matrix: ndarray (K, K) orthogonal rotation matrix.
n_iter: number of iterations performed.

Return type:

factor_mimicking_portfolios(returns, characteristics, n_quantiles=5)[source]¶

Build factor-mimicking portfolios via long-short quantile sorts.

For each characteristic column, assets are sorted into quantiles at each time step. The factor-mimicking return is the difference between the mean return of the top quantile and the mean return of the bottom quantile (long top, short bottom).

Parameters:

returns (DataFrame) – DataFrame of asset returns (T, N).
characteristics (DataFrame) – DataFrame of asset characteristics (T, N) or (N,) for a static sort. If a single row / Series, the same sort is applied to every period. If multiple columns, each column is treated as a separate characteristic, with returns as a single panel.
n_quantiles (int, default: 5) – Number of quantile buckets (default 5).

Return type:

Returns:

DataFrame of factor-mimicking portfolio returns (T, K) where K is the number of characteristics.

risk_factor_decomposition(portfolio_returns, factor_returns)[source]¶

Decompose portfolio risk into factor risk and idiosyncratic risk.

Runs an OLS regression of portfolio_returns on factor_returns and computes the variance attributable to each factor and the residual (idiosyncratic) variance.

Parameters:

portfolio_returns (Series | ndarray) – Portfolio return series (T,).
factor_returns (DataFrame | ndarray) – Factor return matrix (T, K).

Returns:

total_variance: total variance of portfolio returns.
factor_variance: variance explained by the factor model.
idiosyncratic_variance: residual variance.
factor_risk_share: fraction of total variance from factors.
idiosyncratic_risk_share: fraction of total variance from idiosyncratic risk.
betas: regression coefficients (loadings) on each factor.
factor_marginal_contributions: variance contribution of each factor.

Return type:

factor_correlation(factor_returns)[source]¶

Compute the correlation matrix of factors with significance tests.

For each pair of factors, the Pearson correlation and a two-sided p-value (testing H0: rho = 0) are computed.

Parameters:

factor_returns (DataFrame | ndarray) – Factor return matrix (T, K).

Returns:

correlation: ndarray (K, K) correlation matrix.
p_values: ndarray (K, K) of p-values for each pair.

Return type:

common_factors(returns_list, n_components=3)[source]¶

Find common factors shared across multiple asset classes.

Extracts PCA factors from each asset class independently, then performs canonical correlation analysis on the stacked factor scores to identify shared latent factors.

Parameters:

returns_list (list[DataFrame | ndarray]) – List of returns matrices, one per asset class. Each has shape (T, N_i) and all must share the same number of time observations T.
n_components (int, default: 3) – Number of PCA components to extract per asset class before cross-analysis.

Returns:

individual_factors: list of factor arrays, one per asset class.
common_factor_scores: ndarray (T, n_components) of shared factor scores obtained from a second-level PCA on concatenated individual factors.
cross_correlations: correlation matrix between the individual factor sets.

Return type:

fama_french_factors(returns, characteristics, n_quantiles=3)[source]¶

Construct Fama-French style factors from a cross-section of returns.

Sorts assets into portfolios based on each characteristic at each time period, then computes long-short (top-minus-bottom quantile) factor returns. This replicates the standard methodology used by Fama and French to construct SMB, HML, and related factors.

When to use:

To create custom factors from firm characteristics (e.g., book-to-market, momentum, profitability, investment).
To replicate or extend the Fama-French factor zoo.
To test whether a new characteristic is a priced risk factor (construct the factor, then test its risk premium via fama_macbeth).

Mathematical formulation:

For each period t and each characteristic c:

Sort assets into q quantiles based on the characteristic.
Compute the equal-weighted mean return for the top and bottom quantile.
Factor return = mean(top quantile returns) - mean(bottom quantile returns).

This is a zero-cost, long-short portfolio that isolates the return premium associated with the characteristic.

How to interpret:

A consistently positive factor return means that assets with high values of the characteristic outperform those with low values (and vice versa for negative).
The t-statistic of the mean factor return tests whether the premium is significantly different from zero.
Standard Fama-French uses terciles (3 groups) for the main sort and independent double sorts for intersections (e.g., size x value).

Parameters:

returns (DataFrame) – DataFrame of asset returns (T, N) with a DatetimeIndex and asset names as columns.
characteristics (DataFrame) – DataFrame of asset characteristics (T, N) with the same index and columns as returns, or a DataFrame with (N, K) where N assets are the index and K characteristics are the columns (static sort).
n_quantiles (int, default: 3) – Number of quantile buckets (default 3 = terciles, matching Fama-French convention).

Return type:

Returns:

DataFrame of factor returns (T, K) where K is the number of characteristics.

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> T, N = 100, 30
>>> dates = pd.bdate_range("2020-01-01", periods=T)
>>> assets = [f"s{i}" for i in range(N)]
>>> ret = pd.DataFrame(rng.normal(0, 0.02, (T, N)), index=dates, columns=assets)
>>> bm = pd.DataFrame(rng.uniform(0.5, 3.0, (T, N)), index=dates, columns=assets)
>>> ff = fama_french_factors(ret, bm)
>>> ff.shape[0] == T
True

See also

factor_mimicking_portfolios: General factor-mimicking portfolio: construction.

factor_exposure: Regress returns on constructed factors.

factor_exposure(returns, factor_returns)[source]¶

Regress returns on factor returns to estimate factor exposures (betas).

For each asset (or a single portfolio), runs an OLS regression of returns on factor returns (with intercept) and reports the factor betas, t-statistics, and R-squared.

When to use:

To estimate a portfolio’s or asset’s exposure to known factors (e.g., market, size, value, momentum).
As the first step in factor-based risk decomposition.
To validate that a factor-neutral strategy truly has zero exposure to target factors.

Mathematical formulation:

For each asset i:

\[r_i = \alpha_i + \sum_{k=1}^K \beta_{ik} f_k + \epsilon_i\]

The betas are the OLS coefficients on the factor returns.

How to interpret:

beta > 0: positive exposure to the factor (moves with it).
beta = 0: no exposure.
|t_stat| > 2: the exposure is statistically significant.
R_squared: fraction of return variance explained by the factor model. Higher R-squared means the model is a good fit.

Parameters:

returns (DataFrame | Series) – Asset returns. A Series for a single asset or a DataFrame (T, N) for multiple assets.
factor_returns (DataFrame) – DataFrame of factor returns (T, K).

Returns:

alpha, beta_<factor_name> for each factor, t_<factor_name> for each factor’s t-statistic, and r_squared.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> T = 200
>>> mkt = pd.Series(rng.normal(0, 0.01, T), name="MKT")
>>> smb = pd.Series(rng.normal(0, 0.005, T), name="SMB")
>>> factors = pd.DataFrame({"MKT": mkt, "SMB": smb})
>>> ret = 1.2 * mkt + 0.5 * smb + rng.normal(0, 0.003, T)
>>> result = factor_exposure(pd.Series(ret, name="fund"), factors)
>>> abs(result.loc["fund", "beta_MKT"] - 1.2) < 0.3
True

See also

factor_loadings: Lower-level loadings estimation. factor_risk_decomposition: Risk decomposition from exposures.

factor_risk_decomposition(returns, factor_returns)[source]¶

Decompose total risk into systematic (factor) and idiosyncratic components.

Runs a factor model regression, then separates the total variance of the return series into the portion explained by the factors (systematic risk) and the unexplained residual (idiosyncratic risk).

When to use:

To understand what fraction of a portfolio’s risk comes from common factor exposures vs. stock-specific bets.
For risk budgeting: allocate risk limits to systematic and idiosyncratic components.
To evaluate diversification: a well-diversified portfolio has low idiosyncratic risk relative to total risk.

Mathematical formulation:

\[\text{Var}(r) = \beta' \Sigma_f \beta + \sigma^2_\epsilon\]

where \beta is the vector of factor exposures, \Sigma_f is the factor covariance matrix, and \sigma^2_\epsilon is the idiosyncratic variance.

The R-squared of the regression gives the systematic risk share:

\[R^2 = 1 - \frac{\sigma^2_\epsilon}{\text{Var}(r)}\]

Parameters:

returns (Series) – Return series for a single asset or portfolio.
factor_returns (DataFrame) – DataFrame of factor returns (T, K).

Returns:

systematic_risk: variance explained by factors.
idiosyncratic_risk: residual variance.
total_risk: total return variance.
R_squared: fraction of risk that is systematic.
betas: factor exposure coefficients.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> mkt = rng.normal(0, 0.01, 200)
>>> ret = pd.Series(1.2 * mkt + rng.normal(0, 0.005, 200))
>>> factors = pd.DataFrame({"MKT": mkt})
>>> result = factor_risk_decomposition(ret, factors)
>>> result["R_squared"] > 0.3
True

See also

risk_factor_decomposition: Lower-level decomposition with: marginal contributions.

factor_exposure: Factor exposure estimation.

Descriptive Statistics¶

Descriptive statistics for financial return and price series.

summary_stats(returns)[source]¶

Compute summary statistics for a return series.

Parameters:: returns (Series) – Simple return series.
Return type:: dict
Returns:: Dictionary with mean, std, skew, kurtosis, min, max, and count.

annualized_return(returns, periods_per_year=252)[source]¶

Compute annualized return from a simple return series.

Parameters:

returns (Series) – Simple return series.
periods_per_year (int, default: 252) – Number of periods per year (252 for daily).

Return type:

Returns:

Annualized return as a float.

annualized_volatility(returns, periods_per_year=252)[source]¶

Compute annualized volatility from a simple return series.

Parameters:

returns (Series) – Simple return series.
periods_per_year (int, default: 252) – Number of periods per year (252 for daily).

Return type:

Returns:

Annualized volatility as a float.

max_drawdown(prices)[source]¶

Compute maximum drawdown from a price series.

Parameters:: prices (Series) – Price series (not returns).
Return type:: float
Returns:: Maximum drawdown as a negative float (e.g., -0.25 for 25% drawdown).

calmar_ratio(returns, periods_per_year=252)[source]¶

Compute the Calmar ratio (annualized return / max drawdown).

Parameters:

returns (Series) – Simple return series.
periods_per_year (int, default: 252) – Number of periods per year.

Return type:

Returns:

Calmar ratio as a float.

omega_ratio(returns, threshold=0.0)[source]¶

Compute the Omega ratio.

The Omega ratio is the probability-weighted ratio of gains versus losses relative to a threshold.

Parameters:

returns (Series) – Simple return series.
threshold (float, default: 0.0) – Return threshold (default 0).

Return type:

Returns:

Omega ratio as a float.

rolling_sharpe(returns, window=60, risk_free_rate=0.0, periods_per_year=252)[source]¶

Compute the rolling Sharpe ratio over a moving window.

When to use:

To monitor strategy performance stability over time.
To detect regime changes in risk-adjusted returns (e.g., a strategy that worked pre-2020 but degraded post-2020).
To compare two strategies’ time-varying risk-adjusted performance.

Mathematical formulation:

For each window of length w:

\[\text{Sharpe}_t = \frac{\bar{r}_t - r_f}{\sigma_t} \cdot \sqrt{P}\]

where \bar{r}_t and \sigma_t are the rolling mean and standard deviation of returns, r_f is the per-period risk-free rate, and P is the annualisation factor (e.g., 252 for daily).

How to interpret:

Sharpe > 1.0: good risk-adjusted performance (annualised).
Sharpe > 2.0: very strong.
Sharpe < 0.0: losing money on a risk-adjusted basis.
Large swings in rolling Sharpe indicate unstable performance.

Parameters:

returns (Series) – Simple return series.
window (int, default: 60) – Rolling window size in periods (default 60, roughly 3 months of daily data).
risk_free_rate (float, default: 0.0) – Per-period risk-free rate (default 0.0).
periods_per_year (int, default: 252) – Annualisation factor (252 for daily data).

Return type:

Returns:

Rolling Sharpe ratio as a pd.Series. First window - 1 values are NaN.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> ret = pd.Series(np.random.normal(0.001, 0.02, 252))
>>> rs = rolling_sharpe(ret, window=60)
>>> rs.dropna().shape[0]
193

See also

annualized_volatility: Annualised standard deviation. calmar_ratio: Drawdown-based risk-adjusted return.

rolling_drawdown(returns, window=60)[source]¶

Compute the rolling maximum drawdown over a moving window.

For each time step, the maximum drawdown is computed using only the most recent window observations. This provides a time-varying measure of downside risk.

When to use:

To monitor the worst-case loss over a recent period.
To detect periods of elevated tail risk that may not show up in rolling volatility.
As an input to risk overlays that tighten exposure when recent drawdowns are deep.
To compare the downside risk profile of different strategies over time.

Mathematical formulation:

For each window [t - w + 1, t], compute the cumulative return series from the window’s returns, find the peak, and measure the maximum drop from peak to trough:

\[\text{MDD}_t = \min_{s \in [t-w+1, t]} \frac{P_s - \max_{u \le s} P_u}{\max_{u \le s} P_u}\]

How to interpret:

Values are negative (or zero). More negative = deeper drawdown.
A rolling drawdown of -0.10 means the portfolio lost 10% from its peak within the window.
Compare to static max drawdown to see if the worst period is concentrated or distributed.

Parameters:

returns (Series) – Simple return series.
window (int, default: 60) – Rolling window size in periods (default 60).

Return type:

Returns:

Rolling maximum drawdown as a pd.Series. First window - 1 values are NaN.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> ret = pd.Series(np.random.normal(0, 0.02, 252))
>>> rd = rolling_drawdown(ret, window=60)
>>> (rd.dropna() <= 0).all()
True

See also

max_drawdown: Full-sample maximum drawdown. calmar_ratio: Return / drawdown ratio.

return_attribution(portfolio_weights, benchmark_weights, portfolio_returns, benchmark_returns)[source]¶

Decompose portfolio excess return using the Brinson-Fachler model.

When to use:

To explain why a portfolio outperformed or underperformed its benchmark.
To separate the contribution of top-down asset allocation decisions from bottom-up security selection.
For reporting to investors or risk committees.

Mathematical formulation:

For each asset/sector i:

Allocation effect: (w^P_i - w^B_i) * (r^B_i - r^B_total)
Selection effect: w^B_i * (r^P_i - r^B_i)
Interaction effect: (w^P_i - w^B_i) * (r^P_i - r^B_i)
Total active return: sum(allocation + selection + interaction) = r^P_total - r^B_total

How to interpret:

Positive allocation: the portfolio overweighted sectors that outperformed the benchmark.
Positive selection: within each sector, the portfolio held better-performing securities.
The interaction term captures the joint effect.
The three components sum to the total excess return.

Parameters:

portfolio_weights (Series) – Portfolio weights per asset/sector.
benchmark_weights (Series) – Benchmark weights per asset/sector (same index).
portfolio_returns (Series) – Portfolio returns per asset/sector.
benchmark_returns (Series) – Benchmark returns per asset/sector (same index).

Returns:

allocation: total allocation effect (float).
selection: total selection effect (float).
interaction: total interaction effect (float).
total_excess: total excess return (float).
detail: DataFrame with per-asset breakdown.

Return type:

Example

>>> import pandas as pd
>>> pw = pd.Series({"Tech": 0.4, "Fin": 0.3, "Health": 0.3})
>>> bw = pd.Series({"Tech": 0.3, "Fin": 0.4, "Health": 0.3})
>>> pr = pd.Series({"Tech": 0.05, "Fin": 0.02, "Health": 0.03})
>>> br = pd.Series({"Tech": 0.04, "Fin": 0.03, "Health": 0.03})
>>> result = return_attribution(pw, bw, pr, br)
>>> abs(result["total_excess"] - (pw @ pr - bw @ br)) < 1e-10
True

See also

risk_contribution: Per-asset risk decomposition.

risk_contribution(weights, cov_matrix)[source]¶

Compute marginal risk contributions per asset.

Risk contribution measures how much each asset contributes to the total portfolio risk (standard deviation). This is the foundation of risk-parity and risk-budgeting portfolio construction.

When to use:

To understand where portfolio risk comes from.
To build risk-parity portfolios where each asset contributes equally to total risk.
For risk monitoring: detect when a single position dominates portfolio risk.
To compare intended risk budgets with realised risk allocation.

Mathematical formulation:

The marginal contribution to risk (MCR) for asset i is:

\[\text{MCR}_i = w_i \cdot \frac{(\Sigma w)_i}{\sigma_p}\]

where \Sigma is the covariance matrix, w is the weight vector, and \sigma_p = \sqrt{w' \Sigma w} is the portfolio standard deviation.

The marginal contributions sum to the total portfolio risk:

\[\sum_i \text{MCR}_i = \sigma_p\]

How to interpret:

Values are in the same units as portfolio standard deviation.
Each value represents the portion of total portfolio risk attributable to that asset.
Negative contributions are possible for assets that hedge overall portfolio risk.

Parameters:

weights (Series | ndarray) – Portfolio weights (1-D array or Series).
cov_matrix (DataFrame | ndarray) – Covariance matrix (2-D array or DataFrame).

Return type:

Returns:

pd.Series of marginal risk contributions per asset.

Example

>>> import pandas as pd, numpy as np
>>> w = pd.Series({"A": 0.5, "B": 0.3, "C": 0.2})
>>> cov = pd.DataFrame(
...     np.diag([0.04, 0.09, 0.01]),
...     index=["A", "B", "C"], columns=["A", "B", "C"],
... )
>>> rc = risk_contribution(w, cov)
>>> abs(rc.sum() - np.sqrt(w.values @ cov.values @ w.values)) < 1e-10
True

See also

return_attribution: Return decomposition (Brinson model). shrunk_covariance: Better covariance input for risk contributions.

Regression¶

Regression models for financial econometrics.

Provides OLS, WLS, rolling OLS, Fama-MacBeth cross-sectional regression, and Newey-West HAC-robust regression – the standard toolkit for empirical asset pricing, factor modelling, and return attribution.

How to choose:

OLS (ols): the starting point. Estimates the linear relationship between a dependent variable and regressors. Assumes homoskedastic, serially uncorrelated errors.
WLS (wls): use when observation reliability varies (e.g., weight by inverse variance, or give more weight to recent data). Common for heteroskedastic financial returns.
Newey-West OLS (newey_west_ols): use when residuals are both heteroskedastic and autocorrelated. Standard errors are HAC-robust, so t-statistics and p-values are reliable even when the OLS error assumptions fail (which they usually do in finance).
Rolling OLS (rolling_ols): use to track time-varying coefficients (e.g., evolving beta, hedge ratio). Essential for detecting parameter instability.
Fama-MacBeth (fama_macbeth): the standard for estimating risk premia in cross-sectional asset pricing. Two-pass procedure that handles the errors-in-variables problem.

References

Fama & MacBeth (1973), “Risk, Return, and Equilibrium: Empirical Tests”
Newey & West (1987), “A Simple, Positive Semi-Definite, Heteroskedasticity and Autocorrelation Consistent Covariance Matrix”

ols(y, X, add_constant=True)[source]¶

Ordinary least squares regression.

When to use:

Mathematical formulation:

y = X * beta + epsilon beta_hat = (X’X)^{-1} X’y

How to interpret:

coefficients[0] is the intercept (alpha in CAPM).
coefficients[1:] are the factor loadings (betas).
t_stats and p_values: test H0: beta_i = 0. Reject if |t| > 2 (roughly, p < 0.05).
r_squared: fraction of variance explained. For CAPM on individual stocks, R^2 of 0.10-0.30 is typical.
residuals: unexplained returns. Check for autocorrelation, heteroskedasticity, and normality.

Parameters:

y (Series | ndarray) – Dependent variable (e.g., asset returns).
X (DataFrame | ndarray) – Independent variables (e.g., factor returns). Can be 1-D (single factor) or 2-D (multiple factors).
add_constant (bool, default: True) – Whether to add an intercept column to X.

Return type:

Returns:

Dictionary with coefficients (array), t_stats (array), p_values (array), r_squared, adj_r_squared, and residuals (array).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = np.random.normal(0.0005, 0.01, 252)
>>> stock = 0.001 + 1.2 * market + np.random.normal(0, 0.005, 252)
>>> result = ols(stock, market)
>>> abs(result["coefficients"][1] - 1.2) < 0.3  # beta near 1.2
True

See also

wls: Weighted least squares. newey_west_ols: OLS with HAC-robust standard errors. rolling_ols: Time-varying coefficient estimation.

rolling_ols(y, X, window=60, add_constant=True)[source]¶

Rolling window OLS regression.

When to use:

Use rolling OLS when: - You suspect the relationship between variables changes over

time (e.g., a stock’s beta increasing during crises).

You need time-varying hedge ratios for pairs trading.
You want to validate that a full-sample OLS result is stable.

For a more adaptive approach, consider Kalman regression (wraquant.regimes.kalman_regression), which estimates time-varying coefficients with a state-space model rather than fixed windows.

How to interpret:

coefficients: DataFrame of rolling betas. Plot to see how each coefficient evolves. Large swings indicate parameter instability.
r_squared: rolling R-squared. A declining R-squared suggests the model’s explanatory power is deteriorating.

Parameters:

y (Series) – Dependent variable series.
X (DataFrame | Series) – Independent variable(s). A Series is treated as a single regressor; a DataFrame may contain multiple regressors.
window (int, default: 60) – Rolling window size (e.g., 60 for ~3 months of daily data). Shorter windows are more responsive but noisier.
add_constant (bool, default: True) – Whether to add an intercept column.

Return type:

Returns:

Dictionary with coefficients (DataFrame of rolling betas, NaN before the first full window) and r_squared (Series of rolling R-squared values).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = pd.Series(np.random.normal(0, 0.01, 252))
>>> stock = 1.2 * market + np.random.normal(0, 0.005, 252)
>>> result = rolling_ols(stock, market, window=60)
>>> result["coefficients"].dropna().iloc[-1, 1]  # beta estimate
1.2...

See also

ols: Full-sample OLS regression. wraquant.regimes.kalman_regression: Kalman-filter-based

time-varying regression.

wls(y, X, weights, add_constant=True)[source]¶

Weighted least squares regression.

When to use:

Use WLS when: - You know or suspect heteroskedasticity (non-constant error

variance). For example, high-cap stocks have less noisy returns than micro-caps.

You want to give more weight to recent observations (exponentially decaying weights for time-series regression).
You have grouped data where some groups are more precisely measured than others.

If you do not know the weights, use OLS with Newey-West (HAC) standard errors (newey_west_ols) instead.

Mathematical formulation:

Minimises sum_i w_i * (y_i - X_i * beta)^2

Equivalent to OLS on the transformed system: sqrt(w_i) * y_i = sqrt(w_i) * X_i * beta + epsilon_i

How to interpret:

Parameters:

y (Series | ndarray) – Dependent variable.
X (DataFrame | ndarray) – Independent variables.
weights (Series | ndarray) – Observation weights (higher weight = more influence). Common choices: inverse variance, exponential decay, or sample size per group.
add_constant (bool, default: True) – Whether to add an intercept column.

Return type:

Returns:

Dictionary with coefficients, t_stats, p_values, r_squared, adj_r_squared, and residuals.

Example

>>> import numpy as np
>>> np.random.seed(42)
>>> X = np.random.randn(100, 2)
>>> y = X @ [1.5, -0.5] + np.random.randn(100) * (1 + np.abs(X[:, 0]))
>>> w = 1.0 / (1 + np.abs(X[:, 0]))  # inverse heteroskedasticity
>>> result = wls(y, X, weights=w)
>>> len(result["coefficients"])
3

See also

ols: Unweighted OLS (assumes homoskedasticity). newey_west_ols: OLS with robust standard errors.

fama_macbeth(panel_y, panel_X)[source]¶

Fama-MacBeth two-pass cross-sectional regression.

When to use:

Use Fama-MacBeth when: - You want to estimate risk premia for factors (e.g., “does

size, value, or momentum have a positive risk premium?”).

You have panel data: many assets observed over many time periods.
You want t-statistics that properly account for cross- sectional correlation (unlike pooled OLS, which overstates significance).

For time-series factor models (e.g., CAPM alpha of a single fund), use plain OLS instead.

Mathematical formulation:

Pass 1: For each time period t, run cross-sectional OLS:: r_{i,t} = gamma_{0,t} + gamma_{1,t} * beta_{i,t} + e_{i,t}
Pass 2: Average the cross-sectional slopes over time:: gamma_k = mean(gamma_{k,t}) t_k = gamma_k / (std(gamma_{k,t}) / sqrt(T))

How to interpret:

risk_premia: average slope coefficients. A positive, statistically significant risk premium means the factor is priced.
t_stats: Fama-MacBeth t-statistics. |t| > 2 suggests the premium is statistically significant.
gamma_series: period-by-period slopes. Plot to see time variation in risk premia.
r_squared: average explanatory power of the factors in the cross-section.

Parameters:

panel_y (DataFrame) – DataFrame of returns with shape (T, N) where T is the number of time periods and N is the number of assets.
panel_X (DataFrame | dict[str, DataFrame]) – Factor exposures. Either a single DataFrame with the same shape as panel_y (single factor), or a dictionary mapping factor names to DataFrames of exposures.

Return type:

Returns:

Dictionary with risk_premia (array), t_stats (array), r_squared (float), and gamma_series (DataFrame of period-by-period slope coefficients).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> T, N = 120, 30
>>> betas = np.random.randn(T, N)
>>> returns = pd.DataFrame(0.005 * betas + np.random.randn(T, N) * 0.01)
>>> exposures = pd.DataFrame(betas)
>>> result = fama_macbeth(returns, exposures)
>>> len(result["risk_premia"])
2

See also

ols: Single time-series OLS regression. newey_west_ols: OLS with HAC-robust standard errors.

References

Fama & MacBeth (1973), “Risk, Return, and Equilibrium: Empirical Tests”

newey_west_ols(y, X, max_lags=None, add_constant=True)[source]¶

OLS regression with Newey-West HAC standard errors.

When to use:

Use Newey-West whenever: - The regression involves overlapping returns (e.g., monthly

returns sampled from daily data) which introduce mechanical autocorrelation.

You suspect GARCH-type volatility clustering in the residuals (heteroskedasticity + autocorrelation).
You want valid inference without specifying the exact form of heteroskedasticity (unlike WLS which requires known weights).

Newey-West is strictly better than plain OLS for inference in financial data. The point estimates are the same; only the standard errors, t-statistics, and p-values differ.

Mathematical formulation:

The HAC variance estimator is:: V_HAC = (X’X)^{-1} S (X’X)^{-1}

where S = sum_{j=-L}^{L} w(j) * Gamma_j is the kernel- weighted sum of autocovariance matrices of X * epsilon, using the Bartlett kernel w(j) = 1 - |j|/(L+1).

How to interpret:

Parameters:

y (Series | ndarray) – Dependent variable.
X (DataFrame | ndarray) – Independent variables.
max_lags (int | None, default: None) – Maximum number of lags for the Newey-West kernel. When None, uses floor(4 * (T/100)^(2/9)) (Andrews 1991 rule of thumb).
add_constant (bool, default: True) – Whether to add an intercept column.

Return type:

Returns:

Dictionary with coefficients, t_stats, p_values, r_squared, adj_r_squared, residuals, and hac_se (HAC standard errors).

Example

>>> import numpy as np
>>> np.random.seed(42)
>>> X = np.random.randn(200, 2)
>>> y = X @ [1.0, 0.5] + np.random.randn(200) * 0.5
>>> result = newey_west_ols(y, X)
>>> len(result["hac_se"])
3

See also

ols: OLS with classical (non-robust) standard errors. wls: Weighted least squares for known heteroskedasticity.

References

Newey & West (1987), “A Simple, Positive Semi-Definite, Heteroskedasticity and Autocorrelation Consistent Covariance Matrix”

Correlation¶

Correlation and covariance estimation for financial data.

Accurate correlation and covariance estimation is fundamental to portfolio construction, risk management, and factor modelling. The sample covariance matrix is a poor estimator when the number of assets (p) is comparable to or exceeds the number of observations (T) – a common situation in finance. This module provides shrinkage estimators that regularise the covariance matrix for more stable and better- conditioned estimates.

Key concepts:

Shrinkage blends the noisy sample covariance with a structured target (identity, diagonal, or constant correlation) to reduce estimation error.
The optimal shrinkage intensity balances bias (too much shrinkage) against variance (too little shrinkage).
All methods here produce positive semi-definite matrices, which is required for downstream use in optimisation.

References

Ledoit & Wolf (2004), “A well-conditioned estimator for large- dimensional covariance matrices”
Chen, Wiesel, Eldar & Hero (2010), “Shrinkage Algorithms for MMSE Covariance Estimation” (OAS)

correlation_matrix(returns, method='pearson')[source]¶

Compute a correlation matrix from asset returns.

When to use:

"pearson" (default): standard linear correlation. Use for most portfolio and factor analyses.
"spearman": rank correlation. Use when you suspect non-linear but monotonic relationships, or when returns have heavy tails / outliers.
"kendall": concordance-based. More robust than Spearman for small samples. Also connects naturally to copula models (Kendall’s tau has a direct relationship to copula parameters).

Parameters:

returns (DataFrame) – DataFrame of asset returns (columns = assets).
method (str, default: 'pearson') – Correlation method – "pearson", "spearman", or "kendall".

Return type:

Returns:

Correlation matrix as a DataFrame (p x p, symmetric, diagonal entries = 1.0).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame(np.random.randn(100, 3), columns=["A", "B", "C"])
>>> corr = correlation_matrix(returns)
>>> corr.shape
(3, 3)

See also

shrunk_covariance: Regularised covariance estimation. rolling_correlation: Time-varying pairwise correlation. wraquant.risk.copulas.rank_correlation: Kendall/Spearman for

copula analysis.

shrunk_covariance(returns, method='ledoit_wolf')[source]¶

Compute a shrinkage-estimated covariance matrix.

When to use:

Always prefer shrinkage over the raw sample covariance for portfolio optimisation. The improvement is largest when p/T is close to or exceeds 1 (e.g., 500 stocks with 252 daily observations).

"ledoit_wolf" (default): analytically optimal shrinkage toward a structured target. Best general-purpose choice. Automatically determines the optimal shrinkage intensity.
"oas" (Oracle Approximating Shrinkage): assumes the underlying distribution is Gaussian and computes the oracle- approximating shrinkage intensity. Slightly better than Ledoit-Wolf when normality holds.
"basic": simple shrinkage toward the diagonal with a fixed (non-optimal) shrinkage coefficient. Use only as a baseline.

Mathematical formulation:

Sigma_shrunk = (1 - alpha) * S + alpha * F

where S is the sample covariance, F is the shrinkage target (e.g., identity or diagonal), and alpha is the shrinkage intensity (0 = no shrinkage, 1 = full shrinkage to target).

How to interpret:

Parameters:

returns (DataFrame) – DataFrame of asset returns (columns = assets).
method (str, default: 'ledoit_wolf') – Shrinkage method – "ledoit_wolf" (default), "oas", or "basic".

Return type:

Returns:

Shrunk covariance matrix as a DataFrame (p x p).

Raises:

ValueError – If method is not recognized.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame(np.random.randn(100, 5), columns=list("ABCDE"))
>>> cov = shrunk_covariance(returns, method="ledoit_wolf")
>>> cov.shape
(5, 5)

See also

correlation_matrix: Correlation (standardised covariance). wraquant.ml.preprocessing.denoised_correlation: Random matrix

theory-based denoising.

References

Ledoit & Wolf (2004), “A well-conditioned estimator for large-dimensional covariance matrices”
Chen et al. (2010), “Shrinkage Algorithms for MMSE Covariance Estimation”

rolling_correlation(x, y, window)[source]¶

Compute rolling Pearson correlation between two series.

When to use:

Use rolling correlation to: - Monitor diversification benefit over time (correlation

rising toward 1.0 means diversification is eroding).

Detect correlation regime changes for pairs trading or hedging ratio adjustment.
Validate the stationarity assumption of portfolio optimisation inputs.

How to interpret:

Values near +1.0: strong positive co-movement (little diversification benefit).
Values near 0.0: approximately uncorrelated.
Values near -1.0: strong negative co-movement (excellent diversification or natural hedge).
Spikes toward +1.0 during sell-offs are typical (“correlation goes to 1 in a crisis”).

Parameters:

x (Series) – First return series.
y (Series) – Second return series (same index).
window (int) – Rolling window size (e.g., 60 for ~3 months of daily data).

Return type:

Returns:

Rolling Pearson correlation series. First window - 1 values are NaN.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> x = pd.Series(np.random.randn(200))
>>> y = pd.Series(0.5 * x + np.random.randn(200) * 0.5)
>>> rc = rolling_correlation(x, y, window=60)
>>> rc.dropna().iloc[0] > 0
True

See also

correlation_matrix: Full cross-asset correlation matrix. wraquant.risk.dcc.rolling_correlation_dcc: DCC-GARCH-based

dynamic correlation.

partial_correlation(data)[source]¶

Compute the partial correlation matrix, controlling for all other variables.

When to use:

Mathematical formulation:

For each pair (i, j), regress both X_i and X_j on all remaining variables, then compute the Pearson correlation of the residuals:

\[\rho_{ij \cdot \text{rest}} = \text{corr}(e_i, e_j)\]

where e_i = X_i - \hat{X}_i is the residual from regressing X_i on all other columns.

Equivalently, partial correlations can be obtained from the inverse of the correlation matrix (the precision matrix):

\[\rho_{ij \cdot \text{rest}} = -\frac{P_{ij}}{\sqrt{P_{ii} P_{jj}}}\]

where P = R^{-1} is the precision matrix.

How to interpret:

Values near 0 indicate no direct linear relationship once shared drivers are removed.
A large drop from raw correlation to partial correlation signals that the association is mostly indirect (mediated).
The diagonal is always 1.0.

Parameters:: data (DataFrame) – DataFrame with columns as variables (assets, factors, etc.) and rows as observations. Must have at least 3 columns.
Return type:: DataFrame
Returns:: Partial correlation matrix as a DataFrame (p x p, symmetric, diagonal = 1.0).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> z = np.random.randn(200)
>>> data = pd.DataFrame({
...     "A": z + np.random.randn(200) * 0.3,
...     "B": z + np.random.randn(200) * 0.3,
...     "C": np.random.randn(200),
... })
>>> pcorr = partial_correlation(data)
>>> pcorr.shape
(3, 3)

See also

correlation_matrix: Standard (marginal) correlation matrix. mutual_information: Non-linear dependence measure.

distance_correlation(x, y)[source]¶

Compute the Brownian distance correlation between two variables.

When to use:

You suspect a nonlinear relationship that Pearson/Spearman will miss (e.g., a U-shaped or threshold relationship).
You want a single-number summary of any type of dependence.
You need a test statistic for independence that is consistent against all alternatives with finite first moments.

Mathematical formulation:

Compute the pairwise Euclidean distance matrices a_{kl} = |X_k - X_l| and b_{kl} = |Y_k - Y_l|.
Double-center each matrix: A_{kl} = a_{kl} - \bar{a}_{k\cdot} - \bar{a}_{\cdot l} + \bar{a}_{\cdot\cdot}
Distance covariance squared: \text{dCov}^2(X, Y) = \frac{1}{n^2} \sum_{k,l} A_{kl} B_{kl}
Distance correlation: \text{dCor}(X, Y) = \frac{\text{dCov}(X, Y)}{\sqrt{\text{dVar}(X) \cdot \text{dVar}(Y)}}

How to interpret:

0.0: independence (no dependence of any kind).
1.0: perfect dependence (deterministic relationship).
Values between 0 and 1 indicate partial dependence.
Distance correlation >= |Pearson correlation|, so it always detects at least as much dependence.

Parameters:

x (Series | ndarray) – First variable (1-D array or Series).
y (Series | ndarray) – Second variable (1-D array or Series, same length).

Return type:

Returns:

Distance correlation as a float in [0, 1].

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 200)
>>> y = x ** 2 + rng.normal(0, 0.3, 200)  # nonlinear
>>> dcor = distance_correlation(x, y)
>>> dcor > 0.3  # detects nonlinear dependence
True

References

Szekely, G. J., Rizzo, M. L. & Bakirov, N. K. (2007). “Measuring and testing dependence by correlation of distances.” Annals of Statistics, 35(6), 2769-2794.

See also

correlation_matrix: Linear (Pearson) correlation. mutual_information: Information-theoretic dependence measure.

kendall_tau(x, y)[source]¶

Compute Kendall’s tau-b rank correlation coefficient with p-value.

When to use:

When data has heavy tails or outliers that distort Pearson correlation.
For small samples where Spearman is less reliable.
When you need a rank-based measure that connects naturally to copula parameters (Kendall’s tau has a one-to-one mapping to the parameter of many copula families).

Mathematical formulation:

\[\tau_b = \frac{C - D}{\sqrt{(C + D + T_x)(C + D + T_y)}}\]

where C = concordant pairs, D = discordant pairs, T_x and T_y = pairs tied only on x or y.

How to interpret:

+1: perfect concordance (monotonically increasing relationship).
-1: perfect discordance (monotonically decreasing).
0: no ordinal association.
|tau| > 0.3 is generally considered a moderate association in financial data.
The p-value tests H0: tau = 0 (independence).

Parameters:

x (Series | ndarray) – First variable (1-D array or Series).
y (Series | ndarray) – Second variable (1-D array or Series, same length).

Returns:

tau: Kendall’s tau-b statistic.
p_value: two-sided p-value for H0: tau = 0.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 100)
>>> y = 0.8 * x + rng.normal(0, 0.5, 100)
>>> result = kendall_tau(x, y)
>>> result["tau"] > 0
True
>>> result["p_value"] < 0.05
True

See also

correlation_matrix: Pearson/Spearman/Kendall full matrix. distance_correlation: Nonlinear dependence measure.

mutual_information(x, y, n_bins=20, method='binning')[source]¶

Estimate mutual information between two continuous variables.

When to use:

Feature selection for ML-based trading models: MI identifies features with any predictive signal, not just linear ones.
Comparing the information content of different alpha signals.
Measuring the ``true’’ dependence between returns and macro indicators that may have complex, non-monotonic relationships.

Mathematical formulation:

\[I(X; Y) = \sum_{x} \sum_{y} p(x, y) \log \frac{p(x, y)}{p(x) p(y)}\]

How to interpret:

MI = 0: independence (knowing X tells you nothing about Y).
MI > 0: some dependence exists.
MI is measured in nats (when using natural log) and is non-negative.
There is no upper bound in general, but normalised MI (MI / sqrt(H(X)*H(Y))) can be used for comparisons.

Parameters:

x (Series | ndarray) – First continuous variable.
y (Series | ndarray) – Second continuous variable (same length).
n_bins (int, default: 20) – Number of bins for discretisation ("binning" method). More bins capture finer structure but need more data.
method (str, default: 'binning') – Estimation method – "binning" (default) or "kde".

Return type:

Returns:

Estimated mutual information in nats (>= 0).

Raises:

ValueError – If method is not recognized.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 500)
>>> y = x + rng.normal(0, 0.5, 500)
>>> mi = mutual_information(x, y)
>>> mi > 0
True

See also

distance_correlation: Another non-linear dependence measure. correlation_matrix: Linear dependence only.

correlation_significance(x, y, method='pearson', confidence=0.95)[source]¶

Test whether the correlation between two variables is significantly non-zero.

Computes the sample correlation, performs a t-test for the null hypothesis H0: rho = 0, and constructs a confidence interval using Fisher’s z-transformation.

When to use:

To confirm that an observed correlation is statistically significant and not just sampling noise.
To obtain confidence intervals for reporting correlation estimates with uncertainty.
Before using a correlation estimate in portfolio construction or risk models — insignificant correlations may be unreliable.

Mathematical formulation:

Test statistic:

\[t = r \sqrt{\frac{n - 2}{1 - r^2}}\]

which follows a t-distribution with n - 2 degrees of freedom under H0.

Confidence interval via Fisher z-transformation:

\[z = \text{arctanh}(r), \quad SE = \frac{1}{\sqrt{n - 3}}\]

The interval [z - z_{\alpha/2} \cdot SE, z + z_{\alpha/2} \cdot SE] is back-transformed via tanh() to the correlation scale.

Parameters:

x (Series | ndarray) – First variable.
y (Series | ndarray) – Second variable (same length).
method (str, default: 'pearson') – Correlation method – "pearson" (default) or "spearman".
confidence (float, default: 0.95) – Confidence level for the interval (default 0.95).

Returns:

r: sample correlation coefficient.
t_stat: t-test statistic.
p_value: two-sided p-value for H0: rho = 0.
ci_lower: lower bound of the confidence interval.
ci_upper: upper bound of the confidence interval.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 100)
>>> y = 0.5 * x + rng.normal(0, 1, 100)
>>> result = correlation_significance(x, y)
>>> result["p_value"] < 0.05
True
>>> result["ci_lower"] < result["r"] < result["ci_upper"]
True

See also

correlation_matrix: Compute correlations without significance test. kendall_tau: Rank correlation with p-value.

minimum_spanning_tree_correlation(corr_matrix)[source]¶

Compute the minimum spanning tree (MST) of a correlation matrix.

When to use:

To visualise market structure and identify clusters of related assets (sectors, factor groups, regimes).
As input to hierarchical risk parity (HRP) portfolio construction (Lopez de Prado, 2016).
To detect changes in market structure over time by comparing MSTs across different periods.
To reduce dimensionality of the correlation matrix for network-based analysis.

Mathematical formulation:

The correlation matrix is converted to a distance matrix:

\[d_{ij} = \sqrt{2(1 - \rho_{ij})}\]

This metric satisfies the triangle inequality and maps perfect correlation (rho = 1) to zero distance and zero correlation (rho = 0) to distance sqrt(2).

Prim’s or Kruskal’s algorithm is then applied to find the MST of the complete weighted graph.

How to interpret:

Parameters:: corr_matrix (DataFrame) – Correlation matrix as a DataFrame (p x p, symmetric).
Return type:: DataFrame
Returns:: Adjacency matrix of the MST as a DataFrame (p x p). Non-zero entries indicate edges in the tree, with values equal to the distance sqrt(2 * (1 - rho)).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame(np.random.randn(100, 4), columns=list("ABCD"))
>>> corr = returns.corr()
>>> mst = minimum_spanning_tree_correlation(corr)
>>> mst.shape
(4, 4)
>>> (mst.values > 0).sum()  # MST has p-1 edges, each counted twice
6

References

Mantegna, R. N. (1999). “Hierarchical structure in financial markets.”
Lopez de Prado, M. (2016). “Building diversified portfolios that outperform out-of-sample.”

See also

correlation_matrix: Compute the input correlation matrix. shrunk_covariance: Regularised covariance for more stable MSTs.

Distributions¶

Distribution fitting and tail analysis for financial data.

fit_distribution(data, dist='norm')[source]¶

Fit a parametric distribution to data.

Parameters:

data (Series) – Data series to fit.
dist (str, default: 'norm') – Name of a scipy.stats distribution (e.g., "norm", "t", "lognorm").

Return type:

Returns:

Dictionary with params (tuple of fitted parameters), ks_statistic, and ks_pvalue from a Kolmogorov-Smirnov goodness-of-fit test.

Raises:

AttributeError – If dist is not a valid scipy distribution.

tail_ratio(returns, quantile=0.05)[source]¶

Compute the tail ratio (right tail / left tail).

A tail ratio > 1 indicates a fatter right tail (more extreme gains) relative to the left tail.

Parameters:

returns (Series) – Return series.
quantile (float, default: 0.05) – Quantile for tail measurement (default 5%).

Return type:

Returns:

Tail ratio as a float.

hurst_exponent(data)[source]¶

Estimate the Hurst exponent via rescaled range (R/S) analysis.

The Hurst exponent characterises the long-term memory of a series:

H < 0.5: mean-reverting
H = 0.5: random walk
H > 0.5: trending / persistent

Parameters:: data (Series) – Time series (prices or returns).
Return type:: float
Returns:: Estimated Hurst exponent as a float.

fit_stable_distribution(data)[source]¶

Fit a stable (Levy) distribution to data.

Uses scipy’s levy_stable distribution to estimate the four parameters: alpha (stability), beta (skewness), loc, and scale.

Parameters:

data (Series | ndarray) – Data array or series.

Returns:

alpha: stability parameter (0, 2].
beta: skewness parameter [-1, 1].
loc: location parameter.
scale: scale parameter.
ks_statistic: Kolmogorov-Smirnov statistic.
ks_pvalue: KS test p-value.

Return type:

tail_index(data, method='hill', threshold_quantile=0.9)[source]¶

Estimate the tail index of a distribution.

The tail index characterises the heaviness of distribution tails. A finite tail index indicates power-law tails (Pareto-like).

Parameters:

data (Series | ndarray) – Data array or series.
method (str, default: 'hill') – Estimation method. One of "hill" (Hill estimator), "pickands" (Pickands estimator), or "moment" (moment estimator of Dekkers-Einmahl-de Haan).
threshold_quantile (float, default: 0.9) – Quantile above which tail observations are used (default 0.9, i.e. top 10%).

Returns:

tail_index: estimated tail index (xi).
method: method used.
n_tail: number of observations in the tail.

Return type:

Raises:

ValueError – If method is not one of the supported estimators.

qqplot_data(data, dist='norm')[source]¶

Generate quantile-quantile plot data.

Computes theoretical and sample quantiles for constructing a Q-Q plot against a reference distribution.

Parameters:

data (Series | ndarray) – Data array or series.
dist (str, default: 'norm') – Name of a scipy.stats distribution to use as the theoretical reference (default "norm").

Returns:

theoretical_quantiles: array of theoretical quantiles.
sample_quantiles: array of ordered sample values.
slope: slope of the best-fit line through the Q-Q plot.
intercept: intercept of the best-fit line.

Return type:

jarque_bera(data)[source]¶

Perform the Jarque-Bera test for normality.

Tests the null hypothesis that the data is normally distributed, based on sample skewness and kurtosis.

Parameters:

data (Series | ndarray) – Data array or series.

Returns:

statistic: Jarque-Bera test statistic.
p_value: p-value of the test.
skewness: sample skewness.
kurtosis: sample excess kurtosis.

Return type:

kolmogorov_smirnov(data, dist='norm')[source]¶

Perform the Kolmogorov-Smirnov goodness-of-fit test.

Tests the null hypothesis that data was drawn from the specified distribution. The distribution parameters are first estimated via MLE.

Parameters:

data (Series | ndarray) – Data array or series.
dist (str, default: 'norm') – Name of a scipy.stats distribution (default "norm").

Returns:

statistic: KS test statistic.
p_value: p-value of the test.
dist: distribution name tested.
params: fitted distribution parameters.

Return type:

anderson_darling(data, dist='norm')[source]¶

Perform the Anderson-Darling goodness-of-fit test.

Parameters:

data (Series | ndarray) – Data array or series.
dist (str, default: 'norm') – Distribution to test against. Supported values depend on scipy.stats.anderson and include "norm", "expon", "logistic", "gumbel", "gumbel_l", "gumbel_r".

Returns:

statistic: Anderson-Darling test statistic.
critical_values: array of critical values for each significance level.
significance_levels: corresponding significance levels (%).

Return type:

Example

>>> import numpy as np
>>> data = np.random.default_rng(42).normal(0, 1, 1000)
>>> anderson_darling(data)

kernel_density_estimate(data, n_points=200, bandwidth='scott')[source]¶

Kernel density estimation using Gaussian kernels.

Parameters:

data (ndarray | Series) – Sample data (1D array of returns or prices).
n_points (int, default: 200) – Number of evaluation points for the density curve.
bandwidth (str | float, default: 'scott') – Bandwidth method (“scott”, “silverman”) or float.

Returns:

x – Evaluation points (n_points,).
density – Estimated density values (n_points,).
bandwidth – Bandwidth used.
mode – Location of peak density.
cdf – Cumulative distribution values (n_points,).

Return type:

dict[str, ndarray]

Example

>>> from wraquant.stats.distributions import kernel_density_estimate
>>> kde = kernel_density_estimate(returns, n_points=500)
>>> var_95 = kde['x'][np.searchsorted(kde['cdf'], 0.05)]

best_fit_distribution(data, candidates=None)[source]¶

Rank candidate distributions by goodness of fit.

Parameters:

data (Series | ndarray) – Data array or series.
candidates (list[str] | None, default: None) – List of scipy.stats distribution names to test. Defaults to ["norm", "t", "skewnorm", "gennorm", "nct", "johnsonsu"] – a set well-suited for financial returns.

Returns:

distribution, params, ks_statistic, ad_statistic, aic, sorted by AIC (ascending).

Return type:

Example

>>> import numpy as np
>>> data = np.random.default_rng(42).standard_t(df=5, size=1000)
>>> best_fit_distribution(data)

Cointegration¶

Cointegration tests and pairs trading utilities for financial data.

engle_granger(y1, y2, max_lag=None)[source]¶

Engle-Granger two-step cointegration test.

Regresses y1 on y2 via OLS, then tests the residuals for a unit root using the Augmented Dickey-Fuller test.

Parameters:

y1 (Series) – First price series.
y2 (Series) – Second price series.
max_lag (int | None, default: None) – Maximum number of lags for the ADF test. When None, adfuller selects lags automatically via AIC.

Return type:

Returns:

Dictionary with statistic (ADF test statistic), p_value, is_cointegrated (at 5 % significance), hedge_ratio (OLS slope coefficient), and residuals (pd.Series).

johansen(data, det_order=0, k_ar_diff=1)[source]¶

Johansen cointegration test for multiple time series.

Requires the timeseries optional dependency group (provides statsmodels.tsa.vector_ar.vecm).

Parameters:

data (DataFrame) – DataFrame of price series (columns = assets).
det_order (int, default: 0) – Deterministic term order. -1 for no deterministic term, 0 for constant, 1 for linear trend.
k_ar_diff (int, default: 1) – Number of lagged differences in the model.

Return type:

Returns:

Dictionary with trace_stats (array of trace statistics), eigenvalues, critical_values_95 (95 % critical values), and n_cointegrating (number of cointegrating relationships at the 5 % level).

half_life(spread)[source]¶

Estimate the half-life of mean reversion for a spread series.

Fits an OLS regression of the change in spread on the lagged spread level: delta_spread = phi * spread_lag + eps. The half-life is -log(2) / log(1 + phi).

Parameters:: spread (Series) – Spread (or residual) series.
Return type:: float
Returns:: Half-life in the same time units as the spread index. Returns float('inf') when the spread is not mean-reverting.

spread(y1, y2, hedge_ratio=None)[source]¶

Compute the spread between two price series.

When hedge_ratio is None it is estimated via OLS.

Parameters:

y1 (Series) – First price series (the dependent variable).
y2 (Series) – Second price series (the independent variable).
hedge_ratio (float | None, default: None) – Explicit hedge ratio. If None, the ratio is estimated from the data using OLS.

Return type:

Returns:

Spread series (y1 - hedge_ratio * y2).

zscore_signal(spread, window=20)[source]¶

Compute a rolling z-score of the spread for trading signals.

Parameters:

spread (Series) – Spread series.
window (int, default: 20) – Rolling window size for mean and standard deviation.

Return type:

Returns:

Rolling z-score series.

hedge_ratio(y1, y2, method='ols')[source]¶

Estimate the hedge ratio between two price series.

Parameters:

y1 (Series) – First (dependent) price series.
y2 (Series) – Second (independent) price series.
method (str, default: 'ols') – Estimation method — "ols" for ordinary least squares or "tls" for total least squares (orthogonal regression).

Return type:

Returns:

Hedge ratio as a float.

Raises:

ValueError – If method is not recognized.

pairs_backtest_signals(spread, entry_z=2.0, exit_z=0.5)[source]¶

Generate pairs trading signals based on z-score thresholds.

The strategy goes short the spread when z-score > entry_z, goes long when z-score < -entry_z, and exits when the z-score crosses back inside [-exit_z, exit_z].

Parameters:

spread (Series) – Spread series (raw, not z-scored).
entry_z (float, default: 2.0) – Z-score threshold for entry (absolute value).
exit_z (float, default: 0.5) – Z-score threshold for exit (absolute value).

Return type:

Returns:

Signal series with values in {-1, 0, 1}. 1 means long the spread, -1 means short, 0 means flat.

find_cointegrated_pairs(prices_df, significance=0.05)[source]¶

Scan a DataFrame of price series and find all cointegrated pairs.

For each pair of columns, the Engle-Granger cointegration test is applied. Pairs with a p-value below significance are returned.

Parameters:

prices_df (DataFrame) – DataFrame of price series (columns = asset names).
significance (float, default: 0.05) – Significance level for the cointegration test.

Return type:

list[tuple]

Returns:

List of tuples (asset1, asset2, p_value, hedge_ratio) for each cointegrated pair, sorted by p-value ascending.

Statistical Tests¶

Statistical hypothesis tests for financial data.

test_normality(data, method='jarque_bera')[source]¶

Test whether a series is normally distributed.

Parameters:

data (Series) – Data series to test.
method (str, default: 'jarque_bera') – Test method — "jarque_bera" (default), "shapiro", or "dagostino".

Return type:

Returns:

Dictionary with statistic, p_value, and is_normal (at 5% significance level).

Raises:

ValueError – If method is not recognized.

test_stationarity(data, method='adf')[source]¶

Test whether a time series is stationary.

Parameters:

data (Series) – Time series to test.
method (str, default: 'adf') – Test method — "adf" (Augmented Dickey-Fuller, default) or "kpss".

Return type:

Returns:

Dictionary with statistic, p_value, and is_stationary (at 5% significance level).

Raises:

ValueError – If method is not recognized.

test_autocorrelation(data, nlags=10)[source]¶

Ljung-Box test for autocorrelation.

Parameters:

data (Series) – Time series to test.
nlags (int, default: 10) – Number of lags to test.

Return type:

Returns:

Dictionary with statistic (at max lag), p_value, is_autocorrelated (at 5% significance), and the full results DataFrame.

shapiro_wilk(data)[source]¶

Shapiro-Wilk test for normality.

When to use:

When you have fewer than 2000 observations and need a reliable normality assessment (e.g., validating assumptions before parametric VaR, calibrating option pricing models).
As a complement to Jarque-Bera: Shapiro-Wilk catches departures in the center of the distribution that JB (which focuses on moments 3 and 4) may miss.
For validating regression residuals before using t-based confidence intervals.

Mathematical formulation:

\[W = \frac{\left(\sum_{i=1}^n a_i x_{(i)}\right)^2}{\sum_{i=1}^n (x_i - \bar{x})^2}\]

where x_{(i)} are the order statistics and a_i are tabulated constants derived from the expected values and covariance matrix of order statistics from a normal distribution.

How to interpret:

W is in (0, 1]. Values near 1 indicate normality.
Reject normality if p_value < 0.05.
For financial returns, rejection is typical (fat tails), confirming that Gaussian-based risk measures are unreliable.

Parameters:

data (Series | ndarray) – Data series or 1-D array. Sample size should be between 3 and 5000 (scipy limitation).

Returns:

statistic: Shapiro-Wilk W statistic.
p_value: p-value for H0: data is normally distributed.
is_normal: bool, True if p > 0.05.

Return type:

Example

>>> import numpy as np
>>> data = np.random.default_rng(42).normal(0, 1, 200)
>>> result = shapiro_wilk(data)
>>> result["is_normal"]
True

durbin_watson(residuals)[source]¶

Durbin-Watson test for first-order autocorrelation in residuals.

When to use:

After fitting any OLS regression to time-series data (e.g., CAPM beta estimation, factor models). Autocorrelated residuals mean standard errors are biased and t-statistics are unreliable.
As a diagnostic before deciding whether to use Newey-West (HAC) standard errors.
For model validation: significant autocorrelation suggests a missing variable or incorrect functional form.

Mathematical formulation:

\[DW = \frac{\sum_{t=2}^T (e_t - e_{t-1})^2}{\sum_{t=1}^T e_t^2}\]

How to interpret:

DW ≈ 2.0: no autocorrelation.
DW < 2.0: positive autocorrelation (residuals tend to have the same sign as their predecessor).
DW > 2.0: negative autocorrelation.
Rule of thumb: DW < 1.5 or DW > 2.5 indicates significant autocorrelation. For precise inference, compare to the Durbin-Watson tables for dL and dU critical values.

Parameters:

residuals (Series | ndarray) – Regression residuals (1-D array or Series).

Returns:

statistic: Durbin-Watson statistic (range [0, 4]).
interpretation: string describing the result.

Return type:

Example

>>> import numpy as np
>>> residuals = np.random.default_rng(42).normal(0, 1, 100)
>>> result = durbin_watson(residuals)
>>> 1.5 < result["statistic"] < 2.5  # no autocorrelation
True

See also

test_autocorrelation: Ljung-Box test for higher-order autocorrelation.

breusch_pagan(residuals, exog)[source]¶

Breusch-Pagan Lagrange Multiplier test for heteroskedasticity.

Tests whether the variance of regression residuals depends on the values of the independent variables. If heteroskedasticity is present, OLS standard errors are biased and inference is invalid.

When to use:

After OLS regression on financial data where the volatility of returns (and hence residuals) may depend on market conditions, firm size, or other regressors.
To decide between OLS and WLS, or whether to use White/HC robust standard errors.
For validating GARCH model residuals: after fitting a GARCH model, the standardized residuals should be homoskedastic.

Mathematical formulation:

Regress squared residuals e^2 on the original regressors.
The LM statistic is n * R^2 from this auxiliary regression.
Under H0 (homoskedasticity), LM ~ chi^2(k) where k is the number of regressors.

How to interpret:

Low p-value (< 0.05): reject H0, heteroskedasticity is present. Use robust standard errors or WLS.
High p-value: no evidence of heteroskedasticity. OLS inference is valid.

Parameters:

residuals (Series | ndarray) – OLS regression residuals (1-D array or Series).
exog (DataFrame | ndarray) – Design matrix of independent variables used in the original regression (should include constant if one was used).

Returns:

lm_stat: Lagrange Multiplier statistic.
p_value: p-value from chi-squared distribution.
f_stat: F-statistic variant.
f_p_value: p-value from F-distribution.
is_heteroskedastic: bool, True if p_value < 0.05.

Return type:

Example

>>> import numpy as np, statsmodels.api as sm
>>> rng = np.random.default_rng(42)
>>> X = rng.normal(0, 1, (200, 2))
>>> X = sm.add_constant(X)
>>> y = X @ [1, 0.5, -0.3] + rng.normal(0, 1, 200)
>>> from wraquant.stats.regression import ols
>>> resid = ols(y, X, add_constant=False)["residuals"]
>>> result = breusch_pagan(resid, X)
>>> "lm_stat" in result
True

See also

white_test: More general heteroskedasticity test. durbin_watson: Test for autocorrelation instead.

white_test(residuals, exog)[source]¶

White’s test for heteroskedasticity.

When to use:

When you want a comprehensive heteroskedasticity diagnostic that does not assume the variance is a linear function of regressors (which Breusch-Pagan assumes).
When the Breusch-Pagan test fails to reject but you still suspect non-linear heteroskedasticity.
Note: White’s test has lower power than BP when BP’s assumptions are correct, and requires more observations because it estimates more parameters.

Mathematical formulation:

Parameters:

residuals (Series | ndarray) – OLS regression residuals (1-D array or Series).
exog (DataFrame | ndarray) – Design matrix of independent variables (should include constant if used in the original regression).

Returns:

lm_stat: White LM statistic.
p_value: p-value from chi-squared distribution.
f_stat: F-statistic variant.
f_p_value: p-value from F-distribution.
is_heteroskedastic: bool, True if p_value < 0.05.

Return type:

Example

>>> import numpy as np, statsmodels.api as sm
>>> rng = np.random.default_rng(42)
>>> X = rng.normal(0, 1, (200, 2))
>>> X = sm.add_constant(X)
>>> # Heteroskedastic errors: variance depends on X
>>> y = X @ [1, 0.5, -0.3] + rng.normal(0, 1, 200) * (1 + np.abs(X[:, 1]))
>>> from wraquant.stats.regression import ols
>>> resid = ols(y, X, add_constant=False)["residuals"]
>>> result = white_test(resid, X)
>>> "lm_stat" in result
True

See also

breusch_pagan: Simpler but less general heteroskedasticity test.

chow_test(y, X, break_point, add_constant=True)[source]¶

Chow test for structural break at a known break point.

When to use:

To test whether a known event (e.g., a policy announcement, market crash, regulatory change) caused a structural change in the relationship between variables.
As a diagnostic for rolling regression: if the Chow test rejects stability, rolling or regime-switching models are warranted.
To validate that a backtested model’s parameters are stable across in-sample and out-of-sample periods.

Mathematical formulation:

Fit the regression on the full sample, sub-sample 1 (before break), and sub-sample 2 (after break). The F-statistic is:

\[F = \frac{(\text{RSS}_{\text{full}} - \text{RSS}_1 - \text{RSS}_2) / k}{(\text{RSS}_1 + \text{RSS}_2) / (n - 2k)}\]

where k is the number of parameters and n is the total sample size.

How to interpret:

Large F-stat (small p-value < 0.05): reject the null of stable coefficients. A structural break is detected.
Small F-stat: no evidence of a break. The relationship appears stable across the two sub-periods.

Parameters:

y (Series | ndarray) – Dependent variable (1-D array or Series).
X (DataFrame | ndarray) – Independent variables.
break_point (int) – Index (0-based row number) at which to split the sample. Must be at least k + 1 from either end.
add_constant (bool, default: True) – Whether to add an intercept to X.

Returns:

f_stat: Chow F-statistic.
p_value: p-value from the F-distribution.
break_detected: bool, True if p_value < 0.05.

Return type:

Raises:

ValueError – If break_point is too close to the endpoints.

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> X = rng.normal(0, 1, (200, 1))
>>> y = np.concatenate([
...     X[:100] @ [1.0] + rng.normal(0, 0.5, 100),
...     X[100:] @ [3.0] + rng.normal(0, 0.5, 100),
... ])
>>> result = chow_test(y, X, break_point=100)
>>> result["break_detected"]
True

variance_inflation_factor(X)[source]¶

Compute the Variance Inflation Factor (VIF) for each feature.

When to use:

Before running any multiple regression (OLS, factor model, Fama-MacBeth) to check for multicollinearity.
When regression coefficients have unexpected signs or large standard errors despite significant F-statistics.
As a feature selection diagnostic in ML pipelines.

Mathematical formulation:

For each feature X_j, regress it on all other features and compute:

\[\text{VIF}_j = \frac{1}{1 - R_j^2}\]

where R_j^2 is the R-squared from regressing X_j on the remaining features.

How to interpret:

VIF = 1: no collinearity.
VIF < 5: low collinearity, generally acceptable.
5 <= VIF < 10: moderate collinearity, warrants attention.
VIF >= 10: severe collinearity. The coefficient is poorly estimated. Consider removing the feature, combining features, or using regularization (ridge regression).

Parameters:: X (DataFrame) – DataFrame of independent variables (each column is a feature). Do not include an intercept/constant column.
Return type:: Series
Returns:: pd.Series of VIF values indexed by feature name.

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> x1 = rng.normal(0, 1, 200)
>>> x2 = x1 + rng.normal(0, 0.1, 200)  # nearly collinear with x1
>>> x3 = rng.normal(0, 1, 200)
>>> X = pd.DataFrame({"x1": x1, "x2": x2, "x3": x3})
>>> vif = variance_inflation_factor(X)
>>> vif["x1"] > 10  # collinear pair
True

See also

ols: OLS regression where VIF diagnostics are needed.

Factor Analysis¶

Statistical factor analysis for financial data.

Provides PCA-based factor extraction, factor loadings estimation, varimax rotation, factor-mimicking portfolios, and risk decomposition. All implementations use pure numpy/scipy.

pca_factors(returns, n_components=3, method='svd')[source]¶

Extract statistical factors from a returns matrix using PCA.

Parameters:

returns (DataFrame) – DataFrame of asset returns with shape (T, N) where T is the number of observations and N is the number of assets.
n_components (int, default: 3) – Number of principal components to retain.
method (str, default: 'svd') – Decomposition method. "svd" (default) uses numpy.linalg.svd on the demeaned returns; "eig" uses eigendecomposition of the covariance matrix.

Returns:

factors: DataFrame of extracted factors (T, n_components).
loadings: DataFrame of factor loadings (N, n_components).
explained_variance: array of variance explained by each component.
explained_variance_ratio: array of fraction of total variance explained.

Return type:

factor_loadings(returns, factors)[source]¶

Compute factor loadings by regressing each asset’s returns on factors.

For each asset column in returns, an OLS regression is run against factors (with intercept) to obtain betas (loadings), alphas, and R-squared values.

Parameters:

returns (DataFrame | ndarray) – Asset returns matrix (T, N).
factors (DataFrame | ndarray) – Factor returns matrix (T, K).

Returns:

loadings: ndarray (N, K) of factor betas.
alphas: ndarray (N,) of intercepts.
r_squared: ndarray (N,) of regression R-squared values.
residuals: ndarray (T, N) of regression residuals.

Return type:

scree_plot_data(returns)[source]¶

Return eigenvalues and explained variance ratios for a scree plot.

Parameters:

returns (DataFrame | ndarray) – Asset returns matrix (T, N).

Returns:

eigenvalues: array of eigenvalues in descending order.
explained_variance_ratio: array of fraction of total variance explained by each component.
cumulative_variance_ratio: cumulative sum of explained variance ratios.

Return type:

varimax_rotation(loadings, max_iter=500, tol=1e-06)[source]¶

Apply varimax rotation to a factor loadings matrix.

Varimax maximises the sum of variances of squared loadings within each factor, producing a simpler (more interpretable) structure.

Parameters:

loadings (ndarray) – Factor loadings matrix (N, K) where N is the number of variables and K is the number of factors.
max_iter (int, default: 500) – Maximum number of iterations.
tol (float, default: 1e-06) – Convergence tolerance on the rotation criterion change.

Returns:

rotated_loadings: ndarray (N, K) of rotated loadings.
rotation_matrix: ndarray (K, K) orthogonal rotation matrix.
n_iter: number of iterations performed.

Return type:

factor_mimicking_portfolios(returns, characteristics, n_quantiles=5)[source]¶

Build factor-mimicking portfolios via long-short quantile sorts.

Parameters:

returns (DataFrame) – DataFrame of asset returns (T, N).
characteristics (DataFrame) – DataFrame of asset characteristics (T, N) or (N,) for a static sort. If a single row / Series, the same sort is applied to every period. If multiple columns, each column is treated as a separate characteristic, with returns as a single panel.
n_quantiles (int, default: 5) – Number of quantile buckets (default 5).

Return type:

Returns:

DataFrame of factor-mimicking portfolio returns (T, K) where K is the number of characteristics.

risk_factor_decomposition(portfolio_returns, factor_returns)[source]¶

Decompose portfolio risk into factor risk and idiosyncratic risk.

Runs an OLS regression of portfolio_returns on factor_returns and computes the variance attributable to each factor and the residual (idiosyncratic) variance.

Parameters:

portfolio_returns (Series | ndarray) – Portfolio return series (T,).
factor_returns (DataFrame | ndarray) – Factor return matrix (T, K).

Returns:

total_variance: total variance of portfolio returns.
factor_variance: variance explained by the factor model.
idiosyncratic_variance: residual variance.
factor_risk_share: fraction of total variance from factors.
idiosyncratic_risk_share: fraction of total variance from idiosyncratic risk.
betas: regression coefficients (loadings) on each factor.
factor_marginal_contributions: variance contribution of each factor.

Return type:

factor_correlation(factor_returns)[source]¶

Compute the correlation matrix of factors with significance tests.

For each pair of factors, the Pearson correlation and a two-sided p-value (testing H0: rho = 0) are computed.

Parameters:

factor_returns (DataFrame | ndarray) – Factor return matrix (T, K).

Returns:

correlation: ndarray (K, K) correlation matrix.
p_values: ndarray (K, K) of p-values for each pair.

Return type:

common_factors(returns_list, n_components=3)[source]¶

Find common factors shared across multiple asset classes.

Extracts PCA factors from each asset class independently, then performs canonical correlation analysis on the stacked factor scores to identify shared latent factors.

Parameters:

returns_list (list[DataFrame | ndarray]) – List of returns matrices, one per asset class. Each has shape (T, N_i) and all must share the same number of time observations T.
n_components (int, default: 3) – Number of PCA components to extract per asset class before cross-analysis.

Returns:

individual_factors: list of factor arrays, one per asset class.
common_factor_scores: ndarray (T, n_components) of shared factor scores obtained from a second-level PCA on concatenated individual factors.
cross_correlations: correlation matrix between the individual factor sets.

Return type:

fama_french_factors(returns, characteristics, n_quantiles=3)[source]¶

Construct Fama-French style factors from a cross-section of returns.

When to use:

To create custom factors from firm characteristics (e.g., book-to-market, momentum, profitability, investment).
To replicate or extend the Fama-French factor zoo.
To test whether a new characteristic is a priced risk factor (construct the factor, then test its risk premium via fama_macbeth).

Mathematical formulation:

For each period t and each characteristic c:

Sort assets into q quantiles based on the characteristic.
Compute the equal-weighted mean return for the top and bottom quantile.
Factor return = mean(top quantile returns) - mean(bottom quantile returns).

This is a zero-cost, long-short portfolio that isolates the return premium associated with the characteristic.

How to interpret:

A consistently positive factor return means that assets with high values of the characteristic outperform those with low values (and vice versa for negative).
The t-statistic of the mean factor return tests whether the premium is significantly different from zero.
Standard Fama-French uses terciles (3 groups) for the main sort and independent double sorts for intersections (e.g., size x value).

Parameters:

returns (DataFrame) – DataFrame of asset returns (T, N) with a DatetimeIndex and asset names as columns.
characteristics (DataFrame) – DataFrame of asset characteristics (T, N) with the same index and columns as returns, or a DataFrame with (N, K) where N assets are the index and K characteristics are the columns (static sort).
n_quantiles (int, default: 3) – Number of quantile buckets (default 3 = terciles, matching Fama-French convention).

Return type:

Returns:

DataFrame of factor returns (T, K) where K is the number of characteristics.

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> T, N = 100, 30
>>> dates = pd.bdate_range("2020-01-01", periods=T)
>>> assets = [f"s{i}" for i in range(N)]
>>> ret = pd.DataFrame(rng.normal(0, 0.02, (T, N)), index=dates, columns=assets)
>>> bm = pd.DataFrame(rng.uniform(0.5, 3.0, (T, N)), index=dates, columns=assets)
>>> ff = fama_french_factors(ret, bm)
>>> ff.shape[0] == T
True

See also

factor_mimicking_portfolios: General factor-mimicking portfolio: construction.

factor_exposure: Regress returns on constructed factors.

factor_exposure(returns, factor_returns)[source]¶

Regress returns on factor returns to estimate factor exposures (betas).

For each asset (or a single portfolio), runs an OLS regression of returns on factor returns (with intercept) and reports the factor betas, t-statistics, and R-squared.

When to use:

To estimate a portfolio’s or asset’s exposure to known factors (e.g., market, size, value, momentum).
As the first step in factor-based risk decomposition.
To validate that a factor-neutral strategy truly has zero exposure to target factors.

Mathematical formulation:

For each asset i:

\[r_i = \alpha_i + \sum_{k=1}^K \beta_{ik} f_k + \epsilon_i\]

The betas are the OLS coefficients on the factor returns.

How to interpret:

beta > 0: positive exposure to the factor (moves with it).
beta = 0: no exposure.
|t_stat| > 2: the exposure is statistically significant.
R_squared: fraction of return variance explained by the factor model. Higher R-squared means the model is a good fit.

Parameters:

returns (DataFrame | Series) – Asset returns. A Series for a single asset or a DataFrame (T, N) for multiple assets.
factor_returns (DataFrame) – DataFrame of factor returns (T, K).

Returns:

alpha, beta_<factor_name> for each factor, t_<factor_name> for each factor’s t-statistic, and r_squared.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> T = 200
>>> mkt = pd.Series(rng.normal(0, 0.01, T), name="MKT")
>>> smb = pd.Series(rng.normal(0, 0.005, T), name="SMB")
>>> factors = pd.DataFrame({"MKT": mkt, "SMB": smb})
>>> ret = 1.2 * mkt + 0.5 * smb + rng.normal(0, 0.003, T)
>>> result = factor_exposure(pd.Series(ret, name="fund"), factors)
>>> abs(result.loc["fund", "beta_MKT"] - 1.2) < 0.3
True

See also

factor_loadings: Lower-level loadings estimation. factor_risk_decomposition: Risk decomposition from exposures.

factor_risk_decomposition(returns, factor_returns)[source]¶

Decompose total risk into systematic (factor) and idiosyncratic components.

Runs a factor model regression, then separates the total variance of the return series into the portion explained by the factors (systematic risk) and the unexplained residual (idiosyncratic risk).

When to use:

To understand what fraction of a portfolio’s risk comes from common factor exposures vs. stock-specific bets.
For risk budgeting: allocate risk limits to systematic and idiosyncratic components.
To evaluate diversification: a well-diversified portfolio has low idiosyncratic risk relative to total risk.

Mathematical formulation:

\[\text{Var}(r) = \beta' \Sigma_f \beta + \sigma^2_\epsilon\]

where \beta is the vector of factor exposures, \Sigma_f is the factor covariance matrix, and \sigma^2_\epsilon is the idiosyncratic variance.

The R-squared of the regression gives the systematic risk share:

\[R^2 = 1 - \frac{\sigma^2_\epsilon}{\text{Var}(r)}\]

Parameters:

returns (Series) – Return series for a single asset or portfolio.
factor_returns (DataFrame) – DataFrame of factor returns (T, K).

Returns:

systematic_risk: variance explained by factors.
idiosyncratic_risk: residual variance.
total_risk: total return variance.
R_squared: fraction of risk that is systematic.
betas: factor exposure coefficients.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> rng = np.random.default_rng(42)
>>> mkt = rng.normal(0, 0.01, 200)
>>> ret = pd.Series(1.2 * mkt + rng.normal(0, 0.005, 200))
>>> factors = pd.DataFrame({"MKT": mkt})
>>> result = factor_risk_decomposition(ret, factors)
>>> result["R_squared"] > 0.3
True

See also

risk_factor_decomposition: Lower-level decomposition with: marginal contributions.

factor_exposure: Factor exposure estimation.

Factor Models¶

Factor models and attribution for asset pricing.

fama_french_regression(returns, factors_df)[source]¶

Regress asset returns on Fama-French factors.

Parameters:

returns (Series) – Asset return series.
factors_df (DataFrame) – DataFrame of factor returns. Columns are factor names. An optional RF column is the risk-free rate.

Return type:

Returns:

Dictionary with alpha (intercept), betas (dict mapping factor name to coefficient), t_stats (dict mapping name to t-statistic), p_values (dict), and r_squared.

factor_attribution(returns, factor_returns)[source]¶

Decompose returns into factor contributions and specific return.

Runs a regression of returns on factor_returns and attributes the mean return to each factor.

Parameters:

returns (Series) – Asset return series.
factor_returns (DataFrame) – DataFrame of factor return series.

Return type:

Returns:

Dictionary with factor_contributions (dict mapping factor name to its average contribution), specific_return (mean residual return), total_return (mean of returns), and r_squared.

information_coefficient(predictions, returns)[source]¶

Compute the information coefficient (Spearman rank correlation).

The IC measures the predictive power of a signal: the rank correlation between the cross-sectional predictions and subsequent realised returns.

Parameters:

predictions (Series | ndarray) – Predicted values or signals.
returns (Series | ndarray) – Subsequent realised returns.

Return type:

Returns:

Spearman rank correlation coefficient (between -1 and 1).

quantile_analysis(predictions, returns, n_quantiles=5)[source]¶

Analyse returns by prediction quantile.

Sorts observations into quantiles based on predictions and computes summary statistics for each quantile bucket.

Parameters:

predictions (Series) – Predicted values or signals.
returns (Series) – Subsequent realised returns.
n_quantiles (int, default: 5) – Number of quantile buckets (default 5 = quintiles).

Return type:

Returns:

DataFrame indexed by quantile (1 = lowest, n_quantiles = highest) with columns mean_return, std_return, hit_rate (fraction of positive returns), and count.

Dependence¶

Advanced dependence measures for financial data.

Beyond linear correlation, financial risk management requires measures that capture tail dependence, nonlinear relationships, and the full dependence structure between variables. This module provides tools for tail dependence estimation, copula-based dependence modelling, rank correlation matrices, and concordance analysis.

Key concepts:

Tail dependence quantifies the probability that two variables jointly experience extreme values. This is critical for portfolio risk during crises, where correlations spike.
Copulas separate the marginal distributions from the dependence structure, allowing flexible modelling of joint distributions.
Rank correlation is robust to monotonic transformations and outliers, making it suitable for heavy-tailed financial data.
Concordance index measures the agreement between two orderings, useful for model validation and survival analysis.

References

Joe, H. (2014). Dependence Modeling with Copulas.
McNeil, A. J., Frey, R. & Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques and Tools.
Harrell, F. E. et al. (1996). “Multivariable prognostic models.”

tail_dependence_coefficient(x, y, threshold=0.05)[source]¶

Estimate upper and lower tail dependence coefficients from empirical data.

When to use:

To quantify the risk of joint extreme losses in a portfolio.
To assess whether diversification benefits disappear during market crashes (asymmetric tail dependence).
To select the appropriate copula family: Gaussian copulas have zero tail dependence, while Clayton (lower) and Gumbel (upper) copulas can model it.
To compare the tail behaviour of different asset pairs.

Mathematical formulation:

The upper tail dependence coefficient is:

\[\lambda_U = \lim_{q \to 1} P(Y > F_Y^{-1}(q) \mid X > F_X^{-1}(q))\]

The lower tail dependence coefficient is:

\[\lambda_L = \lim_{q \to 0} P(Y \le F_Y^{-1}(q) \mid X \le F_X^{-1}(q))\]

We estimate these empirically using the rank-transformed data (pseudo-observations) and counting joint exceedances.

How to interpret:

lambda = 0: no tail dependence (e.g., Gaussian copula). Diversification holds in the tails.
lambda > 0: positive tail dependence. Extreme events tend to happen together.
lambda_L > lambda_U: lower tail dependence is stronger than upper (common in equity markets — crashes are more contagious than rallies).
Values are bounded in [0, 1].

Parameters:

x (Series | ndarray) – First variable (1-D array or Series).
y (Series | ndarray) – Second variable (1-D array or Series, same length).
threshold (float, default: 0.05) – Quantile threshold for defining “extreme” (default 0.05, meaning the top/bottom 5%).

Returns:

upper_lambda: estimated upper tail dependence coefficient.
lower_lambda: estimated lower tail dependence coefficient.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> # Gaussian data has zero tail dependence
>>> x = rng.normal(0, 1, 5000)
>>> y = 0.7 * x + rng.normal(0, 0.71, 5000)
>>> result = tail_dependence_coefficient(x, y)
>>> 0 <= result["upper_lambda"] <= 1
True
>>> 0 <= result["lower_lambda"] <= 1
True

References

Joe, H. (2014). Dependence Modeling with Copulas, Ch. 2.
McNeil et al. (2015). Quantitative Risk Management, Ch. 7.

See also

copula_selection: Fit copulas that model tail dependence. rank_correlation_matrix: Rank-based dependence matrix.

copula_selection(x, y)[source]¶

Fit multiple copula families and select the best by AIC.

When to use:

To model joint distributions for portfolio risk (e.g., joint simulation of asset returns for VaR/CVaR).
To capture tail dependence or asymmetric dependence that Gaussian models miss.
To select the appropriate copula for bivariate analysis before using it in a larger risk framework.

Copula families fitted:

Gaussian: symmetric, zero tail dependence.
Student-t (approximated): symmetric, positive tail dependence in both tails.
Clayton: lower tail dependence (joint crashes).
Gumbel: upper tail dependence (joint rallies).

Mathematical formulation:

A copula C(u, v) is a joint CDF on [0,1]^2 whose marginals are uniform. By Sklar’s theorem, any joint distribution can be written as:

\[F(x, y) = C(F_X(x), F_Y(y))\]

Each copula family has a parameter theta that controls the strength and shape of dependence. We estimate theta by maximum likelihood on the pseudo-observations.

Parameters:

x (Series | ndarray) – First variable (1-D array or Series).
y (Series | ndarray) – Second variable (1-D array or Series, same length).

Returns:

best_copula: name of the best-fitting copula.
all_fits: DataFrame with columns copula, parameter, log_likelihood, aic, sorted by AIC ascending.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> x = rng.normal(0, 1, 500)
>>> y = 0.7 * x + rng.normal(0, 0.71, 500)
>>> result = copula_selection(x, y)
>>> isinstance(result["all_fits"], pd.DataFrame)
True
>>> result["best_copula"] in ["gaussian", "student_t", "clayton", "gumbel"]
True

References

Joe, H. (2014). Dependence Modeling with Copulas.
Nelsen, R. B. (2006). An Introduction to Copulas.

See also

tail_dependence_coefficient: Empirical tail dependence. rank_correlation_matrix: Rank-based dependence matrix.

rank_correlation_matrix(data, method='spearman')[source]¶

Compute the Spearman rank correlation matrix.

When to use:

When the relationship between variables is monotonic but not necessarily linear.
When data contains outliers that would distort Pearson correlation.
For copula parameter estimation (Spearman’s rho has a direct relationship to many copula parameters).
As a robustness check on Pearson correlation: if the two differ substantially, the relationship may be nonlinear.

Mathematical formulation:

\[\rho_S(X, Y) = \text{Pearson}(\text{rank}(X), \text{rank}(Y))\]

Equivalently:

\[\rho_S = 1 - \frac{6 \sum d_i^2}{n(n^2 - 1)}\]

where d_i is the difference in ranks.

How to interpret:

Same scale as Pearson: [-1, 1].
rho_S > rho_P: the relationship is stronger in the ranks than in the raw values (concave/convex relationship).
rho_S ≈ rho_P: the relationship is approximately linear.

Parameters:

data (DataFrame) – DataFrame with columns as variables and rows as observations.
method (str, default: 'spearman') – Rank correlation method – "spearman" (default) or "kendall".

Return type:

Returns:

Rank correlation matrix as a DataFrame (p x p).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> data = pd.DataFrame(np.random.randn(100, 4), columns=list("ABCD"))
>>> rcm = rank_correlation_matrix(data)
>>> rcm.shape
(4, 4)

See also

correlation_matrix: Pearson correlation matrix. partial_correlation: Partial correlation controlling for others.

concordance_index(predicted, observed)[source]¶

Compute Harrell’s concordance index (C-index).

When to use:

To evaluate the discriminatory power of a predictive model (e.g., a credit scoring model, a default probability model, or an alpha signal).
When you care about ordinal accuracy (ranking) rather than calibration (magnitude).
As an alternative to AUC for continuous outcomes (the C-index generalises AUC to continuous data).

Mathematical formulation:

\[C = \frac{\text{concordant pairs}}{\text{concordant + discordant pairs}}\]

A pair (i, j) is concordant if the predicted and observed orderings agree: (\hat{y}_i > \hat{y}_j \text{ and } y_i > y_j) or (\hat{y}_i < \hat{y}_j \text{ and } y_i < y_j).

How to interpret:

C = 0.5: random (no predictive power).
C = 1.0: perfect concordance (perfect ranking).
C < 0.5: worse than random (model has the sign wrong).
C > 0.7: generally considered acceptable in finance.
C > 0.8: strong discriminatory power.

Parameters:

predicted (Series | ndarray) – Predicted values or scores.
observed (Series | ndarray) – Observed/actual values.

Return type:

Returns:

Concordance index as a float in [0, 1].

Example

>>> import numpy as np
>>> rng = np.random.default_rng(42)
>>> predicted = rng.normal(0, 1, 100)
>>> observed = predicted + rng.normal(0, 0.5, 100)
>>> c = concordance_index(predicted, observed)
>>> c > 0.7  # good concordance
True

References

See also

kendall_tau: Rank correlation (related to C-index).

Robust Statistics¶

Robust statistical methods for financial data.

Standard statistical measures (mean, std, covariance) are sensitive to outliers, fat tails, and data contamination – all common in financial data. This module provides robust alternatives that remain reliable under such conditions.

mad(data, scale='normal')[source]¶

Compute the Median Absolute Deviation (MAD).

MAD is a robust measure of dispersion. Unlike standard deviation, it is not influenced by a few extreme values, making it ideal for financial return distributions with fat tails.

Parameters:

data (Series | ndarray) – Data series or array.
scale (str, default: 'normal') – Scaling factor. Use "normal" (default) so that the result is consistent with standard deviation for normally distributed data. Use None for the raw MAD.

Return type:

Returns:

MAD as a float.

Example

>>> import pandas as pd
>>> returns = pd.Series([0.01, 0.02, -0.01, -0.05, 0.10])
>>> mad(returns)

winsorize(data, lower=0.05, upper=0.05)[source]¶

Cap extreme values at given percentiles (Winsorization).

Winsorization limits extreme values to reduce the influence of outliers without removing observations. This is preferable to trimming when you want to keep the same sample size.

Parameters:

data (Series | ndarray) – Data series or array.
lower (float, default: 0.05) – Fraction to clip on the lower tail (default 5%).
upper (float, default: 0.05) – Fraction to clip on the upper tail (default 5%).

Return type:

Series | ndarray

Returns:

Winsorized data, same type as input.

Example

>>> import pandas as pd
>>> returns = pd.Series([0.01, 0.02, -0.50, 0.03, 0.80])
>>> winsorize(returns, lower=0.1, upper=0.1)

trimmed_mean(data, proportiontocut=0.05)[source]¶

Compute the trimmed mean, excluding extreme observations.

The trimmed mean removes a fraction of the highest and lowest values before computing the average. Use it when the mean is distorted by outliers (e.g., flash-crash returns).

Parameters:

data (Series | ndarray) – Data series or array.
proportiontocut (float, default: 0.05) – Fraction to cut from each tail (default 5%).

Return type:

Returns:

Trimmed mean as a float.

Example

>>> import numpy as np
>>> data = np.array([1, 2, 3, 4, 100])
>>> trimmed_mean(data, proportiontocut=0.2)

trimmed_std(data, proportiontocut=0.05)[source]¶

Compute the standard deviation after trimming extreme values.

Combines trimming with standard deviation computation for a measure of dispersion that is less sensitive to outliers than standard std but retains more information than MAD.

Parameters:

data (Series | ndarray) – Data series or array.
proportiontocut (float, default: 0.05) – Fraction to cut from each tail (default 5%).

Return type:

Returns:

Trimmed standard deviation as a float.

Example

>>> import numpy as np
>>> data = np.array([1, 2, 3, 4, 100])
>>> trimmed_std(data, proportiontocut=0.2)

robust_zscore(data)[source]¶

Compute robust z-scores using median and MAD.

Parameters:: data (Series | ndarray) – Data series or array.
Return type:: Series
Returns:: Robust z-scores as a pd.Series.

Example

>>> import pandas as pd
>>> returns = pd.Series([0.01, 0.02, -0.01, -0.05, 0.50])
>>> robust_zscore(returns)

robust_covariance(data, support_fraction=None)[source]¶

Estimate a robust covariance matrix via Minimum Covariance Determinant.

Parameters:

data (DataFrame) – DataFrame of asset returns (columns = assets).
support_fraction (float | None, default: None) – Fraction of data to use in support (default None lets sklearn choose).

Returns:

covariance: robust covariance matrix (np.ndarray).
location: robust location estimate (np.ndarray).
support_fraction: fraction of data used.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> returns = pd.DataFrame(np.random.randn(100, 3), columns=['A', 'B', 'C'])
>>> robust_covariance(returns)

huber_mean(data, delta=1.5, max_iter=50, tol=1e-08)[source]¶

Compute the Huber M-estimator of location.

Parameters:

data (Series | ndarray) – Data series or array.
delta (float, default: 1.5) – Threshold parameter controlling robustness. Smaller values give more robustness (closer to median). Default 1.5 is a standard choice.
max_iter (int, default: 50) – Maximum number of IRLS iterations.
tol (float, default: 1e-08) – Convergence tolerance for the location estimate.

Return type:

Returns:

Huber location estimate as a float.

Example

>>> import numpy as np
>>> data = np.array([1, 2, 3, 4, 100])
>>> huber_mean(data, delta=1.5)

outlier_detection(data, method='mad', threshold=3.0)[source]¶

Flag outliers using a robust detection method.

Outlier detection is critical in finance for identifying data errors (bad ticks), extreme events, or contaminated observations before computing risk metrics.

Parameters:

data (Series | ndarray) – Data series or array.
method (Literal['mad', 'iqr', 'grubbs'], default: 'mad') –
Detection method: - "mad": Median Absolute Deviation (default). Flag

points whose robust z-score exceeds threshold. Best general-purpose choice for financial data.
- "iqr": Interquartile Range. Flag points outside [Q1 - threshold*IQR, Q3 + threshold*IQR]. Classic box-plot method.
- "grubbs": Grubbs’ test for a single outlier. Tests whether the most extreme value is an outlier assuming approximate normality.
threshold (float, default: 3.0) – Sensitivity parameter (default 3.0). For MAD this is the z-score cutoff; for IQR it is the multiplier.

Returns:

outliers: boolean array (True = outlier).
n_outliers: count of flagged outliers.
method: method used.

Return type: