Risk Management (`wraquant.risk`)¶

Returns:

Annualized Sharpe ratio as a float.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> daily_returns = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> sr = sharpe_ratio(daily_returns, risk_free=0.04)
>>> isinstance(sr, float)
True

Caveats:

Assumes returns are IID and normally distributed; both assumptions are violated in practice.
Penalises upside volatility equally with downside volatility; use sortino_ratio if you only care about downside risk.
Annualisation via sqrt(N) is only exact for IID returns.

See also

sortino_ratio: Uses downside deviation only. information_ratio: Measures alpha relative to a benchmark. wraquant.backtest.tearsheet.comprehensive_tearsheet: Full report including Sharpe. wraquant.stats.descriptive.rolling_sharpe: Time-varying Sharpe ratio.

References

Sharpe (1966), “Mutual Fund Performance”
Bailey & Lopez de Prado (2012), “The Sharpe Ratio Efficient Frontier”

sortino_ratio(returns, risk_free=0.0, periods_per_year=252)[source]¶

Annualized Sortino ratio (downside risk only).

The Sortino ratio replaces total standard deviation with downside deviation – the standard deviation of negative excess returns only. This avoids penalising strategies for upside volatility, making it more appropriate for asymmetric return distributions (which most equity strategies exhibit).

When to use:

Prefer Sortino over Sharpe when returns are skewed or when you only care about downside risk. Particularly useful for options strategies, trend-following, and any strategy with convex payoffs.

Mathematical formulation:

Sortino = (mean(r - r_f) / DD) * sqrt(N)

where DD = sqrt(mean(min(r - r_f, 0)^2)) is the downside deviation (computed as the second lower partial moment).

How to interpret:

Values follow a similar scale to Sharpe, but Sortino is typically higher because the denominator (downside deviation) is smaller than total standard deviation. A Sortino of 2.0 is roughly equivalent to a Sharpe of 1.5 for normally distributed returns.

Parameters:

returns (Series) – Simple return series.
risk_free (float, default: 0.0) – Annual risk-free rate.
periods_per_year (int, default: 252) – Trading periods per year.

Return type:

Returns:

Annualized Sortino ratio. Returns inf when there are no negative excess returns and the mean is positive, and 0.0 when the mean is non-positive.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> daily_returns = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> sr = sortino_ratio(daily_returns, risk_free=0.04)
>>> sr > 0
True

See also

sharpe_ratio: Uses total standard deviation. max_drawdown: Peak-to-trough loss measure.

References

Sortino & van der Meer (1991), “Downside Risk”
Sortino & Satchell (2001), “Managing Downside Risk in Financial Markets”

information_ratio(returns, benchmark)[source]¶

Information ratio (active return / tracking error).

The information ratio measures the average active return (alpha) relative to the variability of that active return (tracking error). It answers: “is the manager consistently adding value, or is the alpha noisy and unreliable?”

When to use:

Use IR when evaluating a strategy relative to a benchmark. Sharpe measures absolute risk-adjusted return; IR measures relative risk-adjusted return. Most relevant for active fund managers benchmarked against an index.

Mathematical formulation:

IR = mean(r_p - r_b) / std(r_p - r_b)

where r_p is the portfolio return and r_b is the benchmark return. This version is not annualized; multiply by sqrt(periods_per_year) to annualize.

How to interpret:

IR < 0: underperforming the benchmark.
0 < IR < 0.3: modest skill; hard to distinguish from luck.
0.3 < IR < 0.5: good active management.
IR > 0.5: exceptional; sustained alpha generation.
IR > 1.0: very rare and likely indicates short sample bias.

Parameters:

returns (Series) – Portfolio return series.
benchmark (Series) – Benchmark return series (same frequency and index).

Return type:

Returns:

Information ratio (not annualized) as a float.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> portfolio = pd.Series(np.random.normal(0.0006, 0.01, 252))
>>> benchmark = pd.Series(np.random.normal(0.0004, 0.009, 252))
>>> ir = information_ratio(portfolio, benchmark)
>>> isinstance(ir, float)
True

See also

sharpe_ratio: Absolute risk-adjusted return. hit_ratio: Fraction of positive-return periods.

References

Grinold & Kahn (2000), “Active Portfolio Management”

max_drawdown(prices)[source]¶

Maximum drawdown from a price series.

Maximum drawdown is the largest peak-to-trough decline in the price series. It measures the worst historical loss an investor would have experienced if they bought at the peak and sold at the trough.

When to use:

Use max drawdown as a “worst case” risk measure. It is more intuitive than VaR for communicating tail risk to non-technical stakeholders. Often used as a hard constraint in portfolio optimisation (e.g., “max drawdown must stay below 15%”).

Mathematical formulation:

MDD = min_t ( (P_t - max_{s<=t} P_s) / max_{s<=t} P_s )

The result is a negative number (or zero if the price series only goes up).

How to interpret:

MDD = -0.10 means the worst peak-to-trough loss was 10%.
MDD = -0.50 means the strategy lost half its value at worst.
For S&P 500 since 1928, the worst MDD was about -0.86 (1929-1932). Post-2000, the worst was about -0.57 (GFC).
A strategy with a Sharpe of 1.0 and max drawdown of -0.30 has a Calmar ratio of about 0.33.

Parameters:: prices (Series) – Price or equity curve series (not returns). Must be positive values.
Return type:: float
Returns:: Maximum drawdown as a negative float (e.g., -0.25 for a 25% drawdown). Returns 0.0 if the series is monotonically increasing.

Example

>>> import pandas as pd
>>> prices = pd.Series([100, 110, 105, 95, 108, 102])
>>> mdd = max_drawdown(prices)
>>> round(mdd, 4)
-0.1364

See also

sharpe_ratio: Risk-adjusted return measure. sortino_ratio: Downside-only risk-adjusted return. wraquant.backtest.tearsheet.comprehensive_tearsheet: Full report with drawdowns. wraquant.stats.descriptive.rolling_drawdown: Time-varying drawdown.

References

Magdon-Ismail & Atiya (2004), “Maximum Drawdown”

hit_ratio(returns)[source]¶

Fraction of positive return periods (win rate).

The hit ratio measures how often the strategy produces a positive return. It is the simplest measure of consistency and is often used alongside payoff ratio (average win / average loss) to characterise a strategy’s profile.

When to use:

Use hit ratio for quick strategy diagnostics. A trend-following strategy typically has a low hit ratio (30-45%) but large average wins. A mean-reversion strategy typically has a high hit ratio (55-70%) but smaller average wins. Neither is inherently better – what matters is the product of hit ratio and payoff ratio.

How to interpret:

0.50 = coin flip; no directional edge.
0.55 = statistically meaningful edge on daily data.
0.60+ = strong edge; verify you are not overfitting.
< 0.40 does not mean a bad strategy if the avg win >> avg loss.

Parameters:: returns (Series) – Simple return series.
Return type:: float
Returns:: Hit ratio as a float between 0 and 1.

Example

>>> import pandas as pd
>>> returns = pd.Series([0.01, -0.005, 0.008, -0.003, 0.012])
>>> hit_ratio(returns)
0.6

See also

sharpe_ratio: Risk-adjusted return measure. information_ratio: Alpha per unit of tracking error.

treynor_ratio(returns, benchmark, risk_free=0.0, periods_per_year=252)[source]¶

Treynor ratio: excess return per unit of systematic (beta) risk.

The Treynor ratio is similar to the Sharpe ratio but uses beta (systematic risk) rather than total standard deviation in the denominator. It is the appropriate performance measure for well-diversified portfolios where specific risk has been diversified away.

When to use:

Use Treynor for comparing portfolios that are part of a larger diversified portfolio (so only systematic risk matters). If the portfolio is the investor’s entire wealth, use Sharpe instead.

Mathematical formulation:

Treynor = (R_p - R_f) / beta_p

where R_p is annualized portfolio return, R_f is the risk-free rate, and beta_p is the portfolio’s beta vs the benchmark.

Parameters:

returns (Series) – Portfolio return series.
benchmark (Series) – Benchmark return series (same frequency and index).
risk_free (float, default: 0.0) – Annual risk-free rate.
periods_per_year (int, default: 252) – Trading periods per year.

Return type:

Returns:

Treynor ratio as a float. Higher is better. Negative indicates the portfolio underperformed the risk-free rate.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> portfolio = 1.2 * market + np.random.normal(0.0001, 0.005, 252)
>>> tr = treynor_ratio(portfolio, market, risk_free=0.04)
>>> isinstance(tr, float)
True

See also

sharpe_ratio: Total risk-adjusted return. jensens_alpha: Excess return above CAPM prediction.

References

Treynor (1965), “How to Rate Management of Investment Funds”

m_squared(returns, benchmark, risk_free=0.0, periods_per_year=252)[source]¶

M-squared (Modigliani-Modigliani) risk-adjusted performance.

M-squared leverages or deleverages the portfolio to match the benchmark’s volatility, then measures the excess return at that risk level. The result is in return units (basis points / percent), making it easier to interpret than the dimensionless Sharpe ratio.

When to use:: Use M-squared when you want to compare two portfolios with different risk levels on a common scale. M-squared answers: “if both portfolios had the same risk as the benchmark, which would earn more?”
Mathematical formulation:: M^2 = SR_p * sigma_b + R_f

Parameters:

returns (Series) – Portfolio return series.
benchmark (Series) – Benchmark return series.
risk_free (float, default: 0.0) – Annual risk-free rate.
periods_per_year (int, default: 252) – Trading periods per year.

Return type:

Returns:

M-squared as an annualized return (float). Positive means the portfolio outperforms the benchmark on a risk-adjusted basis.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> portfolio = pd.Series(np.random.normal(0.0006, 0.01, 252))
>>> benchmark = pd.Series(np.random.normal(0.0004, 0.009, 252))
>>> m2 = m_squared(portfolio, benchmark, risk_free=0.04)
>>> isinstance(m2, float)
True

See also

sharpe_ratio: Dimensionless risk-adjusted return.

References

Modigliani & Modigliani (1997), “Risk-Adjusted Performance”

jensens_alpha(returns, benchmark, risk_free=0.0, periods_per_year=252)[source]¶

Jensen’s alpha: excess return above CAPM-predicted return.

Jensen’s alpha measures the average return of the portfolio in excess of what the Capital Asset Pricing Model (CAPM) predicts given the portfolio’s beta. A positive alpha indicates that the manager generated return beyond what the market risk exposure would explain.

Mathematical formulation:

alpha = R_p - [R_f + beta * (R_m - R_f)]

where R_p is the portfolio return, R_m is the benchmark return, and beta is the portfolio’s beta to the benchmark.

Parameters:

returns (Series) – Portfolio return series.
benchmark (Series) – Benchmark return series.
risk_free (float, default: 0.0) – Annual risk-free rate.
periods_per_year (int, default: 252) – Trading periods per year.

Return type:

Returns:

Annualized Jensen’s alpha as a float.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> stock = 0.0002 + 1.0 * market + np.random.normal(0, 0.005, 252)
>>> alpha = jensens_alpha(stock, market, risk_free=0.04)
>>> isinstance(alpha, float)
True

See also

treynor_ratio: Risk-adjusted return using beta. appraisal_ratio: Alpha per unit of residual risk.

References

Jensen (1968), “The Performance of Mutual Funds in the Period 1945-1964”

appraisal_ratio(returns, benchmark, risk_free=0.0, periods_per_year=252)[source]¶

Appraisal ratio: Jensen’s alpha per unit of residual risk.

The appraisal ratio (also called the Treynor-Black appraisal ratio) measures the manager’s alpha relative to the risk taken to achieve it (residual/idiosyncratic volatility). A high appraisal ratio means the manager generates alpha efficiently.

Mathematical formulation:

AR = alpha / sigma_epsilon

where alpha is Jensen’s alpha and sigma_epsilon is the standard deviation of the regression residuals (annualized).

Parameters:

returns (Series) – Portfolio return series.
benchmark (Series) – Benchmark return series.
risk_free (float, default: 0.0) – Annual risk-free rate.
periods_per_year (int, default: 252) – Trading periods per year.

Return type:

Returns:

Appraisal ratio as a float. Higher is better.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> stock = 0.0003 + 1.0 * market + np.random.normal(0, 0.005, 252)
>>> ar = appraisal_ratio(stock, market, risk_free=0.04)
>>> isinstance(ar, float)
True

See also

jensens_alpha: The numerator of the appraisal ratio. information_ratio: Alpha per unit of tracking error.

References

Treynor & Black (1973), “How to Use Security Analysis to Improve Portfolio Selection”

capture_ratios(returns, benchmark)[source]¶

Up-capture and down-capture ratios.

Capture ratios measure how much of the benchmark’s up and down movements the portfolio captures. An ideal portfolio has high up-capture (>100%) and low down-capture (<100%).

When to use:

Use capture ratios for: - Evaluating defensive vs aggressive positioning: a portfolio

with 90% up-capture and 70% down-capture is defensively positioned and will outperform in bear markets.

Manager selection: compare capture ratios across funds.
Style analysis: growth funds typically have high up-capture and high down-capture; value funds often have lower both.

Mathematical formulation:

Up-capture = (mean(r_p | r_b > 0) / mean(r_b | r_b > 0)) * 100 Down-capture = (mean(r_p | r_b < 0) / mean(r_b | r_b < 0)) * 100

Capture ratio = up-capture / down-capture

Parameters:

returns (Series) – Portfolio return series.
benchmark (Series) – Benchmark return series (same frequency and index).

Returns:

up_capture (float) – Percentage of benchmark’s up movements captured (>100 = amplified).
down_capture (float) – Percentage of benchmark’s down movements captured (<100 = dampened losses).
capture_ratio (float) – up_capture / down_capture. Values > 1 indicate the portfolio adds value through asymmetric participation.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> benchmark = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> portfolio = 0.8 * benchmark + np.random.normal(0.0001, 0.005, 252)
>>> caps = capture_ratios(portfolio, benchmark)
>>> caps["up_capture"] > 0
True

See also

treynor_ratio: Systematic risk-adjusted return.

value_at_risk(returns, confidence=0.95, method='historical')[source]¶

Estimate Value-at-Risk (VaR).

VaR answers the question: “With X% confidence, what is the maximum loss I should expect over one period?” More precisely, VaR is the (1 - confidence) quantile of the return distribution, flipped to a positive loss number.

When to use:

Use VaR for regulatory reporting, margin calculations, and setting position limits. Choose historical VaR when you have enough data (>500 observations) and want to capture fat tails without distributional assumptions. Choose parametric VaR when data is scarce or when you need analytical sensitivities (e.g., delta-normal VaR for a derivatives book).

Mathematical formulation:

Historical: VaR_alpha = -quantile(returns, 1 - alpha) Parametric: VaR_alpha = -(mu + sigma * Phi^{-1}(1 - alpha))

where alpha is the confidence level, mu and sigma are the sample mean and standard deviation, and Phi^{-1} is the standard normal inverse CDF.

How to interpret:

A 95% daily VaR of 0.02 means: “on 95% of days, the portfolio loses less than 2%. On the remaining 5% of days, the loss exceeds 2%.” VaR says nothing about how much worse the loss can be beyond the threshold – that is what CVaR captures.

Parameters:

returns (Series) – Simple return series (e.g., daily percentage changes).
confidence (float, default: 0.95) – Confidence level (e.g., 0.95 for 95%, 0.99 for 99%). Basel III uses 0.99; internal risk management often uses 0.95.
method (str, default: 'historical') –
Estimation method: - "historical" – empirical quantile (non-parametric,

default). No distributional assumption; captures fat tails.
- "parametric" – Gaussian assumption. Smooth but underestimates tail risk for leptokurtic returns.

Return type:

Returns:

VaR as a positive float representing the loss threshold. For example, 0.025 means a 2.5% loss at the given confidence level.

Raises:

ValueError – If method is not recognized.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0, 0.01, 1000))
>>> var_95 = value_at_risk(returns, confidence=0.95)
>>> var_95 > 0
True

Caveats:

VaR is not sub-additive: the VaR of a portfolio can exceed the sum of individual VaRs. Use conditional_var for a coherent measure.
Historical VaR is sensitive to the sample window; recent crises dominate short windows.
Parametric VaR severely underestimates tail risk for fat- tailed distributions (equities, credit).

See also

conditional_var: Expected loss beyond the VaR threshold (CVaR). garch_var: GARCH-based time-varying VaR using conditional volatility. wraquant.vol.models.garch_fit: Fit GARCH model for conditional vol. wraquant.risk.stress.stress_test_returns: Scenario-based analysis.

References

Jorion (2006), “Value at Risk: The New Benchmark”
Basel Committee on Banking Supervision (2019), “Minimum capital requirements for market risk”

conditional_var(returns, confidence=0.95, method='historical')[source]¶

Estimate Conditional VaR (Expected Shortfall / CVaR).

CVaR answers: “given that the loss exceeds VaR, what is the expected loss?” It captures the severity of tail losses, not just their threshold. Unlike VaR, CVaR is a coherent risk measure (Artzner et al. 1999) – it satisfies sub-additivity, meaning the CVaR of a portfolio is at most the sum of individual CVaRs.

When to use:

CVaR is preferred over VaR for: - Portfolio optimisation (mean-CVaR optimisation is convex). - Regulatory capital under Basel IV / FRTB. - Any situation where you care about tail severity, not just

tail frequency.

Use historical CVaR with long samples (>1000 obs) and parametric CVaR when you need smooth gradients or have short data.

Mathematical formulation:

Historical: CVaR_alpha = -mean(returns | returns <= VaR_quantile) Parametric: CVaR_alpha = -(mu - sigma * phi(z_alpha) / (1 - alpha))

where z_alpha = Phi^{-1}(1 - alpha), phi is the standard normal PDF, and Phi is the CDF.

How to interpret:

A 95% daily CVaR of 0.035 means: “on the worst 5% of days, the average loss is 3.5%.” CVaR is always >= VaR at the same confidence level. For normal distributions, 95% CVaR is about 1.25x the 95% VaR. For fat-tailed distributions, the ratio is much larger – this ratio itself is a useful diagnostic of tail heaviness.

Parameters:

returns (Series) – Simple return series.
confidence (float, default: 0.95) – Confidence level (e.g., 0.95 for 95%).
method (str, default: 'historical') –
Estimation method: - "historical" – mean of returns in the tail

(default). Non-parametric; captures fat tails.
- "parametric" – Gaussian formula. Smooth but underestimates tail risk for heavy-tailed distributions.

Return type:

Returns:

CVaR as a positive float representing the expected tail loss. For example, 0.035 means an expected loss of 3.5% in the tail.

Raises:

ValueError – If method is not recognized.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0, 0.01, 1000))
>>> cvar = conditional_var(returns, confidence=0.95)
>>> var = value_at_risk(returns, confidence=0.95)
>>> cvar >= var  # CVaR is always >= VaR
True

See also

value_at_risk: The VaR threshold itself. wraquant.risk.monte_carlo.importance_sampling_var: Variance-

reduced tail estimation.

References

Artzner et al. (1999), “Coherent Measures of Risk”
Rockafellar & Uryasev (2000), “Optimization of Conditional Value-at-Risk”

garch_var(returns, alpha=0.05, vol_model='GARCH', p=1, q=1, dist='normal', horizon=1)[source]¶

Value at Risk using GARCH conditional volatility.

Combines GARCH volatility forecasting with parametric VaR, producing time-varying risk estimates that adapt to current market conditions. Superior to static VaR in volatile markets.

The conditional VaR at time t is:: VaR_t = -mu + sigma_t * z_alpha

where sigma_t is the GARCH-forecasted volatility and z_alpha is the quantile of the fitted error distribution.

Parameters:

returns (Series | ndarray) – Return series.
alpha (float, default: 0.05) – Significance level (0.05 = 95% VaR).
vol_model (str, default: 'GARCH') – GARCH variant (“GARCH”, “EGARCH”, “GJR”).
p (int, default: 1) – GARCH lag order.
q (int, default: 1) – ARCH lag order.
dist (str, default: 'normal') – Error distribution (“normal”, “t”, “skewt”).
horizon (int, default: 1) – Forecast horizon in periods.

Returns:

var (pd.Series) – Time-varying VaR series.
cvar (pd.Series) – Time-varying CVaR/ES series.
conditional_vol (pd.Series) – GARCH conditional volatility.
breaches (pd.Series) – Boolean where actual loss exceeded VaR.
breach_rate (float) – Fraction of breaches (should be ~alpha).
garch_params (dict) – Fitted GARCH parameters.

Return type:

Example

>>> from wraquant.risk.var import garch_var
>>> result = garch_var(returns, alpha=0.05, vol_model="GJR", dist="t")
>>> print(f"Breach rate: {result['breach_rate']:.3f} (target: 0.050)")

greeks_var(portfolio_greeks, spot, vol, rf=0.0, dt=0.003968253968253968, spot_shock=0.01, vol_shock=0.01, n_scenarios=10000, confidence=0.95, seed=None)[source]¶

VaR approximation using portfolio Greeks (delta-gamma-vega).

Approximates portfolio P&L using a second-order Taylor expansion with Greeks from the price/ module, then estimates VaR from the resulting P&L distribution. This bridges price/ (Greeks computation) and risk/ (VaR estimation).

The P&L approximation is:

PnL approx delta * dS + 0.5 * gamma * dS^2 + vega * d_sigma + theta * dt

where dS and d_sigma are simulated from normal distributions with standard deviations spot_shock * spot and vol_shock respectively.

When to use:: Use delta-gamma-vega VaR for options portfolios where the P&L is nonlinear in the underlying. Standard (delta-only) VaR underestimates risk for portfolios with significant gamma or vega exposure. Full revaluation VaR is more accurate but much slower; this method is a fast approximation.

Parameters:

portfolio_greeks (dict[str, float]) – Dictionary with portfolio-level Greeks. Required keys: 'delta', 'gamma'. Optional keys: 'vega', 'theta'.
spot (float) – Current spot price of the underlying.
vol (float) – Current implied volatility (annualised, e.g. 0.20 for 20%).
rf (float, default: 0.0) – Risk-free rate (annualised). Used for drift in spot dynamics.
dt (float, default: 0.003968253968253968) – Time step as a fraction of a year (default 1/252 for one trading day).
spot_shock (float, default: 0.01) – Standard deviation of the spot return used for simulation (default 0.01 = 1%). Typically set to vol * sqrt(dt) for a realistic one-day shock.
vol_shock (float, default: 0.01) – Standard deviation of the volatility change (default 0.01 = 1 vol point).
n_scenarios (int, default: 10000) – Number of Monte Carlo scenarios (default 10,000).
confidence (float, default: 0.95) – VaR confidence level (default 0.95).
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

'var' (float) – Estimated VaR (positive = loss).
'cvar' (float) – Estimated CVaR / Expected Shortfall.
'mean_pnl' (float) – Mean P&L across scenarios.
'std_pnl' (float) – Standard deviation of P&L.
'delta_component' (float) – VaR contribution from delta.
'gamma_component' (float) – VaR contribution from gamma.
'vega_component' (float) – VaR contribution from vega.
'theta_component' (float) – Deterministic theta P&L.

Return type:

Example

>>> greeks = {'delta': 100, 'gamma': -50, 'vega': 200, 'theta': -10}
>>> result = greeks_var(greeks, spot=100, vol=0.20, seed=42)
>>> result['var'] > 0
True
>>> result['cvar'] >= result['var']
True

Notes

For a single option, compute Greeks with wraquant.price.greeks and pass them here. For a portfolio, sum the Greeks across positions first.

See also

value_at_risk: Standard return-based VaR. wraquant.price.greeks: Compute option Greeks.

portfolio_volatility(weights, cov_matrix)[source]¶

Compute portfolio volatility from weights and covariance matrix.

Parameters:

weights (ndarray) – Asset weight vector (n,).
cov_matrix (ndarray) – Covariance matrix (n, n).

Return type:

Returns:

Portfolio volatility (standard deviation).

risk_contribution(weights, cov_matrix)[source]¶

Compute each asset’s risk contribution to portfolio volatility.

Parameters:

weights (ndarray) – Asset weight vector (n,).
cov_matrix (ndarray) – Covariance matrix (n, n).

Return type:

Returns:

Array of fractional risk contributions that sum to 1.

diversification_ratio(weights, cov_matrix)[source]¶

Compute the diversification ratio.

The diversification ratio is the ratio of the weighted average of individual volatilities to the portfolio volatility. A value of 1 means no diversification benefit.

Parameters:

weights (ndarray) – Asset weight vector (n,).
cov_matrix (ndarray) – Covariance matrix (n, n).

Return type:

Returns:

Diversification ratio (>= 1).

component_var(weights, returns, alpha=0.05)[source]¶

Component Value-at-Risk: per-asset contribution to portfolio VaR.

Decomposes portfolio VaR into additive per-asset contributions using the Euler (marginal) decomposition. The sum of component VaRs equals the portfolio VaR. This tells you where the tail risk is concentrated.

When to use:

Use component VaR for: - Identifying which assets dominate portfolio tail risk. - Setting per-asset risk limits. - Reporting risk contributions to portfolio managers and risk

committees.

Mathematical formulation:

Component VaR_i = w_i * (partial VaR / partial w_i)

Under the delta-normal approximation: CVaR_i = w_i * (Sigma @ w)_i / sigma_p * VaR_p

Parameters:

weights (ndarray) – Portfolio weight vector (n_assets,).
returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
alpha (float, default: 0.05) – Significance level (0.05 = 95% VaR).

Return type:

Returns:

pd.Series of per-asset VaR contributions, indexed by asset names. Sum equals the portfolio VaR.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "A": np.random.normal(0.0005, 0.01, 252),
...     "B": np.random.normal(0.0003, 0.015, 252),
... })
>>> weights = np.array([0.6, 0.4])
>>> cvar = component_var(weights, returns, alpha=0.05)
>>> cvar.sum() > 0  # total VaR is positive
True

See also

marginal_var: Sensitivity of VaR to weight changes. incremental_var: VaR change from adding/removing an asset.

marginal_var(weights, cov, alpha=0.05)[source]¶

Marginal VaR: sensitivity of portfolio VaR to weight changes.

Marginal VaR measures how much portfolio VaR changes for a small (infinitesimal) change in the weight of each asset. It is the gradient of portfolio VaR with respect to weights.

When to use:

Use marginal VaR for: - Position sizing: assets with high marginal VaR should have

smaller positions.

Optimisation: marginal VaR should be equal across assets at the optimal portfolio (risk parity condition).
Hedging: the hedge ratio is proportional to the marginal VaR.

Mathematical formulation:

Marginal VaR_i = dVaR/dw_i = z_alpha * (Sigma @ w)_i / sigma_p

Parameters:

weights (ndarray) – Portfolio weight vector (n_assets,).
cov (ndarray) – Covariance matrix (n_assets x n_assets).
alpha (float, default: 0.05) – Significance level (0.05 = 95% VaR).

Return type:

Returns:

np.ndarray of marginal VaR values per asset.

Example

>>> import numpy as np
>>> cov = np.array([[0.0004, 0.0001], [0.0001, 0.0009]])
>>> weights = np.array([0.6, 0.4])
>>> mvar = marginal_var(weights, cov, alpha=0.05)
>>> len(mvar) == 2
True

See also

component_var: Additive VaR decomposition (weight * marginal VaR).

incremental_var(weights, returns, alpha=0.05)[source]¶

Incremental VaR: change in portfolio VaR from adding each asset.

For each asset, computes the difference between the portfolio VaR with and without that asset (reallocating its weight proportionally to remaining assets). This measures the discrete impact of each position on tail risk.

When to use:: Use incremental VaR when deciding whether to add or remove a position. Unlike marginal VaR (which is an infinitesimal measure), incremental VaR captures the full nonlinear impact including diversification effects.

Parameters:

weights (ndarray) – Portfolio weight vector (n_assets,).
returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
alpha (float, default: 0.05) – Significance level (0.05 = 95% VaR).

Return type:

Returns:

np.ndarray of incremental VaR values per asset. Positive means adding the asset increases portfolio VaR (adds risk); negative means it reduces VaR (diversification benefit).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "A": np.random.normal(0.001, 0.01, 252),
...     "B": np.random.normal(0.0005, 0.008, 252),
... })
>>> weights = np.array([0.6, 0.4])
>>> ivar = incremental_var(weights, returns, alpha=0.05)
>>> len(ivar) == 2
True

See also

component_var: Euler-based additive decomposition. marginal_var: Infinitesimal sensitivity.

risk_budgeting(cov, target_risk=None)[source]¶

Find portfolio weights that achieve target risk contributions.

Risk budgeting finds the weights such that each asset’s risk contribution (Euler decomposition) matches a target budget. With equal targets (default), this is the risk parity portfolio.

When to use:

Use risk budgeting for: - Risk parity portfolio construction (equal risk contribution). - Custom risk allocation (e.g., 60% risk from equities, 40%

from bonds, regardless of capital allocation).

Avoiding concentration: risk-budgeted portfolios avoid overweighting high-volatility assets.

Mathematical formulation:

Find w such that: w_i * (Sigma @ w)_i / sigma_p = b_i * sigma_p

where b_i is the target risk budget (sum to 1).

Parameters:

cov (ndarray) – Covariance matrix (n x n).
target_risk (ndarray | None, default: None) – Target risk contribution vector (sums to 1). If None, uses equal risk contributions (1/n for each asset).

Returns:

weights (np.ndarray) – Optimal portfolio weights.
risk_contributions (np.ndarray) – Achieved risk contributions (should match target).
portfolio_vol (float) – Portfolio volatility.
converged (bool) – Whether the optimiser converged.

Return type:

Example

>>> import numpy as np
>>> cov = np.array([[0.04, 0.006], [0.006, 0.01]])
>>> result = risk_budgeting(cov)
>>> np.allclose(result["risk_contributions"], 0.5, atol=0.05)
True

See also

wraquant.risk.portfolio.risk_contribution: Compute risk: contributions for given weights.

wraquant.opt.portfolio: Full portfolio optimisation suite.

References

Maillard, Roncalli & Teiletche (2010), “The Properties of Equally Weighted Risk Contribution Portfolios”

concentration_ratio(weights, cov)[source]¶

Herfindahl concentration ratio of risk contributions.

Measures how concentrated portfolio risk is across assets using the Herfindahl-Hirschman Index (HHI) of risk contributions. An equally risk-contributed portfolio has HHI = 1/n (minimum concentration).

When to use:

Use concentration ratio to: - Detect hidden risk concentrations even when capital weights

look diversified. A portfolio with equal weights can still have concentrated risk if one asset is much more volatile.

Monitor risk concentration over time.
Compare portfolios: lower concentration ratio = more diversified risk.

Mathematical formulation:

CR = sum(rc_i^2) where rc_i is asset i’s fractional risk contribution (sum to 1.0).

CR = 1/n for equal risk contribution; CR = 1.0 for single-asset.

Parameters:

weights (ndarray) – Portfolio weight vector (n_assets,).
cov (ndarray) – Covariance matrix (n_assets x n_assets).

Return type:

Returns:

Herfindahl concentration ratio between 1/n and 1.0.

Example

>>> import numpy as np
>>> cov = np.array([[0.04, 0.0], [0.0, 0.04]])
>>> weights = np.array([0.5, 0.5])
>>> cr = concentration_ratio(weights, cov)
>>> abs(cr - 0.5) < 0.01  # equal vol + equal weight -> equal risk
True

See also

diversification_ratio: Alternative diversification metric.

tracking_error(returns, benchmark)[source]¶

Active risk metrics relative to a benchmark.

Tracking error (TE) is the standard deviation of the active return (portfolio return minus benchmark return). It measures how much the portfolio’s performance deviates from the benchmark.

When to use:: Use tracking error for: - Index tracking: target TE < 50bp for passive strategies. - Active management: typical TE of 2-8% for active equity funds. - Risk budgeting: allocate TE budget across portfolio managers.

Parameters:

returns (Series) – Portfolio return series.
benchmark (Series) – Benchmark return series (same frequency and index).

Returns:

tracking_error (float) – Annualized tracking error.
information_ratio (float) – Annualized active return / tracking error.
active_return (float) – Annualized mean active return.
max_active_drawdown (float) – Worst cumulative active return drawdown.
active_return_std (float) – Daily active return standard deviation (non-annualized).

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> portfolio = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> benchmark = pd.Series(np.random.normal(0.0004, 0.009, 252))
>>> result = tracking_error(portfolio, benchmark)
>>> result["tracking_error"] > 0
True

See also

active_share: Weight-based difference from benchmark. wraquant.risk.metrics.information_ratio: Simpler IR calculation.

active_share(weights, benchmark_weights)[source]¶

Active share: weight-based deviation from benchmark.

Active share measures how different a portfolio’s holdings are from its benchmark. It is computed as half the sum of absolute weight differences.

When to use:

Use active share to classify portfolio management style: - Active share < 20%: closet indexer (charging active fees for

passive exposure).

20-60%: moderate active.
60-80%: genuinely active.
> 80%: concentrated active or different investment universe.

Mathematical formulation:

Active Share = (1/2) * sum_i |w_i - w_bench_i|

Parameters:

weights (ndarray) – Portfolio weight vector.
benchmark_weights (ndarray) – Benchmark weight vector (same length).

Return type:

Returns:

Active share as a float between 0 and 1.

Example

>>> import numpy as np
>>> portfolio = np.array([0.4, 0.3, 0.2, 0.1])
>>> benchmark = np.array([0.25, 0.25, 0.25, 0.25])
>>> as_ = active_share(portfolio, benchmark)
>>> 0 <= as_ <= 1
True

See also

tracking_error: Return-based deviation from benchmark.

References

Cremers & Petajisto (2009), “How Active Is Your Fund Manager?”

rolling_beta(returns, benchmark, window=60)[source]¶

Rolling OLS beta of asset returns against a benchmark.

Computes the ordinary least squares regression slope of returns on benchmark over a rolling window. This is the standard approach for tracking how an asset’s market sensitivity evolves over time.

When to use:

Use rolling beta to: - Monitor regime changes in market exposure (beta rising during

sell-offs indicates contagion).

Calibrate dynamic hedging ratios (e.g., beta-hedge a long position with index futures).
Detect structural breaks in a strategy’s factor exposure.

Mathematical formulation:

beta_t = Cov(r, b; t-w:t) / Var(b; t-w:t)

where r is the asset return, b is the benchmark return, and w is the rolling window size.

Parameters:

returns (Series) – Asset return series (e.g., daily simple returns).
benchmark (Series) – Benchmark return series (same frequency and aligned index). Typically a broad market index (S&P 500, MSCI World).
window (int, default: 60) – Rolling window size in periods. 60 trading days (~3 months) is standard for equity beta. Use 120-252 for more stable estimates; use 20-40 for faster-reacting estimates.

Return type:

Returns:

pd.Series of rolling beta values, indexed to match returns. The first window - 1 values are NaN (insufficient data).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> stock = 1.2 * market + np.random.normal(0, 0.005, 252)
>>> beta = rolling_beta(stock, market, window=60)
>>> abs(beta.iloc[-1] - 1.2) < 0.3
True

See also

ewma_beta: Exponentially weighted alternative (no fixed window). conditional_beta: Separate up/down market betas.

References

Ang & Chen (2007), “Asymmetric Correlations of Equity Portfolios”

blume_adjusted_beta(raw_beta)[source]¶

Blume-adjusted beta (mean-reversion adjustment).

Blume (1971) documented that betas regress toward 1.0 over time. The adjustment applies the empirical relationship:

adjusted_beta = 0.33 + 0.67 * raw_beta

This is the adjustment used by Bloomberg and most commercial risk systems when reporting “adjusted beta.”

When to use:: Use Blume adjustment for forward-looking beta estimates (e.g., cost of equity in CAPM). The raw historical beta is a biased predictor of future beta; the Blume adjustment reduces this bias.

Parameters:: raw_beta (float) – Historical OLS beta estimate.
Return type:: float
Returns:: Adjusted beta as a float, shrunk toward 1.0.

Example

>>> blume_adjusted_beta(1.5)
1.335
>>> blume_adjusted_beta(0.5)
0.665
>>> blume_adjusted_beta(1.0)
1.0

See also

vasicek_adjusted_beta: Bayesian shrinkage with uncertainty. rolling_beta: Source of the raw beta input.

References

Blume (1971), “On the Assessment of Risk”, Journal of Finance
Blume (1975), “Betas and Their Regression Tendencies”

vasicek_adjusted_beta(raw_beta, cross_sectional_mean=1.0, raw_se=0.2, prior_se=0.3)[source]¶

Vasicek Bayesian shrinkage beta adjustment.

Combines the sample beta with a prior (typically the cross-sectional mean beta of 1.0) using a precision-weighted average. Assets with imprecise beta estimates (high standard error) are shrunk more toward the prior.

When to use:

Use Vasicek adjustment when you have an estimate of beta’s standard error (e.g., from OLS regression). It is more principled than Blume’s fixed-weight adjustment because the shrinkage intensity adapts to estimation uncertainty.

Mathematical formulation:

adjusted_beta = (prior_se^2 / (prior_se^2 + raw_se^2)) * raw_beta

(raw_se^2 / (prior_se^2 + raw_se^2)) * cross_sectional_mean

Parameters:

raw_beta (float) – Historical OLS beta estimate.
cross_sectional_mean (float, default: 1.0) – Prior mean beta (cross-sectional average). Typically 1.0 for market beta.
raw_se (float, default: 0.2) – Standard error of the raw beta estimate from OLS regression. Higher values cause more shrinkage.
prior_se (float, default: 0.3) – Standard deviation of the cross-sectional beta distribution. Represents uncertainty in the prior.

Return type:

Returns:

Vasicek-adjusted beta as a float.

Example

>>> vasicek_adjusted_beta(1.5, cross_sectional_mean=1.0, raw_se=0.2, prior_se=0.3)
1.3461538461538463
>>> # High SE -> more shrinkage toward prior
>>> vasicek_adjusted_beta(1.5, raw_se=0.5, prior_se=0.3)
1.1323529411764706

See also

blume_adjusted_beta: Simpler fixed-weight adjustment.

References

Vasicek (1973), “A Note on Using Cross-Sectional Information in Bayesian Estimation of Security Betas”

dimson_beta(returns, benchmark, lags=1)[source]¶

Dimson beta for illiquid or thinly traded assets.

Standard OLS beta underestimates the true beta of assets that trade infrequently, because non-synchronous trading introduces measurement error. The Dimson (1979) correction runs a multiple regression of asset returns on contemporaneous and lagged benchmark returns, then sums all coefficients to recover the “true” beta.

When to use:

Use Dimson beta for: - Small-cap and micro-cap stocks with thin trading. - Private equity or real estate benchmarked against a public index. - Emerging market assets with liquidity constraints. A significant difference between total_beta and the contemporaneous beta suggests non-synchronous trading effects.

Mathematical formulation:

r_t = alpha + beta_0 * b_t + beta_1 * b_{t-1} + … + beta_k * b_{t-k} + eps

Dimson beta = sum(beta_0, beta_1, …, beta_k)

Parameters:

returns (Series) – Asset return series.
benchmark (Series) – Benchmark return series (same frequency and index).
lags (int, default: 1) – Number of lagged benchmark terms to include. 1 is standard for daily data; use 2-3 for very illiquid assets.

Returns:

total_beta (float) – Sum of all lag coefficients (the Dimson-adjusted beta).
lag_betas (list[float]) – Individual coefficients for each lag (index 0 = contemporaneous).
alpha (float) – Regression intercept.
r_squared (float) – R-squared of the multiple regression.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> benchmark = pd.Series(np.random.normal(0.0005, 0.01, 300))
>>> # Illiquid asset: reacts with a lag
>>> returns = 0.5 * benchmark + 0.4 * benchmark.shift(1).fillna(0) + \\
...     np.random.normal(0, 0.005, 300)
>>> result = dimson_beta(returns, benchmark, lags=1)
>>> result["total_beta"] > result["lag_betas"][0]
True

See also

rolling_beta: Standard OLS beta (assumes synchronous trading). ewma_beta: Exponentially weighted beta.

References

Dimson (1979), “Risk Measurement When Shares are Subject to Infrequent Trading”, Journal of Financial Economics

conditional_beta(returns, benchmark)[source]¶

Conditional (asymmetric) beta: separate up-market and down-market betas.

Standard beta assumes symmetric sensitivity to the benchmark. In practice, many assets have higher beta in down markets than up markets (the “leverage effect” and flight-to-quality dynamics). Conditional beta splits the regression into up-market days (benchmark > 0) and down-market days (benchmark <= 0).

When to use:

Use conditional beta to: - Assess downside protection: an asset with low downside beta

and high upside beta is a desirable portfolio component.

Detect asymmetric risk exposure: if downside_beta >> upside_beta, the asset amplifies losses more than gains.
Evaluate hedge fund or options-like payoff profiles.

Mathematical formulation:

Up-market: r_t = alpha_up + beta_up * b_t + eps, for b_t > 0 Down-market: r_t = alpha_down + beta_down * b_t + eps, for b_t <= 0

Parameters:

returns (Series) – Asset return series.
benchmark (Series) – Benchmark return series (same frequency and index).

Returns:

upside_beta (float) – Beta in up-market periods.
downside_beta (float) – Beta in down-market periods.
beta_asymmetry (float) – downside_beta - upside_beta. Positive means the asset is more sensitive to down moves.
n_up (int) – Number of up-market observations.
n_down (int) – Number of down-market observations.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = pd.Series(np.random.normal(0.0005, 0.01, 500))
>>> stock = 1.0 * market + np.random.normal(0, 0.005, 500)
>>> result = conditional_beta(stock, market)
>>> isinstance(result["upside_beta"], float)
True

See also

rolling_beta: Time-varying beta (not conditional on direction). ewma_beta: Exponentially weighted beta.

References

Pettengill, Sundaram & Mathur (1995), “The Conditional Relation Between Beta and Returns”
Ang & Chen (2002), “Asymmetric Correlations of Equity Portfolios”

ewma_beta(returns, benchmark, halflife=60)[source]¶

Exponentially weighted moving average (EWMA) beta.

Uses exponentially weighted covariance and variance to compute a time-varying beta that adapts to recent market conditions faster than a fixed rolling window. More recent observations receive exponentially higher weight.

When to use:

Use EWMA beta when you need a smooth, responsive beta estimate that adapts quickly to regime changes. Compared to rolling beta: - EWMA has no “cliff effect” (old observations do not drop out

abruptly).

EWMA adapts faster to structural breaks (smaller halflife).
EWMA is smoother (no window-edge artifacts).

Mathematical formulation:

beta_t = EWCov(r, b; lambda) / EWVar(b; lambda)

where lambda = 1 - exp(-ln(2) / halflife) is the decay factor.

Parameters:

returns (Series) – Asset return series.
benchmark (Series) – Benchmark return series (same frequency and index).
halflife (int, default: 60) – Decay halflife in periods. 60 days is standard. Shorter halflife (20-30) reacts faster but is noisier. Longer halflife (90-120) is smoother but lags.

Return type:

Returns:

pd.Series of EWMA beta values.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> stock = 1.3 * market + np.random.normal(0, 0.005, 252)
>>> beta = ewma_beta(stock, market, halflife=60)
>>> abs(beta.iloc[-1] - 1.3) < 0.4
True

See also

rolling_beta: Fixed-window alternative. conditional_beta: Direction-dependent beta.

References

RiskMetrics Technical Document (1996), J.P. Morgan

factor_risk_model(returns, factors)[source]¶

Regress asset returns on factors and decompose total risk.

Fits a multivariate OLS regression of returns on the provided factor returns, then decomposes total variance into the portion explained by factors (systematic risk) and the residual (specific/idiosyncratic risk).

When to use:

Use this function when you have a set of candidate factors (market, value, momentum, macro variables) and want to understand how much of the return variation they explain. The factor_risk / specific_risk split guides hedging decisions: hedge systematic risk with factor instruments; accept specific risk if you believe in the alpha.

Mathematical formulation:

r_t = alpha + B * f_t + eps_t

Total variance = B’ * Sigma_f * B + sigma_eps^2 Factor risk share = B’ * Sigma_f * B / Total variance Specific risk share = sigma_eps^2 / Total variance

Parameters:

returns (Series | DataFrame) – Asset return series (pd.Series for one asset, pd.DataFrame for multiple assets – uses first column).
factors (DataFrame) – DataFrame of factor returns with columns as factor names and a compatible index.

Returns:

betas (dict[str, float]) – Factor loadings (regression coefficients). Positive beta = positive exposure.
alpha (float) – Regression intercept (excess return not explained by factors).
factor_risk (float) – Fraction of total variance explained by factors (0 to 1).
specific_risk (float) – Fraction of total variance from idiosyncratic sources (1 - factor_risk).
r_squared (float) – R-squared of the regression.
residual_vol (float) – Annualized volatility of residuals.
contributions (dict[str, float]) – Each factor’s individual contribution to systematic variance.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> mkt = np.random.normal(0.0005, 0.01, 252)
>>> smb = np.random.normal(0, 0.005, 252)
>>> stock = 1.1 * mkt + 0.3 * smb + np.random.normal(0, 0.005, 252)
>>> result = factor_risk_model(
...     pd.Series(stock),
...     pd.DataFrame({"MKT": mkt, "SMB": smb}),
... )
>>> result["factor_risk"] > 0.5
True

See also

statistical_factor_model: PCA-based (no prior on factors). fama_french_regression: Specialised Fama-French interface.

References

Menchero (2011), “The Barra Risk Model Handbook”

statistical_factor_model(returns, n_factors=3)[source]¶

PCA-based statistical factor model with risk decomposition.

Extracts latent factors from the cross-section of asset returns using Principal Component Analysis (PCA). The first principal component typically captures market-wide movements; subsequent components capture sector, style, and other systematic effects.

When to use:: Use statistical factor models when you do not have a prior on which factors drive returns. PCA discovers the dominant sources of covariation. Useful for: - Constructing factor-mimicking portfolios. - Dimensionality reduction before portfolio optimisation. - Identifying hidden risk concentrations.

Parameters:

returns (DataFrame) – DataFrame of asset returns (columns = assets, rows = dates). Should have at least n_factors + 1 columns.
n_factors (int, default: 3) – Number of principal components to extract. 3-5 is typical for equity portfolios.

Returns:

factors (pd.DataFrame) – Extracted factor return series (columns: PC1, PC2, …).
loadings (np.ndarray) – Factor loadings matrix (n_assets x n_factors).
explained_variance (np.ndarray) – Variance explained by each factor.
explained_variance_ratio (np.ndarray) – Fraction of total variance explained by each factor.
cumulative_variance_ratio (np.ndarray) – Cumulative fraction of variance explained.
factor_risk (float) – Total fraction of variance explained by all extracted factors.
specific_risk (float) – Fraction of variance not explained.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = np.random.normal(0, 0.01, 252)
>>> returns = pd.DataFrame({
...     f"asset_{i}": market * (0.5 + i * 0.2) + np.random.normal(0, 0.005, 252)
...     for i in range(5)
... })
>>> result = statistical_factor_model(returns, n_factors=2)
>>> result["factor_risk"] > 0.3
True

See also

factor_risk_model: When you know which factors to use.

References

Connor & Korajczyk (1986), “Performance Measurement with the Arbitrage Pricing Theory”

fama_french_regression(returns, factors_df)[source]¶

Fama-French factor regression with full diagnostics.

Regresses asset returns on named Fama-French factors (e.g., Mkt-RF, SMB, HML, RMW, CMA, Mom). Reports alpha, betas, t-statistics, and R-squared. The alpha represents the return not explained by factor exposures – a positive, statistically significant alpha indicates genuine skill.

When to use:: Use for performance attribution and alpha measurement. The classic 3-factor model (Mkt, SMB, HML) is the minimum; the 5-factor model adds RMW (profitability) and CMA (investment). Add Mom (momentum) for the 6-factor model.

Parameters:

returns (Series) – Asset or portfolio return series (excess of risk-free rate if the factors are excess returns).
factors_df (DataFrame) – DataFrame of factor returns. Column names should be descriptive (e.g., “Mkt-RF”, “SMB”, “HML”).

Returns:

alpha (float) – Jensen’s alpha (intercept).
betas (dict[str, float]) – Factor loadings.
t_stats (dict[str, float]) – t-statistics for each coefficient (including alpha under key “alpha”).
p_values (dict[str, float]) – p-values for each coefficient.
r_squared (float) – R-squared.
adj_r_squared (float) – Adjusted R-squared.
residual_vol (float) – Annualized residual volatility.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> mkt = np.random.normal(0.0005, 0.01, 252)
>>> smb = np.random.normal(0, 0.005, 252)
>>> hml = np.random.normal(0, 0.005, 252)
>>> stock = 0.0001 + 1.1 * mkt + 0.3 * smb - 0.2 * hml + \\
...     np.random.normal(0, 0.003, 252)
>>> factors = pd.DataFrame({"Mkt-RF": mkt, "SMB": smb, "HML": hml})
>>> result = fama_french_regression(pd.Series(stock), factors)
>>> abs(result["betas"]["Mkt-RF"] - 1.1) < 0.2
True

See also

factor_risk_model: General factor regression with risk decomposition.

References

Fama & French (1993), “Common Risk Factors in the Returns on Stocks and Bonds”
Fama & French (2015), “A Five-Factor Asset Pricing Model”

factor_contribution(weights, factor_betas, factor_cov)[source]¶

Decompose portfolio factor risk into per-factor contributions.

Given portfolio weights, a matrix of factor loadings, and the factor covariance matrix, computes how much each factor contributes to total portfolio factor risk (variance).

When to use:

Use after estimating a factor model to understand which factors dominate portfolio risk. This guides factor hedging decisions: if 80% of portfolio risk comes from the market factor, you can hedge with index futures to dramatically reduce risk.

Mathematical formulation:

Portfolio factor variance = w’ * B * Sigma_f * B’ * w

Factor i contribution = w’ * B_i * (Sigma_f * B’ * w)_i / total_var

Parameters:

weights (ndarray) – Portfolio weight vector (n_assets,).
factor_betas (ndarray) – Factor loading matrix (n_assets x n_factors). Each row is an asset’s factor exposures.
factor_cov (ndarray) – Factor covariance matrix (n_factors x n_factors).

Returns:

total_factor_var (float) – Total portfolio factor variance.
total_factor_vol (float) – Square root of factor variance.
factor_contributions (np.ndarray) – Each factor’s variance contribution (sums to total_factor_var).
factor_pct_contributions (np.ndarray) – Percentage contributions (sum to 1.0).

Return type:

Example

>>> import numpy as np
>>> weights = np.array([0.3, 0.3, 0.4])
>>> betas = np.array([[1.0, 0.5], [1.2, -0.3], [0.8, 0.1]])
>>> factor_cov = np.array([[0.0004, 0.00005], [0.00005, 0.0001]])
>>> result = factor_contribution(weights, betas, factor_cov)
>>> result["total_factor_var"] > 0
True

See also

factor_risk_model: Estimate factor betas from return data. statistical_factor_model: Extract latent factors via PCA.

cornish_fisher_var(returns, alpha=0.05)[source]¶

Cornish-Fisher expansion VaR (skewness and kurtosis adjusted).

The Cornish-Fisher expansion modifies the standard normal quantile to account for skewness (S) and excess kurtosis (K) of the return distribution. This produces a more accurate VaR than parametric (Gaussian) VaR for non-normal distributions.

When to use:

Use Cornish-Fisher VaR when: - Returns are detectably non-normal (skewness != 0 or kurtosis != 3). - You want a quick analytical adjustment without fitting a full

distribution (e.g., Student-t or EVT).

The sample is too short for reliable historical VaR but long enough to estimate skewness/kurtosis (>100 observations).

Mathematical formulation:

z_cf = z + (z^2 - 1) * S/6 + (z^3 - 3z) * K/24 - (2z^3 - 5z) * S^2/36

CF-VaR = -(mu + sigma * z_cf)

where z = Phi^{-1}(alpha), S = skewness, K = excess kurtosis.

How to interpret:

Compare cf_var to normal_var. If cf_var > normal_var, the distribution has fatter left tails than normal (typical for equities). The adjustment_factor (cf_var / normal_var) tells you how much the normal VaR underestimates tail risk.

Parameters:

returns (Series) – Simple return series.
alpha (float, default: 0.05) – Significance level (0.05 = 95% VaR).

Returns:

cf_var (float) – Cornish-Fisher adjusted VaR (positive number = loss).
normal_var (float) – Standard parametric (Gaussian) VaR.
z_cf (float) – Adjusted quantile.
z_normal (float) – Standard normal quantile.
skewness (float) – Sample skewness.
excess_kurtosis (float) – Sample excess kurtosis.
adjustment_factor (float) – cf_var / normal_var.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0, 0.01, 1000))
>>> result = cornish_fisher_var(returns, alpha=0.05)
>>> result["cf_var"] > 0
True

See also

wraquant.risk.var.value_at_risk: Historical and parametric VaR. tail_ratio_analysis: Tail shape diagnostics.

References

Cornish & Fisher (1937), “Moments and Cumulants in the Specification of Distributions”
Maillard (2012), “A User’s Guide to the Cornish Fisher Expansion”

expected_shortfall_decomposition(weights, returns, alpha=0.05)[source]¶

Decompose Expected Shortfall (CVaR) into per-asset contributions.

Each asset’s contribution to portfolio ES is computed as its average return on the days when the portfolio return is in the worst alpha tail. These contributions are additive (they sum to total portfolio ES).

When to use:

Use ES decomposition for: - Identifying which assets drive tail losses. - Setting per-asset ES limits. - Comparing tail-risk concentration to normal-market risk

concentration.

Mathematical formulation:

ES_i = w_i * E[r_i | r_p <= VaR_alpha(r_p)]

where r_p = w’ @ r is the portfolio return.

Parameters:

weights (ndarray) – Portfolio weight vector (n_assets,).
returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
alpha (float, default: 0.05) – Significance level (0.05 = worst 5% of days).

Return type:

Returns:

pd.Series of per-asset ES contributions. Sum equals portfolio ES (as a positive number).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "A": np.random.normal(0.0005, 0.01, 500),
...     "B": np.random.normal(0.0003, 0.015, 500),
... })
>>> weights = np.array([0.6, 0.4])
>>> es = expected_shortfall_decomposition(weights, returns, alpha=0.05)
>>> es.sum() > 0  # positive = loss
True

See also

component_var: Euler decomposition of VaR. cornish_fisher_var: Skewness-adjusted VaR.

conditional_drawdown_at_risk(returns, alpha=0.05)[source]¶

Conditional Drawdown at Risk (CDaR).

CDaR is the average of the worst alpha fraction of drawdowns in the return series. It is analogous to CVaR (Expected Shortfall) but operates on drawdowns rather than returns. CDaR is a coherent risk measure and is used in drawdown-constrained portfolio optimisation.

When to use:: Use CDaR when drawdown is a primary risk constraint (e.g., hedge funds with max drawdown mandates). CDaR penalises sustained drawdowns, not just point-in-time losses. A portfolio optimised to minimise CDaR will have better drawdown recovery properties than one optimised for VaR.

Parameters:

returns (Series) – Simple return series.
alpha (float, default: 0.05) – Fraction of worst drawdowns to average (0.05 = worst 5%).

Return type:

Returns:

CDaR as a positive float (e.g., 0.15 = average worst-5% drawdown is 15%).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0.0005, 0.01, 500))
>>> cdar = conditional_drawdown_at_risk(returns, alpha=0.05)
>>> cdar >= 0
True

See also

drawdown_at_risk: Quantile-based drawdown measure (DaR). wraquant.risk.metrics.max_drawdown: Single worst drawdown.

References

Chekhlov, Uryasev & Zabarankin (2005), “Drawdown Measure in Portfolio Optimization”

tail_ratio_analysis(returns)[source]¶

Tail ratio analysis with interpretation.

The tail ratio is the ratio of the right tail (gains) to the absolute value of the left tail (losses) at a given percentile. A ratio > 1 means the distribution has fatter right tails (gains are larger than losses at the extremes). A ratio < 1 means fatter left tails (losses are larger than gains).

When to use:

Use tail ratio analysis to: - Assess payoff asymmetry: trend-following should have tail ratio > 1

(large gains, small frequent losses).

Detect negative skew: mean-reversion and short vol strategies typically have tail ratio < 1.
Compare strategies beyond Sharpe ratio.

Parameters:

returns (Series) – Simple return series.

Returns:

tail_ratio (float) – 95th percentile / abs(5th percentile).
right_tail (float) – 95th percentile return.
left_tail (float) – 5th percentile return.
tail_ratio_99 (float) – 99th/1st percentile ratio.
skewness (float) – Sample skewness.
excess_kurtosis (float) – Sample excess kurtosis.
interpretation (str) – Human-readable assessment.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0, 0.01, 1000))
>>> result = tail_ratio_analysis(returns)
>>> result["tail_ratio"] > 0
True

See also

cornish_fisher_var: Skewness-adjusted VaR.

drawdown_at_risk(returns, alpha=0.05)[source]¶

Drawdown at Risk (DaR): worst alpha-quantile drawdown.

DaR is to drawdowns what VaR is to returns. It is the alpha-percentile of the drawdown distribution – i.e., the drawdown that is exceeded only alpha% of the time.

When to use:: Use DaR when setting drawdown limits: - “With 95% confidence, the drawdown will not exceed DaR.” - Useful for fund prospectuses and investor communications. - More intuitive than VaR for many stakeholders because

drawdowns are easier to understand than daily P&L.

Parameters:

returns (Series) – Simple return series.
alpha (float, default: 0.05) – Significance level (0.05 = 95th percentile drawdown).

Return type:

Returns:

DaR as a positive float (e.g., 0.12 = 12% drawdown).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0.0005, 0.01, 500))
>>> dar = drawdown_at_risk(returns, alpha=0.05)
>>> dar >= 0
True

See also

conditional_drawdown_at_risk: Average of worst drawdowns (CDaR). wraquant.risk.metrics.max_drawdown: Single worst drawdown.

crisis_drawdowns(returns, top_n=5)[source]¶

Identify the top N drawdowns with full lifecycle metrics.

Scans the return series for the largest peak-to-trough drawdowns and reports start date, trough date, recovery date, duration, and magnitude for each.

When to use:

Use crisis drawdowns for: - Investor reporting: show the worst historical losses and

recovery times.

Strategy evaluation: compare drawdown profiles across strategies.
Risk limit calibration: set max drawdown limits based on historical experience.

Parameters:

returns (Series) – Simple return series with a DatetimeIndex (or integer index).
top_n (int, default: 5) – Number of largest drawdowns to return.

Returns:

start – Date the drawdown began (peak).
trough – Date of maximum drawdown.
end – Date the drawdown recovered (or last date if still in drawdown).
drawdown – Magnitude of drawdown (negative number).
days_to_trough – Trading days from start to trough.
days_to_recovery – Trading days from trough to recovery (NaN if not recovered).
total_days – Total drawdown duration.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> idx = pd.bdate_range("2020-01-01", periods=500)
>>> returns = pd.Series(np.random.normal(0.0003, 0.01, 500), index=idx)
>>> dd = crisis_drawdowns(returns, top_n=3)
>>> len(dd) <= 3
True

See also

drawdown_attribution: Which assets caused the drawdowns. wraquant.risk.metrics.max_drawdown: Single worst drawdown.

event_impact(returns, event_dates, window=10)[source]¶

Measure portfolio returns around specific events.

For each event date, extracts the returns in a window before and after the event and computes cumulative return, max drawdown, and volatility within each window.

When to use:

Use event impact analysis for: - Post-mortem: “how did the portfolio react to the Fed rate hike?” - Event studies: systematic analysis of recurring events

(earnings, FOMC, NFP).

Scenario planning: calibrate stress scenarios based on actual event impacts.

Parameters:

returns (Series) – Return series with a DatetimeIndex.
event_dates (list[str]) – List of event date strings (ISO format, e.g., “2020-03-16”). Dates not in the index are matched to the nearest available date.
window (int, default: 10) – Number of trading days before and after the event to analyse.

Returns:

pre_cumulative (float) – Cumulative return in the window before the event.
post_cumulative (float) – Cumulative return in the window after the event.
event_day_return (float) – Return on the event day itself.
pre_vol (float) – Volatility in the pre-event window.
post_vol (float) – Volatility in the post-event window.
total_impact (float) – Cumulative return over the full window (pre + event + post).

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> idx = pd.bdate_range("2020-01-01", periods=252)
>>> returns = pd.Series(np.random.normal(0.0005, 0.01, 252), index=idx)
>>> result = event_impact(returns, ["2020-03-16", "2020-06-15"], window=5)
>>> len(result) >= 1
True

See also

wraquant.risk.stress.historical_stress_test: Replay known crises. crisis_drawdowns: Top drawdown periods.

contagion_analysis(returns_df, crisis_dates)[source]¶

Compare normal vs. crisis-period correlations to detect contagion.

Contagion occurs when correlations increase during stress periods beyond what would be expected from higher volatility alone. This function computes the correlation matrix in normal and crisis periods and tests for statistically significant increases.

When to use:

Use contagion analysis for: - Evaluating diversification reliability: do correlations spike

when you need diversification most?

Stress testing: adjust portfolio correlations based on empirically observed crisis behaviour.
Regime-aware portfolio construction: allocate less to assets that become highly correlated during crises.

Parameters:

returns_df (DataFrame) – Multi-asset return DataFrame with DatetimeIndex.
crisis_dates (tuple[str, str]) – Tuple of (start_date, end_date) strings defining the crisis period.

Returns:

normal_corr (pd.DataFrame) – Correlation matrix during non-crisis period.
crisis_corr (pd.DataFrame) – Correlation matrix during the crisis period.
corr_change (pd.DataFrame) – Change in correlation (crisis - normal).
avg_normal_corr (float) – Average off-diagonal correlation in normal period.
avg_crisis_corr (float) – Average off-diagonal correlation in crisis period.
contagion_detected (bool) – True if average crisis correlation significantly exceeds normal.
n_normal (int) – Number of normal-period observations.
n_crisis (int) – Number of crisis-period observations.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> idx = pd.bdate_range("2019-01-01", periods=500)
>>> returns = pd.DataFrame({
...     "A": np.random.normal(0.0005, 0.01, 500),
...     "B": np.random.normal(0.0003, 0.012, 500),
... }, index=idx)
>>> result = contagion_analysis(returns, ("2020-02-01", "2020-06-01"))
>>> "contagion_detected" in result
True

See also

wraquant.risk.stress.joint_stress_test: Apply correlation shocks.

References

Forbes & Rigobon (2002), “No Contagion, Only Interdependence: Measuring Stock Market Comovements”

drawdown_attribution(returns_df, weights)[source]¶

Attribute portfolio drawdowns to individual asset contributions.

For each point in the drawdown, decomposes the portfolio’s loss from peak into per-asset contributions. This shows which assets are responsible for the drawdown at each point in time.

When to use:

Use drawdown attribution for: - Post-mortem analysis: “which position caused the 2020 drawdown?” - Risk monitoring: track per-asset drawdown contributions in

real time.

Portfolio construction: identify assets that consistently contribute to drawdowns and consider hedging or removing them.

Parameters:

returns_df (DataFrame) – Multi-asset return DataFrame (columns = assets).
weights (ndarray) – Portfolio weight vector aligned with columns.

Returns:

portfolio_dd – Total portfolio drawdown at each point.
One column per asset showing that asset’s contribution to the drawdown.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> idx = pd.bdate_range("2020-01-01", periods=252)
>>> returns = pd.DataFrame({
...     "A": np.random.normal(0.0005, 0.01, 252),
...     "B": np.random.normal(0.0003, 0.015, 252),
... }, index=idx)
>>> weights = np.array([0.6, 0.4])
>>> attr = drawdown_attribution(returns, weights)
>>> "portfolio_dd" in attr.columns
True

See also

crisis_drawdowns: Identify top drawdown periods. wraquant.risk.stress.marginal_stress_contribution: Stress-based

attribution.

monte_carlo_var(returns, weights, n_sims=10000, confidence=0.95)[source]¶

Estimate portfolio VaR via Monte Carlo simulation.

Draws from a multivariate normal distribution fitted to historical returns.

Parameters:

returns (DataFrame) – DataFrame of asset returns (columns = assets).
weights (ndarray) – Portfolio weight vector.
n_sims (int, default: 10000) – Number of simulation draws.
confidence (float, default: 0.95) – Confidence level.

Return type:

Returns:

VaR as a positive float.

stress_test(returns, weights, shocks)[source]¶

Compute the portfolio loss under a stress scenario.

Each entry in shocks maps an asset name to a shocked return. Assets not in shocks use their historical mean.

Parameters:

returns (DataFrame) – DataFrame of asset returns (columns = assets).
weights (ndarray) – Portfolio weight vector.
shocks (dict[str, float]) – Mapping of asset name to shocked return value.

Return type:

Returns:

Portfolio return under the stress scenario (typically negative).

stress_test_returns(returns, scenarios)[source]¶

Apply user-defined additive shock scenarios to a return series.

Each scenario name maps to an additive shift applied uniformly to all returns. The function computes the stressed mean, stressed VaR (5th percentile), and stressed CVaR for every scenario. This is the simplest stress test and is useful for quick what-if analysis.

When to use:: Use when you want to evaluate the impact of a uniform adverse shift in returns. For example, “what if every daily return is 10bp worse due to a funding cost shock?” For more realistic scenario analysis, use historical_stress_test (replays real crises) or joint_stress_test (simultaneous multi-factor shocks).
How to interpret:: Compare stressed_var_95 and stressed_cvar_95 across scenarios to identify which shock level pushes your portfolio into unacceptable loss territory. If a moderate shock (-5%) already produces a severe stressed CVaR, the portfolio has insufficient risk budget.

Parameters:

returns (Series | DataFrame) – Historical return series (Series) or multi-asset returns (DataFrame). For DataFrames, the cross-asset mean is used.
scenarios (dict[str, float]) – Mapping of scenario name to additive return shock (e.g. {"crash": -0.10, "boom": 0.05}). A shock of -0.10 subtracts 10% from every observation.

Returns:

"scenario_results" – dict mapping scenario name to a dict with "stressed_mean", "stressed_var_95" (5th percentile of stressed returns), "stressed_cvar_95" (mean of returns below the 5th percentile).
"base_mean" – mean of the original (unstressed) returns.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> result = stress_test_returns(returns, {"mild": -0.005, "severe": -0.02})
>>> result["scenario_results"]["severe"]["stressed_mean"] < result["base_mean"]
True

See also

historical_stress_test: Replay known crisis periods. vol_stress_test: Scale volatility by multipliers. joint_stress_test: Simultaneous multi-factor shocks.

historical_stress_test(returns, crisis_periods=None)[source]¶

Test portfolio returns against known historical crisis periods.

Replays the portfolio through actual historical crises and reports cumulative return, max drawdown, and mean daily return for each period. This is the most intuitive form of stress testing because the scenarios are real events that stakeholders can relate to.

When to use:

Use historical stress testing for: - Board and regulator presentations (“how would we have

performed in the GFC?”).

Identifying whether the portfolio’s risk profile has improved or deteriorated relative to past crises.
Calibrating position limits against known worst cases.

How to interpret:

Compare cumulative_return and max_drawdown across crises. A portfolio that survived the GFC with only -15% cumulative return has very different tail risk from one that lost -45%. The periods_found list tells you which crises overlap with your data – crises not found are silently skipped.

Built-in crisis periods (used when crisis_periods is None):

GFC 2008: 2008-09-01 to 2009-03-31
COVID 2020: 2020-02-19 to 2020-03-23
Dot-Com 2000: 2000-03-10 to 2002-10-09
Euro Debt 2011: 2011-07-01 to 2011-11-30
Taper Tantrum 2013: 2013-05-22 to 2013-09-05
Volmageddon 2018: 2018-02-02 to 2018-02-08
Flash Crash 2010: 2010-05-06

Parameters:

returns (Series | DataFrame) – Return series with a DatetimeIndex.
crisis_periods (dict[str, tuple[str, str]] | None, default: None) – Mapping of crisis name to (start, end) date strings. Periods not covered by the data are skipped.

Returns:

"crisis_results" – dict mapping crisis name to a dict with "cumulative_return" (compounded return over the crisis), "max_drawdown" (worst peak-to-trough within the crisis), "mean_daily_return", "n_days".
"periods_found" – list of crisis names that overlap with the data.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> idx = pd.bdate_range("2008-01-01", "2009-12-31")
>>> returns = pd.Series(np.random.normal(-0.001, 0.02, len(idx)), index=idx)
>>> result = historical_stress_test(returns)
>>> "gfc_2008" in result["periods_found"]
True

See also

stress_test_returns: User-defined additive shocks. reverse_stress_test: Find scenarios that produce a target loss.

vol_stress_test(returns, vol_shocks=None)[source]¶

Stress test by scaling return volatility with multipliers.

Demeaned returns are scaled by each multiplier, then the mean is re-added. This preserves the expected return while increasing (or decreasing) dispersion. The technique is useful for asking “what happens to VaR and CVaR if volatility doubles?”

When to use:: Use volatility stress tests to: - Assess margin adequacy under elevated vol regimes. - Calibrate dynamic position sizing rules. - Compare the portfolio’s sensitivity to vol scaling

(a diversified portfolio should be less sensitive than a concentrated one).
How to interpret:: The stressed_vol should scale linearly with the multiplier (by construction). The key outputs are stressed_var_95 and stressed_cvar_95: if doubling vol (multiplier 2.0) causes CVaR to more than double, the portfolio has convex (nonlinear) tail exposure.

Parameters:

returns (Series | DataFrame) – Historical return series.
vol_shocks (list[float] | None, default: None) – List of volatility multipliers (e.g. [1.5, 2.0, 3.0]). Defaults to [1.5, 2.0, 2.5, 3.0].

Returns:

"vol_results" – dict mapping multiplier (as string) to "stressed_vol", "stressed_var_95", "stressed_cvar_95", "stressed_mean".
"base_vol" – volatility of the original returns.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0.0005, 0.01, 500))
>>> result = vol_stress_test(returns, vol_shocks=[1.5, 2.0])
>>> result["vol_results"]["2.0"]["stressed_vol"] > result["base_vol"]
True

See also

stress_test_returns: Additive shock scenarios. joint_stress_test: Combined vol, spot, and correlation shocks.

spot_stress_test(prices, spot_shocks=None)[source]¶

Shift spot (price) levels by specified percentage amounts.

Each shock is applied as a multiplicative factor to the final price (e.g., -0.10 means a 10% drop from the last price). This is useful for mark-to-market stress testing of current positions.

When to use:: Use spot stress tests for: - Options portfolio Greeks analysis (delta P&L under spot move). - Margin calculation under adverse spot scenarios. - Reporting to counterparties (“what is my exposure if the

underlying drops 20%?”).
How to interpret:: The shocked_price shows the resulting price level after the shock. For a DataFrame (multi-asset), the same percentage shock is applied to each asset’s last price. price_change is the absolute dollar change.

Parameters:

prices (Series | DataFrame) – Price series or DataFrame of asset prices.
spot_shocks (list[float] | None, default: None) – List of percentage shocks (e.g. [-0.30, -0.20, -0.10, 0.10, 0.20]). Defaults to [-0.30, -0.20, -0.10, -0.05, 0.05, 0.10].

Returns:

"spot_results" – dict mapping shock (as string) to "shocked_price", "price_change", "pct_change".
"base_price" – the last observed price.

Return type:

Example

>>> import pandas as pd
>>> prices = pd.Series([100.0, 102.0, 101.0, 103.0])
>>> result = spot_stress_test(prices, spot_shocks=[-0.10, 0.10])
>>> result["spot_results"]["-0.1"]["shocked_price"]
92.7

See also

vol_stress_test: Scale volatility by multipliers. sensitivity_ladder: P&L sensitivity to a single factor.

sensitivity_ladder(portfolio_returns, factor_returns, shock_range=None)[source]¶

Compute portfolio P&L across a range of factor shocks.

Fits a linear regression of portfolio returns on a single factor, then uses the estimated beta to project the portfolio P&L impact at each shock level. The result is a “ladder” – a table of factor values and corresponding portfolio returns.

When to use:

Use sensitivity ladders to: - Understand how exposed the portfolio is to a single risk

factor (e.g., S&P 500, 10Y yield, oil price).

Construct hedging ratios (the beta tells you how much factor exposure to neutralise).
Present risk to traders and PMs in an intuitive format.

Mathematical formulation:

Step 1: Fit r_p = alpha + beta * r_f + epsilon via OLS. Step 2: For each shock s, estimate P&L = alpha + beta * s.

How to interpret:

The ladder maps each factor shock to the estimated portfolio return. A high beta means the portfolio is very sensitive to the factor. r_squared tells you how much of the portfolio’s variance is explained by this factor; if R^2 is low (<0.3), the ladder is unreliable because other factors dominate.

Parameters:

portfolio_returns (Series) – Portfolio return series.
factor_returns (Series) – Single-factor return series (same index).
shock_range (ndarray | list[float] | None, default: None) – Array of factor shock values (e.g. np.linspace(-0.10, 0.10, 21)). Defaults to np.linspace(-0.10, 0.10, 21).

Returns:

"ladder" – dict mapping shock level (float) to estimated portfolio P&L.
"beta" – regression beta (sensitivity).
"alpha" – regression intercept (return when factor = 0).
"r_squared" – R-squared of the regression.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> factor = pd.Series(np.random.normal(0, 0.01, 252))
>>> portfolio = 0.8 * factor + np.random.normal(0, 0.005, 252)
>>> result = sensitivity_ladder(portfolio, factor)
>>> abs(result["beta"] - 0.8) < 0.2  # beta close to true value
True

See also

spot_stress_test: Direct price-level shocks. joint_stress_test: Multi-factor simultaneous shocks.

reverse_stress_test(returns, target_loss, n_sims=10000, seed=None)[source]¶

Find scenarios that produce at least the specified target loss.

Reverse stress testing inverts the usual question: instead of “what is the loss under scenario X?”, it asks “what scenarios produce a loss of at least Y?” This is a regulatory requirement under ICAAP/SREP and is valuable for identifying the portfolio’s breaking point.

When to use:

Use reverse stress tests when you need to: - Identify the conditions under which the portfolio breaches

a risk limit (e.g., -20% annual loss).

Satisfy regulatory requirements for reverse stress testing.
Understand how “unlikely” a catastrophic loss really is.

How to interpret:

probability is the fraction of simulated paths that hit the target loss. A probability of 0.01 means a 1% chance of the target loss under the fitted normal model. avg_loss and worst_loss characterise the severity of qualifying scenarios. threshold_percentile places the target loss in the simulated distribution (e.g., 2nd percentile means the target is a 1-in-50 event).

Caveats:

The simulation assumes normally distributed returns (fitted from the historical sample). For fat-tailed assets, the true probability of extreme losses is higher than estimated here. Consider using filtered_historical_simulation from the monte_carlo sub-module for more realistic tails.

Parameters:

returns (Series | DataFrame) – Historical return series.
target_loss (float) – Target cumulative loss as a negative number (e.g. -0.20 for a 20% loss).
n_sims (int, default: 10000) – Number of Monte Carlo paths to simulate.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

"scenarios_found" – number of simulated paths that hit the target.
"probability" – estimated probability of hitting the target.
"avg_loss" – mean loss across qualifying scenarios.
"worst_loss" – worst loss observed in qualifying scenarios.
"threshold_percentile" – percentile at which the target loss sits in the simulated distribution.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0.0003, 0.01, 252))
>>> result = reverse_stress_test(returns, target_loss=-0.30, n_sims=5000, seed=42)
>>> result["probability"] >= 0
True

See also

historical_stress_test: Replay known crisis periods. stress_test_returns: User-defined additive shocks.

joint_stress_test(returns, vol_shock=2.0, spot_shock=-0.1, correlation_shock=0.0)[source]¶

Apply combined volatility, spot, and correlation shocks.

Real crises involve simultaneous increases in volatility, drops in asset prices, and spikes in correlation (diversification breaks down when you need it most). This function applies all three shocks simultaneously to produce a stressed covariance matrix and stressed expected returns.

When to use:

Use joint stress tests for: - Portfolio optimisation stress testing: feed the stressed

covariance matrix into a mean-variance optimiser.

Capital adequacy under combined adverse conditions.
Comparing diversification benefits under normal vs. stressed conditions (correlation shock toward 1.0 eliminates diversification).

Procedure:

Scale demeaned returns by vol_shock (volatility multiplier).
Shift the mean by spot_shock (additive level shift).
Blend the correlation matrix toward uniform correlation: stressed_corr = (1 - c) * corr + c * ones_matrix, where c = correlation_shock.

How to interpret:

The stressed covariance matrix (stressed_cov) reflects the combined effect of all three shocks. Pass it to wraquant.opt for stress-aware portfolio construction. Compare stressed_vol / base_vol to verify the vol scaling. Compare stressed_corr to the base correlation to see how diversification degrades.

Parameters:

returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
vol_shock (float, default: 2.0) – Volatility multiplier (e.g. 2.0 = double vol).
spot_shock (float, default: -0.1) – Additive shift to mean returns (e.g. -0.10 = subtract 10% from each asset’s mean return).
correlation_shock (float, default: 0.0) – Blend factor toward perfect correlation. 0 = unchanged, 0.5 = halfway to perfect correlation, 1 = all pairwise correlations set to 1.0.

Returns:

"stressed_mean" – stressed mean return per asset (dict).
"stressed_vol" – stressed volatility per asset (dict).
"stressed_corr" – stressed correlation matrix (ndarray).
"stressed_cov" – stressed covariance matrix (ndarray).
"base_mean" – original mean returns (dict).
"base_vol" – original volatilities (dict).

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "SPY": np.random.normal(0.0005, 0.01, 252),
...     "TLT": np.random.normal(0.0002, 0.005, 252),
... })
>>> result = joint_stress_test(returns, vol_shock=2.0, correlation_shock=0.5)
>>> result["stressed_vol"]["SPY"] > result["base_vol"]["SPY"]
True

See also

vol_stress_test: Volatility scaling only. marginal_stress_contribution: Identify worst-contributing asset.

marginal_stress_contribution(portfolio_weights, returns, scenario)[source]¶

Identify which asset contributes most to portfolio stress loss.

Decomposes the total portfolio loss under a stress scenario into per-asset contributions. This is essential for understanding where the risk is concentrated and deciding which positions to hedge or reduce.

When to use:: Use marginal stress contribution after running a stress test to answer: “which position is killing the portfolio under this scenario?” This guides targeted hedging decisions (e.g., buy puts on the worst-contributing asset).
How to interpret:: asset_contributions shows each asset’s dollar P&L under the scenario (weight_i * scenario_return_i). pct_contributions normalises these to sum to 1.0, showing the fraction of total loss attributable to each asset. worst_asset is the asset with the most negative contribution.

Parameters:

portfolio_weights (ndarray) – Weight vector aligned with returns.columns. Must have the same length as the number of columns.
returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
scenario (dict[str, float]) – Mapping of asset name to shocked return value. Assets not in the scenario use their historical mean return.

Returns:

"total_stress_loss" – portfolio return under the scenario.
"asset_contributions" – dict mapping asset name to its P&L contribution (weight * return).
"pct_contributions" – dict mapping asset name to its percentage contribution to total loss.
"worst_asset" – name of the asset contributing the most loss.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "AAPL": np.random.normal(0.001, 0.02, 100),
...     "MSFT": np.random.normal(0.001, 0.015, 100),
... })
>>> weights = np.array([0.6, 0.4])
>>> scenario = {"AAPL": -0.15, "MSFT": -0.05}
>>> result = marginal_stress_contribution(weights, returns, scenario)
>>> result["worst_asset"]
'AAPL'

See also

joint_stress_test: Generate stressed parameters for scenarios. wraquant.risk.portfolio.risk_contribution: Euler risk

decomposition (non-scenario-based).

correlation_stress(returns, shock_levels=None)[source]¶

Stress test portfolio by increasing pairwise correlations.

Blends the empirical correlation matrix toward perfect correlation at various shock levels, recomputes the covariance matrix, and measures the resulting portfolio volatility. This reveals how much diversification benefit the portfolio loses as correlations rise.

When to use:

Use correlation stress for: - Evaluating diversification fragility: how much does portfolio

risk increase if diversification breaks down?

Regulatory stress testing: correlation breakdown is a standard CCAR scenario.
Risk committee presentations: “if all correlations jump to 0.8, our portfolio vol goes from X% to Y%.”

Parameters:

returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
shock_levels (list[float] | None, default: None) – Blend factors toward perfect correlation. 0 = unchanged, 1 = all correlations set to 1.0. Defaults to [0.0, 0.25, 0.5, 0.75, 1.0].

Returns:

"results" – dict mapping shock level to a dict with "portfolio_vol" (equal-weighted portfolio volatility), "avg_correlation" (mean off-diagonal correlation), "stressed_corr" (the stressed correlation matrix).
"base_vol" – equal-weighted portfolio vol with no shock.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "A": np.random.normal(0, 0.01, 252),
...     "B": np.random.normal(0, 0.012, 252),
...     "C": np.random.normal(0, 0.008, 252),
... })
>>> result = correlation_stress(returns, shock_levels=[0.0, 0.5, 1.0])
>>> result["results"][1.0]["portfolio_vol"] >= result["base_vol"]
True

See also

joint_stress_test: Combined vol, spot, and correlation shocks. vol_stress_test: Volatility-only scaling.

liquidity_stress(returns, volumes=None, liquidity_haircuts=None, portfolio_value=1000000.0)[source]¶

Estimate liquidation cost under adverse market conditions.

Models the cost of unwinding a portfolio under stressed liquidity conditions. If volume data is provided, uses a market-impact model; otherwise, applies user-defined haircuts to each asset.

When to use:: Use liquidity stress for: - Estimating portfolio liquidation costs during crises. - Measuring liquidity-adjusted VaR (LVaR). - Satisfying regulatory requirements for liquidity stress testing

(e.g., SEC Rule 22e-4 for mutual funds).

Parameters:

returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
volumes (DataFrame | None, default: None) – Optional DataFrame of trading volumes (same shape and index as returns). If provided, liquidity cost is estimated as spread * sqrt(position / ADV).
liquidity_haircuts (dict[str, float] | None, default: None) – Optional dict mapping asset name to liquidation cost (e.g., {"AAPL": 0.001, "ILLIQ": 0.05}). If neither volumes nor haircuts are provided, uses volatility as a proxy.
portfolio_value (float, default: 1000000.0) – Total portfolio value for position sizing.

Returns:

"total_cost" – Estimated total liquidation cost ($).
"total_cost_pct" – Cost as a fraction of portfolio value.
"asset_costs" – dict mapping asset to its liquidation cost.
"days_to_liquidate" – estimated days to liquidate if limited to 10% of ADV per day (only if volumes provided).

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "A": np.random.normal(0, 0.01, 252),
...     "B": np.random.normal(0, 0.02, 252),
... })
>>> result = liquidity_stress(returns, portfolio_value=1_000_000)
>>> result["total_cost"] > 0
True

See also

vol_stress_test: Volatility scaling stress test. wraquant.execution.cost: Transaction cost modeling.

scenario_library(returns, scenarios=None)[source]¶

Apply pre-defined crisis scenarios from the built-in library.

Provides a curated set of stress scenarios calibrated to historical crises. Each scenario specifies equity shocks, volatility multipliers, and correlation shocks. The function applies each scenario to the provided returns and reports the stressed portfolio metrics.

When to use:

Use the scenario library for: - Quick stress testing without designing custom scenarios. - Regulatory reporting: standard scenarios that regulators

expect to see.

Benchmarking: compare your portfolio’s sensitivity to well-known crises.

Available scenarios:

"gfc_2008" – Global Financial Crisis
"covid_2020" – COVID-19 crash
"dot_com_2000" – Dot-com bubble burst
"rate_hike_2022" – 2022 rate hiking cycle
"stagflation" – Stagflation scenario
"flash_crash" – Flash crash (intraday)
"em_crisis" – Emerging markets crisis

Parameters:

returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
scenarios (list[str] | None, default: None) – List of scenario names from the library. If None, all scenarios are applied.

Returns:

"scenario_results" – dict mapping scenario name to a dict with "stressed_portfolio_return" (equity shock applied), "stressed_vol" (vol-scaled portfolio volatility), "scenario_params" (the raw scenario parameters).
"available_scenarios" – list of all available scenario names.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "SPY": np.random.normal(0.0005, 0.01, 252),
...     "TLT": np.random.normal(0.0002, 0.005, 252),
... })
>>> result = scenario_library(returns, scenarios=["gfc_2008", "covid_2020"])
>>> "gfc_2008" in result["scenario_results"]
True

See also

historical_stress_test: Replay actual crisis returns. joint_stress_test: Custom multi-factor stress test.

fit_gaussian_copula(returns)[source]¶

Fit a Gaussian copula to multivariate returns.

The Gaussian copula models dependence via a multivariate normal distribution applied to the rank-transformed (uniform) data. It captures linear dependence but has zero tail dependence: in the limit, extreme co-movements become independent.

Mathematical formulation:

C(u_1, …, u_d) = Phi_R(Phi^{-1}(u_1), …, Phi^{-1}(u_d))

where Phi is the standard normal CDF, Phi^{-1} is its inverse (quantile function), and Phi_R is the multivariate normal CDF with correlation matrix R.

The density is:

c(u) = |R|^{-1/2} exp(-1/2 z^T (R^{-1} - I) z)

where z_i = Phi^{-1}(u_i).

Tail dependence:

lambda_L = lambda_U = 0 (for rho < 1)

This means the Gaussian copula predicts that extreme co-movements vanish in the tails — a dangerous assumption for financial risk.

Interpretation:

The correlation matrix R describes the “normal-like” dependence structure of the data.
Key limitation: the Gaussian copula implies that joint extreme events (crashes, rallies) are asymptotically independent. The 2008 crisis demonstrated that assets crash together far more than the Gaussian copula predicts.
If tail_dependence() returns non-negligible lower tail dependence (> 0.1), the Gaussian copula is inappropriate — use fit_t_copula or fit_clayton_copula instead.
The fitted correlation matrix R is NOT the same as the Pearson correlation of raw returns. It is the correlation of the rank-transformed (copula) data.

When to use:

As a baseline for dependence modelling.
When tail dependence is genuinely absent (rare in equities).
For quick simulation of correlated returns.
Compare its AIC/BIC against Student-t to test for tail dependence.

Parameters:

returns (ndarray) – Array of shape (n_obs, n_assets) with raw returns. Each column is one asset’s return series.

Returns:

"correlation" (ndarray) – estimated copula correlation matrix R of shape (d, d). Off-diagonal values measure the normal-copula dependence. Higher values = stronger co-movement.
"copula_type" (str) – "gaussian".

Return type:

Example

>>> import numpy as np
>>> from wraquant.risk.copulas import fit_gaussian_copula, fit_t_copula
>>> rng = np.random.default_rng(42)
>>> returns = rng.multivariate_normal([0, 0], [[1, 0.6], [0.6, 1]], 500)
>>> gauss = fit_gaussian_copula(returns)
>>> print(f"Copula correlation: {gauss['correlation'][0, 1]:.3f}")
>>> # Compare with t-copula to test for tail dependence
>>> t_cop = fit_t_copula(returns)
>>> print(f"t-copula df: {t_cop['df']:.1f} (lower = heavier tails)")

See also

fit_t_copula: Student-t copula with symmetric tail dependence.: Use when df < 30 to capture crash co-movement.

fit_clayton_copula: Lower-tail dependence only (crash contagion). tail_dependence: Empirically check if tail dependence exists. copula_simulate: Generate correlated samples from the fitted copula.

fit_t_copula(returns, df=5.0)[source]¶

Fit a Student-t copula to multivariate returns.

The Student-t copula is the standard choice for equity portfolios because it captures symmetric tail dependence: the tendency of assets to experience extreme co-movements (both crashes and rallies) more often than a Gaussian model would predict.

Mathematical formulation:

C(u_1, …, u_d) = t_{R,df}(t_df^{-1}(u_1), …, t_df^{-1}(u_d))

where t_df is the univariate Student-t CDF with df degrees of freedom, and t_{R,df} is the multivariate t CDF with correlation matrix R.

Tail dependence coefficient (bivariate case):

lambda_L = lambda_U = 2 * t_{df+1}(-sqrt((df+1)(1-rho)/(1+rho)))

Unlike the Gaussian copula, this is strictly positive for finite df, meaning extreme co-movements DO occur.

Interpretation:

df (degrees of freedom) controls tail heaviness:
- df = 3-5: Very heavy tails, strong tail dependence. Typical for equity portfolios during turbulent markets. lambda ~ 0.3-0.5 for rho=0.5.
- df = 10-20: Moderate tails. Appropriate for investment- grade fixed income or diversified portfolios.
- df > 30: Nearly Gaussian (tail dependence < 0.05).
- df -> inf: Converges to Gaussian copula exactly.
Compare df across time periods: if df drops from 15 to 5, tail dependence has increased — crisis contagion is building.
Compare fitted df to Gaussian AIC: if t-copula AIC is much lower, tail dependence is statistically significant.

When to use:

Default choice for equity and credit risk modelling.
VaR/CVaR estimation for multi-asset portfolios.
When you expect both crashes AND rallies to be correlated (symmetric dependence).
If you need asymmetric tail dependence (crashes correlated but rallies independent), use fit_clayton_copula instead.

Parameters:

returns (ndarray) – Array of shape (n_obs, n_assets) with raw returns.
df (float, default: 5.0) – Degrees of freedom. Lower = heavier tails. Typical values: 3-5 for equities, 10-20 for bonds, 5-10 as a safe default. If unsure, start with 5 and compare AIC against df=10,15,20.

Returns:

"correlation" (ndarray) – Copula correlation matrix R.
"df" (float) – Degrees of freedom used.
"copula_type" (str) – "student_t".

Return type:

Example

>>> import numpy as np
>>> from wraquant.risk.copulas import fit_t_copula, tail_dependence
>>> rng = np.random.default_rng(42)
>>> # Simulate fat-tailed correlated returns
>>> returns = rng.multivariate_normal([0, 0], [[1, 0.5], [0.5, 1]], 1000)
>>> result = fit_t_copula(returns, df=5)
>>> print(f"Copula corr: {result['correlation'][0, 1]:.3f}")
>>> print(f"df={result['df']} → heavier tails than Gaussian")
>>> # Check tail dependence
>>> td = tail_dependence(returns)
>>> print(f"Lower tail dep: {td['lower']:.3f}")

Notes

Reference: Demarta, S. & McNeil, A.J. (2005). “The t Copula and Related Copulas.” International Statistical Review, 73, 111-129.

See also

fit_gaussian_copula: No tail dependence (Gaussian = t with df=inf). fit_clayton_copula: Asymmetric — lower tail dependence only. tail_dependence: Empirically estimate lambda_L and lambda_U.

fit_clayton_copula(u, v)[source]¶

Fit a Clayton copula (lower tail dependence) to bivariate data.

The Clayton copula has lower tail dependence but zero upper tail dependence. This means assets modelled with a Clayton copula tend to crash together but rally independently – a realistic pattern for equities and credit, where “risk-on/risk-off” dynamics create asymmetric co-movement.

Interpretation:

theta > 0 controls dependence strength. Higher theta = stronger lower tail dependence. theta -> 0 = independence.
lower_tail_dependence = 2^{-1/theta} is the probability that both assets are in their lower tail simultaneously, given that one is. Values above 0.3 indicate meaningful crash co-movement.
The Clayton copula is the canonical choice for modelling “contagion” and “flight to quality” dynamics.

When to use:

Two equity assets or sectors where you expect crash contagion but independent rallies.
Credit portfolios where defaults cluster.
When tail_dependence shows lower > upper.

Parameters:

u (ndarray) – First marginal uniform observations. Pass raw data and the empirical CDF is applied automatically.
v (ndarray) – Second marginal uniform observations.

Returns:

"theta" – Clayton parameter (> 0).
"copula_type" – "clayton".
"lower_tail_dependence" – 2^{-1/theta}. > 0.3 is meaningful.

Return type:

See also

fit_gumbel_copula: Upper tail dependence (opposite asymmetry). fit_t_copula: Symmetric tail dependence.

fit_gumbel_copula(u, v)[source]¶

Fit a Gumbel copula (upper tail dependence) to bivariate data.

The Gumbel copula has upper tail dependence but zero lower tail dependence. Assets rally together in extreme moves but crash independently. This is less common than the Clayton pattern in equities but can apply to certain commodity pairs (e.g., two correlated energy commodities that spike together during supply shocks).

Interpretation:

theta >= 1 controls dependence strength. theta = 1 is independence; larger theta = stronger upper tail dependence.
upper_tail_dependence = 2 - 2^{1/theta} measures the probability of joint extreme upper co-movement.

When to use:

Commodity pairs with supply-shock co-movement.
When tail_dependence shows upper > lower.
When modelling “melt-up” contagion.

Parameters:

u (ndarray) – First marginal uniform observations.
v (ndarray) – Second marginal uniform observations.

Returns:

"theta" – Gumbel parameter (>= 1).
"copula_type" – "gumbel".
"upper_tail_dependence" – 2 - 2^{1/theta}.

Return type:

See also

fit_clayton_copula: Lower tail dependence (more common in equities). fit_t_copula: Symmetric tail dependence.

fit_frank_copula(u, v)[source]¶

Fit a Frank copula (symmetric dependence, no tail dependence).

The Frank copula is symmetric with zero tail dependence in both directions. It is useful as a benchmark: if neither the Clayton nor Gumbel copula fits significantly better than Frank, then there is no evidence of asymmetric tail dependence.

Interpretation:

theta > 0: positive dependence; theta < 0: negative dependence; |theta| large = strong dependence.
The Frank copula allows negative dependence (unlike Clayton and Gumbel), making it useful for hedging pairs.
No tail dependence: extreme co-movements are modelled as asymptotically independent.

When to use:

As a benchmark against Clayton/Gumbel.
When tail_dependence shows negligible tail dependence.
For pairs with negative dependence (hedging relationships).

Parameters:

u (ndarray) – First marginal uniform observations.
v (ndarray) – Second marginal uniform observations.

Returns:

"theta" – Frank parameter. Positive = positive dependence, negative = negative dependence.
"copula_type" – "frank".

Return type:

See also

fit_gaussian_copula: Multivariate with no tail dependence. fit_clayton_copula: Lower tail dependence. fit_gumbel_copula: Upper tail dependence.

copula_simulate(copula_params, n_sims=10000, copula_type=None, seed=None)[source]¶

Simulate from a fitted copula.

Generates samples from the fitted dependence structure with uniform marginals. To get realistic return scenarios, transform each column through the inverse CDF of the desired marginal distribution (e.g., inverse normal, inverse t, or empirical quantile function).

Interpretation:

Output columns are uniform on (0, 1) – they represent the dependence structure only, not the marginal distributions.
Correlated low values in all columns = a joint crash scenario.
To convert to returns: returns[:, j] = norm.ppf(U[:, j], loc=mu_j, scale=sigma_j) or use empirical quantile functions for non-parametric marginals.
Use for Monte Carlo VaR/CVaR, portfolio simulation, or stress testing.

Parameters:

copula_params (dict[str, Any]) – Output from one of the fit_*_copula functions.
n_sims (int, default: 10000) – Number of samples to draw.
copula_type (str | None, default: None) – Override copula type (defaults to copula_params["copula_type"]).
seed (int | None, default: None) – Random seed for reproducibility.

Return type:

Returns:

Array of shape (n_sims, d) with uniform marginals in (0, 1).

Raises:

ValueError – If the copula type is not recognized.

Example

>>> from scipy.stats import norm
>>> cop = fit_gaussian_copula(returns)
>>> U = copula_simulate(cop, n_sims=10000, seed=42)
>>> # Transform to normal marginals:
>>> sim_returns = norm.ppf(U, loc=mu, scale=sigma)

tail_dependence(u, v, method='empirical', threshold=0.05)[source]¶

Estimate lower and upper tail dependence coefficients.

Tail dependence measures the probability that one variable is in its extreme tail given that the other is. This is the key quantity that copula selection hinges on.

Interpretation:

lower ~ P(V <= q | U <= q): the probability of a joint crash. Values above 0.2-0.3 are economically significant and indicate that the Gaussian copula is inappropriate.
upper ~ P(V > 1-q | U > 1-q): the probability of a joint rally.
If lower >> upper: use Clayton copula (crash contagion).
If upper >> lower: use Gumbel copula (rally contagion).
If both are similar: use Student-t copula (symmetric tails).
If both are near zero: Gaussian or Frank copula is adequate.

Caveat:

Empirical tail dependence estimates are noisy with small samples. Use threshold = 0.10 for more observations per tail (less extreme, more precise) or 0.05 for fewer observations (more extreme, noisier).

Parameters:

u (ndarray) – First marginal uniform observations.
v (ndarray) – Second marginal uniform observations.
method (str, default: 'empirical') – "empirical" (default) for non-parametric estimation.
threshold (float, default: 0.05) – Quantile threshold for tail estimation. 0.05 = 5th/95th percentile (default). 0.10 = more data in the tail estimate but less extreme.

Return type:

Returns:

Dict with "lower" and "upper" tail dependence estimates.

Example

>>> import numpy as np
>>> from scipy.stats import norm
>>> rng = np.random.default_rng(0)
>>> # Simulate from a t-copula (symmetric tail dependence)
>>> z = rng.multivariate_normal([0,0], [[1,0.5],[0.5,1]], 5000)
>>> u = norm.cdf(z[:, 0])
>>> v = norm.cdf(z[:, 1])
>>> td = tail_dependence(u, v)
>>> print(f"Lower: {td['lower']:.2f}, Upper: {td['upper']:.2f}")

rank_correlation(returns, method='both')[source]¶

Compute Kendall’s tau and/or Spearman’s rho rank correlation.

Rank correlations measure monotonic association without assuming linearity. Unlike Pearson correlation, they are invariant to monotonic transformations of the marginals and directly related to copula parameters, making them the correct correlation measure for copula modelling.

Interpretation:

Kendall’s tau: Probability of concordance minus probability of discordance. tau = 0.5 means ~75% of pairs move in the same direction.
Spearman’s rho: Pearson correlation of the ranks. Generally |rho| >= |tau| for the same data.
Both are robust to outliers (unlike Pearson).
Copula relationships: for the Gaussian copula, rho_pearson = 2*sin(pi*tau/6). For Clayton, theta = 2*tau/(1-tau).

When to use:

Always prefer rank correlation over Pearson for copula modelling.
Use Kendall’s tau for small samples (more robust).
Use Spearman’s rho for comparison with Pearson.

Parameters:

returns (ndarray) – Array of shape (n_obs, n_assets).
method (str, default: 'both') – "kendall", "spearman", or "both" (default).

Returns:

"kendall_tau" – Kendall’s tau correlation matrix.
"spearman_rho" – Spearman’s rho correlation matrix.

Return type:

Raises:

ValueError – If method is not recognized.

dcc_garch(returns, p=1, q=1)[source]¶

Fit a DCC-GARCH(p, q) model.

Currently supports p=1, q=1 (DCC(1,1)).

Procedure:

Fit univariate GARCH(1,1) to each return series.
Compute standardized residuals.
Estimate DCC parameters (a, b) via MLE.

Parameters:

returns (ndarray) – Array of shape (T, k) with k asset return series.
p (int, default: 1) – DCC lag order for innovation (currently must be 1).
q (int, default: 1) – DCC lag order for conditional (currently must be 1).

Returns:

"a" – DCC innovation parameter.
"b" – DCC persistence parameter.
"garch_params" – list of per-asset GARCH(1,1) parameter dicts.
"qbar" – unconditional correlation matrix of standardized residuals.
"conditional_vols" – array of per-asset conditional volatilities (T, k).
"std_residuals" – standardized residuals (T, k).

Return type:

rolling_correlation_dcc(returns, window=None)[source]¶

Compute DCC-based time-varying correlation matrices.

Fits DCC-GARCH to the full sample, then extracts the time-varying correlation matrix at each time step.

Parameters:

returns (ndarray) – Array of shape (T, k).
window (int | None, default: None) – Ignored (present for API compatibility with simple rolling-window methods). DCC uses the full sample.

Returns:

"correlations" – array of shape (T, k, k) with the time-varying correlation matrices.
"dcc_model" – the fitted DCC model dict.

Return type:

forecast_correlation(dcc_model, horizon=1)[source]¶

Forecast future correlation matrices from a fitted DCC model.

Uses the mean-reverting property of DCC: E[Q_{T+h}] -> Qbar as h -> inf. For finite horizons the forecasted Q is computed recursively assuming that future innovations have zero outer product (their expectation).

Parameters:

dcc_model (dict[str, Any]) – Output of dcc_garch().
horizon (int, default: 1) – Number of steps ahead to forecast.

Returns:

"forecasted_correlations" – array of shape (horizon, k, k).
"forecasted_covariances" – array of shape (horizon, k, k).

Return type:

conditional_covariance(returns, dcc_params=None)[source]¶

Compute time-varying covariance matrices from DCC-GARCH.

If dcc_params is not supplied, the model is fitted to returns automatically.

Parameters:

returns (ndarray) – Array of shape (T, k).
dcc_params (dict[str, Any] | None, default: None) – Pre-fitted DCC model dict (optional).

Returns:

"covariances" – array of shape (T, k, k) with conditional covariance matrices.
"correlations" – array of shape (T, k, k) with conditional correlation matrices.
"volatilities" – array of shape (T, k) with conditional volatilities.

Return type:

pypfopt_efficient_frontier(expected_returns, cov_matrix)[source]¶

Compute the efficient frontier using PyPortfolioOpt.

Solves for the maximum-Sharpe-ratio portfolio and returns the optimal weights along with performance metrics.

Parameters:

expected_returns (Series) – Expected annual returns per asset.
cov_matrix (DataFrame) – Annualised covariance matrix of asset returns.

Returns:

Dictionary containing:

weights – dict mapping asset names to optimal weights.
expected_return – portfolio expected annual return.
volatility – portfolio expected annual volatility.
sharpe_ratio – portfolio Sharpe ratio.

Return type:

riskfolio_portfolio(returns, method='MV')[source]¶

Optimise a portfolio using Riskfolio-Lib.

Parameters:

returns (DataFrame) – Historical return data. Each column is an asset.
method (str, default: 'MV') –
Risk measure for optimisation:
- 'MV' – minimum variance
- 'CVaR' – conditional value-at-risk
- 'MAD' – mean absolute deviation

Returns:

Dictionary containing:

weights – dict mapping asset names to optimal weights.
method – risk measure used.

Return type:

skfolio_optimize(returns, objective='min_variance')[source]¶

Optimise a portfolio using skfolio.

Parameters:

returns (DataFrame) – Historical return data. Each column is an asset.
objective (str, default: 'min_variance') –
Optimisation objective:
- 'min_variance' – minimum variance portfolio
- 'max_sharpe' – maximum Sharpe ratio portfolio

Returns:

Dictionary containing:

weights – dict mapping asset names to optimal weights.
objective – objective used.

Return type:

copulas_fit(data, copula_type='gaussian')[source]¶

Fit a copula model to multivariate data.

Uses the copulas library to fit a copula and provides methods for sampling and density evaluation.

Parameters:

data (DataFrame) – Multivariate data. Each column is a variable.
copula_type (str, default: 'gaussian') –
Type of copula to fit:
- 'gaussian' – Gaussian copula
- 'vine' – vine copula (via copulas library)

Returns:

Dictionary containing:

copula – fitted copula object.
copula_type – type of copula used.
columns – list of column names from the input data.
n_samples – number of observations used for fitting.

Return type:

copulae_fit(data, family='gaussian')[source]¶

Fit a copula using the copulae library.

Alternative to wraquant’s built-in copula fitting (copulas_fit) and the fit_*_copula functions in wraquant.risk.copulas. The copulae library provides additional families (Joe, AMH) and more robust MLE estimation via IFM (Inference Functions for Margins).

Parameters:

data (ndarray | DataFrame) – Data of shape (n, d). Can be raw observations (marginals are automatically converted to pseudo-observations) or pre-transformed uniform marginals on [0, 1].
family (str, default: 'gaussian') –
Copula family to fit:
- 'gaussian' – Gaussian copula (no tail dependence).
- 'student' – Student-t copula (symmetric tail dependence).
- 'clayton' – Clayton copula (lower tail dependence).
- 'gumbel' – Gumbel copula (upper tail dependence).
- 'frank' – Frank copula (symmetric, no tail dependence).

Returns:

Dictionary containing:

params – fitted copula parameters (structure depends on family).
log_likelihood – log-likelihood of the fitted model.
aic – Akaike information criterion.
bic – Bayesian information criterion (approximate).
fitted_copula – the fitted copula object for further use (sampling, CDF evaluation, etc.).

Return type:

Example

>>> import numpy as np
>>> from wraquant.risk.integrations import copulae_fit
>>> rng = np.random.default_rng(42)
>>> data = rng.normal(0, 1, (200, 3))
>>> result = copulae_fit(data, family="gaussian")
>>> result["log_likelihood"]

Notes

Reference: Joe (2014). Dependence Modeling with Copulas. Chapman & Hall/CRC.

See also

copulas_fit: Alternative using the copulas library.
fit_gaussian_copula: Built-in Gaussian copula implementation.
vine_copula: Vine copula fitting via pyvinecopulib.

vine_copula(data, structure='regular')[source]¶

Fit a vine copula using pyvinecopulib.

Parameters:

data (ndarray | DataFrame) – Multivariate data. Each column is a variable. Values should ideally be on the unit interval (pseudo-observations); if not, the data is rank-transformed automatically.
structure (str, default: 'regular') – Vine structure type. Currently only 'regular' is supported.

Returns:

Dictionary containing:

vinecop – fitted pyvinecopulib.Vinecop object.
structure – vine structure used.
n_vars – number of variables.
loglik – log-likelihood of the fitted model.

Return type:

extreme_value_analysis(data)[source]¶

Perform extreme value analysis using pyextremes.

Fits a Generalized Extreme Value (GEV) distribution to block maxima extracted from the data.

Parameters:

data (Series | ndarray) – Univariate time series of observations (e.g. losses or negative returns).

Returns:

Dictionary containing:

shape – GEV shape parameter (xi).
loc – GEV location parameter (mu).
scale – GEV scale parameter (sigma).
return_levels – dict of return levels for common return periods (10, 50, 100 years).

Return type:

merton_model(equity, debt, vol, rf_rate, maturity)[source]¶

Merton structural credit risk model.

Models firm equity as a European call option on total assets with strike equal to the face value of debt. The key insight: equity holders have a call option on the firm’s assets – they get the upside above the debt level but can walk away (default) if assets fall below debt.

The model iteratively solves for the unobservable asset value and asset volatility from the observable equity value and equity volatility, using the Black-Scholes option pricing relationship.

Interpretation:

distance_to_default (DD): How many standard deviations the firm’s asset value is above the default barrier. DD > 4: very safe. DD 2-4: investment grade. DD 1-2: high yield. DD < 1: distress. Moody’s KMV uses this as the primary input to their EDF (Expected Default Frequency) model.
default_probability: N(-DD), the probability that asset value drifts below debt by maturity. This is a risk-neutral probability – real-world default rates are typically lower.
asset_vol: Implied asset volatility. Higher = more default risk. Asset vol is always lower than equity vol because equity is a levered claim.
credit_spread: The yield premium investors should demand for holding risky debt. Compare to market CDS spreads to detect mispricing.

When to use:

Market-implied default probabilities from daily equity data.
Screening for distressed firms (DD < 2).
As a factor in credit scoring models.
For relative value: compare Merton-implied spreads to market CDS spreads.

Red flags:

DD < 1: firm is in acute distress.
Asset vol > 50%: inputs may be unreliable.
Equity vol is stale or missing: model won’t converge.

Parameters:

equity (float) – Current market value of equity (market cap).
debt (float) – Face value of outstanding debt (the “strike”).
vol (float) – Equity volatility (annualized, e.g., 0.30 for 30%).
rf_rate (float) – Continuous risk-free rate (annualized).
maturity (float) – Time to maturity of debt in years (typically 1).

Returns:

asset_value (float) – Implied total asset value.
asset_vol (float) – Implied asset volatility.
d1, d2 (float) – Black-Scholes d1 and d2.
distance_to_default (float) – d2, number of std devs above the default barrier.
default_probability (float) – N(-d2), risk-neutral probability of default.
credit_spread (float) – Implied credit spread over the risk-free rate (annualized).

Return type:

Example

>>> result = merton_model(equity=50e6, debt=40e6, vol=0.35,
...                       rf_rate=0.04, maturity=1.0)
>>> print(f"DD: {result['distance_to_default']:.2f}")
>>> print(f"PD: {result['default_probability']:.4f}")
>>> print(f"Spread: {result['credit_spread']*10000:.0f} bps")

See also

altman_z_score: Accounting-based bankruptcy predictor. cds_spread: Reduced-form CDS pricing.

Notes

Reference: Merton, R.C. (1974). “On the Pricing of Corporate Debt: The Risk Structure of Interest Rates.” Journal of Finance, 29(2), 449-470.

altman_z_score(working_capital, total_assets, retained_earnings, ebit, market_cap, total_liabilities, sales)[source]¶

Altman Z-Score for bankruptcy prediction.

The Z-Score combines five accounting ratios into a single discriminant score that classifies firms by financial health. Despite being from 1968, it remains one of the most widely used credit screening tools.

Interpretation:

Z > 2.99 (“safe”): Firm is financially healthy. Default probability is very low (< 1% over 2 years).
1.81 <= Z <= 2.99 (“grey zone”): Ambiguous. Firm could go either way. Warrants deeper analysis.
Z < 1.81 (“distress”): High bankruptcy risk. Historically, ~95% of firms that defaulted had Z < 1.81 one year prior.

Component interpretation:

x1 (WC/TA): Liquidity. Negative = current liabilities exceed current assets.
x2 (RE/TA): Cumulative profitability and firm age. Young firms have low retained earnings.
x3 (EBIT/TA): Operating profitability.
x4 (Market Cap/TL): Market leverage.
x5 (Sales/TA): Asset turnover efficiency.

When to use:

Quick screening of a large universe of firms.
As a factor in multi-factor credit models.
For early warning systems.

Limitations:

Designed for publicly traded manufacturing firms.
Does not capture market-implied information (use merton_model for that).
Accounting data can be manipulated.

Parameters:

working_capital (float) – Current assets minus current liabilities.
total_assets (float) – Total assets.
retained_earnings (float) – Cumulative retained earnings.
ebit (float) – Earnings before interest and taxes.
market_cap (float) – Market capitalisation of equity.
total_liabilities (float) – Total liabilities.
sales (float) – Net sales / revenue.

Returns:

z_score: The computed Z-Score.
zone: One of "safe" (Z > 2.99), "grey" (1.81 <= Z <= 2.99), or "distress" (Z < 1.81).
x1 .. x5: Individual component ratios.

Return type:

dict[str, float | str]

default_probability(rating_transitions, horizon)[source]¶

Cumulative default probability from a rating transition matrix.

Computes the probability of eventually defaulting within horizon periods, starting from each non-default rating. This is done by raising the one-period transition matrix to the horizon-th power and reading off the default column.

Interpretation:

The output is a vector where each element is the cumulative default probability for a given starting rating.
AAA will have the smallest PD; CCC the largest.
Compare to historical default rates published by Moody’s or S&P to calibrate.
Use with expected_loss for portfolio credit risk.

When to use:

Converting a rating agency transition matrix into PDs for capital calculations.
Stress testing: modify the transition matrix (increase downgrade probabilities) and re-compute PDs.

The last row/column of the transition matrix is assumed to represent the default (absorbing) state.

Parameters:

rating_transitions (ndarray) – Square transition matrix of shape (n, n) where element [i, j] is the one-period probability of migrating from rating i to rating j. Rows must sum to 1.
horizon (int) – Number of periods to compound over (e.g., 5 for 5-year cumulative PD).

Return type:

Returns:

1-D array of cumulative default probabilities for each non-default rating, length n - 1.

Example

>>> import numpy as np
>>> # Simple 3-state matrix: AAA, BBB, Default
>>> T = np.array([[0.95, 0.04, 0.01],
...               [0.02, 0.90, 0.08],
...               [0.00, 0.00, 1.00]])
>>> pd_5yr = default_probability(T, horizon=5)
>>> print(f"5yr PD from AAA: {pd_5yr[0]:.4f}")
>>> print(f"5yr PD from BBB: {pd_5yr[1]:.4f}")

credit_spread(default_prob, recovery_rate, rf_rate=0.0)[source]¶

Implied credit spread from a default probability.

Converts a default probability and recovery rate into the yield spread that compensates investors for bearing credit risk. This is the theoretical “fair value” spread – compare to market spreads to identify cheap or expensive credit.

The formula: spread = -ln(1 - PD * LGD), where LGD = 1 - R. For small PD, this simplifies to spread ~ PD * LGD.

Interpretation:

Output is annualized as a decimal (0.01 = 100 bps).
Multiply by 10,000 for basis points.
If market spread > model spread: bond is cheap (excess compensation for credit risk).
If market spread < model spread: bond is expensive or the model PD is too high.

Parameters:

default_prob (float) – Annualized probability of default (e.g., 0.02 for 2% annual PD).
recovery_rate (float) – Recovery rate in [0, 1]. Investment grade typically 0.40-0.50; high yield 0.25-0.40.
rf_rate (float, default: 0.0) – Risk-free rate (unused in simple model but accepted for API consistency).

Return type:

Returns:

Annualized credit spread as a decimal fraction. Multiply by 10,000 for basis points.

Example

>>> spread = credit_spread(0.02, 0.40)
>>> print(f"Spread: {spread*10000:.0f} bps")  # ~120 bps

loss_given_default(exposure, recovery_rate)[source]¶

Expected loss given default.

LGD = EAD * (1 - Recovery Rate). This is the dollar amount you expect to lose if the borrower defaults.

Interpretation:

A recovery rate of 0.40 means you recover 40 cents on the dollar; LGD is 60% of exposure.
Recovery rates vary by seniority: secured senior ~65%, unsecured senior ~45%, subordinated ~25%.
Use with expected_loss for Basel II/III capital calculations.

Parameters:

exposure (float) – Exposure at default (EAD) – the amount at risk.
recovery_rate (float) – Expected recovery rate in [0, 1].

Return type:

Returns:

Loss given default = exposure * (1 - recovery_rate).

expected_loss(pd_val, lgd, ead)[source]¶

Expected loss (EL = PD x LGD x EAD).

The expected loss is the central formula of credit risk management and Basel II/III regulatory capital calculation. It represents the average loss you expect from a credit exposure.

Interpretation:

EL is the mean of the loss distribution. It should be covered by pricing (loan margins, bond spreads) rather than capital reserves.
Capital reserves cover the unexpected loss (UL), which is the tail beyond EL.
For a portfolio, EL is additive: sum over all exposures.

Parameters:

pd_val (float) – Probability of default (annualized, e.g., 0.02 for 2%).
lgd (float) – Loss given default as a fraction of EAD (e.g., 0.45 for 45% loss rate).
ead (float) – Exposure at default (dollar amount at risk).

Return type:

Returns:

Expected loss in the same units as ead.

Example

>>> el = expected_loss(pd_val=0.02, lgd=0.45, ead=1_000_000)
>>> print(f"Expected loss: ${el:,.0f}")  # $9,000

cds_spread(default_intensity, recovery_rate, maturity)[source]¶

Fair CDS spread from a constant hazard rate (default intensity).

Computes the breakeven CDS premium by equating the expected protection leg (what the protection seller pays at default) with the expected premium leg (what the protection buyer pays over time).

Under a constant hazard rate model, the fair spread is approximately lambda * (1 - R), but this function uses the exact continuous-time formula with quarterly premium payments for greater accuracy.

Interpretation:

The output is the annualized spread (decimal). Multiply by 10,000 for basis points.
Compare to market CDS spreads to detect relative value.
If model spread > market spread: protection is cheap (market underestimates default risk).
If model spread < market spread: protection is expensive or there is a risk premium.
CDS spreads are approximately equal to bond spreads over swaps (CDS-bond basis ~ 0 in normal markets).

When to use:

Pricing CDS contracts given a calibrated hazard rate.
Calibrating hazard rates from market CDS spreads (invert numerically).
Converting between PD and spread for credit analysis.

Parameters:

default_intensity (float) – Constant hazard rate (lambda), annualized. E.g., 0.02 means a 2% probability of default per year.
recovery_rate (float) – Recovery rate in [0, 1]. Standard assumption is 0.40 for senior unsecured corporate debt.
maturity (float) – CDS maturity in years (standard: 1, 3, 5, 7, 10).

Return type:

Returns:

Annualized CDS spread as a decimal fraction. Multiply by 10,000 for basis points.

Example

>>> spread = cds_spread(0.02, 0.40, 5.0)
>>> print(f"5Y CDS: {spread*10000:.0f} bps")  # ~120 bps

kaplan_meier(durations, event_observed)[source]¶

Kaplan-Meier survival curve estimator.

The Kaplan-Meier (KM) estimator is the standard non-parametric method for estimating the survival function S(t) = P(T > t) from potentially censored data. “Censored” means some subjects have not yet experienced the event at the time of observation.

In finance, this answers questions like:

“What fraction of bonds survive to year 5 without defaulting?”
“How long do hedge funds typically survive before closing?”
“What is the probability a drawdown lasts longer than 60 days?”

Interpretation:

survival: S(t) is the probability of surviving beyond time t. A steep drop indicates a period of high hazard.
variance (Greenwood’s formula): Use to construct 95% confidence bands as S(t) +/- 1.96 * sqrt(variance(t)).
The median survival time is where S(t) first drops below 0.5.
A flat survival curve = low hazard rate (few events).
A curve that drops quickly early = high initial hazard (e.g., new funds failing in the first year).

Parameters:

durations (ndarray) – 1-D array of observed durations (time-to-event or time-to-censoring).
event_observed (ndarray) – 1-D boolean/int array where 1 (True) indicates the event occurred and 0 (False) indicates right-censoring.

Returns:

timeline (ndarray) – Sorted unique event times.
survival (ndarray) – Survival probability S(t) at each event time.
variance (ndarray) – Greenwood’s variance estimate for constructing confidence bands.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> durations = rng.exponential(5.0, 200)  # avg survival 5 years
>>> events = rng.binomial(1, 0.7, 200)     # 70% observed, 30% censored
>>> km = kaplan_meier(durations, events)
>>> print(f"5-year survival: {km['survival'][km['timeline'] <= 5][-1]:.2f}")

See also

nelson_aalen: Cumulative hazard estimator. log_rank_test: Compare two survival curves.

nelson_aalen(durations, event_observed)[source]¶

Nelson-Aalen cumulative hazard estimator.

Estimates the cumulative hazard function H(t), which is related to the survival function by S(t) = exp(-H(t)). While the Kaplan-Meier directly estimates S(t), the Nelson-Aalen estimator is more natural for estimating the hazard rate and for models where the hazard is the primary quantity of interest.

Interpretation:

H(t) represents the accumulated risk up to time t.
The slope of H(t) is the instantaneous hazard rate: steep segments indicate periods of high risk.
A linear H(t) suggests a constant hazard rate (exponential survival).
A concave H(t) suggests a decreasing hazard (survival gets easier over time).
A convex H(t) suggests an increasing hazard (risk accelerates – typical for aging/wear-out or credit deterioration).

Parameters:

durations (ndarray) – 1-D array of observed durations.
event_observed (ndarray) – 1-D boolean/int array indicating event occurrence.

Returns:

timeline (ndarray) – Sorted unique event times.
cumulative_hazard (ndarray) – H(t) at each time.
variance (ndarray) – Variance estimate at each time.

Return type:

See also

kaplan_meier: Direct survival function estimator. hazard_rate: Smoothed instantaneous hazard from Nelson-Aalen.

hazard_rate(durations, event_observed, bandwidth=None)[source]¶

Kernel-smoothed hazard rate estimate.

Applies Epanechnikov kernel smoothing to the Nelson-Aalen increments to produce a smooth instantaneous hazard rate function h(t).

The hazard rate answers: “Given that a subject has survived to time t, what is the instantaneous probability of the event?” This is more informative than the cumulative survival function for understanding when risk is highest.

Interpretation:

A flat hazard rate means constant risk (exponential model).
An increasing hazard means risk accelerates with time (typical for credit deterioration, infrastructure aging).
A decreasing hazard means early failures dominate and survivors become stronger (“infant mortality”).
A bathtub-shaped hazard (decreasing then increasing) is common in reliability engineering.

Parameters:

durations (ndarray) – 1-D array of observed durations.
event_observed (ndarray) – 1-D boolean/int array indicating event occurrence.
bandwidth (float | None, default: None) – Kernel bandwidth. If None, uses Silverman’s rule. Larger bandwidth = smoother curve. Smaller = more detail but noisier.

Returns:

timeline (ndarray) – Evaluation grid.
hazard (ndarray) – Smoothed hazard rate at each point.

Return type:

See also

nelson_aalen: The cumulative hazard from which this is derived. kaplan_meier: Survival function estimator.

cox_partial_likelihood(durations, event_observed, covariates, max_iter=100, tol=1e-09)[source]¶

Cox proportional hazards model via Newton-Raphson.

The Cox PH model is the workhorse of survival regression. It estimates the effect of covariates on the hazard rate without specifying the baseline hazard function (semi-parametric). The model assumes the hazard ratio is constant over time (proportional hazards assumption).

In finance: “Does leverage, profitability, or market beta affect the hazard of default, controlling for other factors?”

Interpretation:

beta[j] is the log hazard ratio for covariate j. exp(beta[j]) > 1 means the covariate increases the hazard (bad for survival). exp(beta[j]) < 1 means it decreases the hazard (protective).
se[j]: Standard error. beta[j] / se[j] gives a z-statistic. |z| > 1.96 is significant at the 5% level.
log_partial_likelihood: Higher (less negative) = better fit. Use for comparing nested models via likelihood ratio tests.

Red flags:

Very large beta (|beta| > 5): possible separation/convergence issues.
n_iter = max_iter: did not converge, results unreliable.
se contains NaN: Hessian is singular, model is degenerate.

Parameters:

durations (ndarray) – 1-D array of observed durations.
event_observed (ndarray) – 1-D boolean/int array indicating event occurrence.
covariates (ndarray) – 2-D array of shape (n_subjects, n_covariates).
max_iter (int, default: 100) – Maximum Newton-Raphson iterations.
tol (float, default: 1e-09) – Convergence tolerance on the gradient norm.

Returns:

beta (ndarray) – Regression coefficients. exp(beta) gives hazard ratios.
se (ndarray) – Standard errors of coefficients.
log_partial_likelihood (float) – Maximised log partial likelihood.
n_iter (int) – Number of iterations to convergence.

Return type:

dict[str, ndarray | float]

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> n = 200
>>> leverage = rng.uniform(0.2, 0.8, n)
>>> durations = rng.exponential(5 / (1 + leverage), n)
>>> events = np.ones(n, dtype=bool)
>>> result = cox_partial_likelihood(durations, events, leverage.reshape(-1, 1))
>>> print(f"Leverage HR: {np.exp(result['beta'][0]):.2f}")

See also

kaplan_meier: Non-parametric survival curve (no covariates). weibull_survival: Parametric survival model.

exponential_survival(lambda_param, t)[source]¶

Exponential survival function S(t) = exp(-lambda * t).

The exponential model assumes a constant hazard rate – the probability of the event in the next instant is the same regardless of how long the subject has already survived. This is the “memoryless” property.

Interpretation:

lambda = 0.1 means roughly a 10% chance of the event per unit of time.
Mean survival time = 1 / lambda.
If the hazard is actually increasing or decreasing over time, this model is too simplistic. Use weibull_survival instead.

Parameters:

lambda_param (float) – Hazard rate (constant, > 0). Higher = faster time to event.
t (float | ndarray) – Time point(s) at which to evaluate.

Return type:

Returns:

Survival probability at each t.

See also

weibull_survival: Generalises exponential with time-varying hazard.

weibull_survival(lambda_param, k, t)[source]¶

Weibull survival function S(t) = exp(-(t / lambda)^k).

The Weibull distribution generalises the exponential by allowing the hazard rate to increase or decrease over time. It is the most commonly used parametric survival model in practice.

Interpretation of the shape parameter k:

k = 1: Constant hazard (reduces to exponential). The event is equally likely at any time.
k < 1: Decreasing hazard (“burn-in”). Early failures are most common; survivors become stronger. Typical for infant mortality in manufactured goods.
k > 1: Increasing hazard (“aging”). Risk increases over time. Typical for credit deterioration in distressed firms or aging infrastructure.

In finance:

k > 1 for time-to-default: firms that have survived a long time in distress become more likely to default (debt maturity approaches, liquidity dries up).
k < 1 for drawdown recovery: if a drawdown has already lasted a long time, recovery becomes more likely (mean reversion kicks in).

Parameters:

lambda_param (float) – Scale parameter (> 0). Larger = longer survival.
k (float) – Shape parameter (> 0). k=1 is exponential. k>1 is increasing hazard. k<1 is decreasing hazard.
t (float | ndarray) – Time point(s) at which to evaluate.

Return type:

Returns:

Survival probability at each t.

See also

exponential_survival: Simplest case (k=1). kaplan_meier: Non-parametric alternative.

log_rank_test(durations1, event1, durations2, event2)[source]¶

Log-rank test comparing two survival curves.

Tests the null hypothesis that the survival functions of two groups are identical. This is the standard test for comparing survival experiences between groups.

In finance: “Do investment-grade bonds have significantly different time-to-default than high-yield bonds?” or “Do value stocks have different drawdown durations than growth stocks?”

Interpretation:

p_value < 0.05: reject H0 – the two groups have significantly different survival experiences.
observed1 >> expected1: Group 1 has more events than expected (worse survival).
observed1 << expected1: Group 1 has fewer events than expected (better survival).
The test is most powerful when the hazard ratio is constant (proportional hazards). For crossing survival curves, the Wilcoxon (Breslow) test may be more appropriate.

Parameters:

durations1 (ndarray) – Durations for group 1.
event1 (ndarray) – Event indicators for group 1.
durations2 (ndarray) – Durations for group 2.
event2 (ndarray) – Event indicators for group 2.

Returns:

test_statistic (float) – Chi-squared statistic (1 df).
p_value (float) – P-value. < 0.05 rejects equality.
observed1 (float) – Total observed events in group 1.
expected1 (float) – Expected events under H0.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> d1 = rng.exponential(5.0, 100)  # group 1: avg survival 5y
>>> d2 = rng.exponential(3.0, 100)  # group 2: avg survival 3y
>>> e1 = np.ones(100, dtype=bool)
>>> e2 = np.ones(100, dtype=bool)
>>> result = log_rank_test(d1, e1, d2, e2)
>>> print(f"p-value: {result['p_value']:.4f}")  # should be small

median_survival_time(durations, event_observed)[source]¶

Median survival time from the Kaplan-Meier estimator.

The median survival time is the smallest time t at which the estimated survival function drops to or below 0.5 – i.e., the time by which half the subjects have experienced the event.

Interpretation:

This is the “half-life” of the population.
More robust than mean survival (which is heavily influenced by censoring and long survivors).
Returns np.inf if the survival curve never reaches 0.5, which happens with heavy censoring or if more than half the subjects never experience the event.

In finance:

“The median time-to-default for CCC-rated firms is 2.3 years.”
“The median drawdown recovery time for equity portfolios is 45 trading days.”

Parameters:

durations (ndarray) – 1-D array of observed durations.
event_observed (ndarray) – 1-D boolean/int array indicating event occurrence.

Return type:

Returns:

Median survival time, or np.inf if the survival curve never drops to 0.5 (too much censoring).

importance_sampling_var(returns, n_sims=10000, target_quantile=0.01, shift=None, seed=None)[source]¶

Estimate VaR via importance sampling.

Shifts the sampling distribution toward the tail to obtain more accurate estimates of extreme quantiles with fewer simulations.

Parameters:

returns (ndarray) – 1-D array of historical returns.
n_sims (int, default: 10000) – Number of Monte Carlo draws.
target_quantile (float, default: 0.01) – Quantile level for VaR (e.g., 0.01 for 1%).
shift (float | None, default: None) – Mean shift for the importance sampling distribution. If None, automatically set to the empirical quantile.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

var: Estimated VaR (positive number = loss).
effective_sample_size: Effective sample size after reweighting, as a fraction of n_sims.

Return type:

antithetic_variates(mu, sigma, n_sims, n_assets=1, seed=None)[source]¶

Generate antithetic variate samples for variance reduction.

Produces n_sims paired samples (original + antithetic) from a normal distribution. The antithetic counterpart mirrors each draw about the mean, reducing variance for monotone functions.

Parameters:

mu (float | ndarray) – Mean(s) of the distribution. Scalar or array of length n_assets.
sigma (float | ndarray) – Standard deviation(s). Scalar or array of length n_assets.
n_sims (int) – Number of pairs to generate (total output = 2 * n_sims rows).
n_assets (int, default: 1) – Number of assets / dimensions.
seed (int | None, default: None) – Random seed for reproducibility.

Return type:

Returns:

Array of shape (2 * n_sims, n_assets) containing the original and antithetic draws interleaved.

stratified_sampling(returns, n_strata=10, n_sims=10000, seed=None)[source]¶

Stratified sampling for VaR estimation.

Divides the probability space into equal strata and draws uniformly within each stratum, ensuring better coverage of the tails.

Parameters:

returns (ndarray) – 1-D array of historical returns (used to fit a normal distribution).
n_strata (int, default: 10) – Number of strata.
n_sims (int, default: 10000) – Total number of draws (distributed evenly across strata).
seed (int | None, default: None) – Random seed for reproducibility.

Return type:

Returns:

1-D array of n_sims stratified samples drawn from the fitted normal distribution.

block_bootstrap(returns, block_size, n_sims=1000, seed=None)[source]¶

Block bootstrap for autocorrelated time series.

Resamples contiguous blocks of the input series to preserve serial dependence (Kunsch 1989, Liu & Singh 1992).

Parameters:

returns (ndarray) – 1-D array of return observations.
block_size (int) – Length of each block.
n_sims (int, default: 1000) – Number of bootstrap replications.
seed (int | None, default: None) – Random seed for reproducibility.

Return type:

Returns:

2-D array of shape (n_sims, len(returns)) where each row is one bootstrap replicate.

stationary_bootstrap(returns, avg_block_size=10.0, n_sims=1000, seed=None)[source]¶

Stationary bootstrap with random block sizes (Politis & Romano 1994).

Block lengths follow a geometric distribution with mean avg_block_size, producing a strictly stationary resampled series.

Parameters:

returns (ndarray) – 1-D array of return observations.
avg_block_size (float, default: 10.0) – Expected block size (must be > 1).
n_sims (int, default: 1000) – Number of bootstrap replications.
seed (int | None, default: None) – Random seed for reproducibility.

Return type:

Returns:

2-D array of shape (n_sims, len(returns)) where each row is one bootstrap replicate.

filtered_historical_simulation(returns, vol_model='ewma', decay=0.94, n_sims=1000, seed=None)[source]¶

Filtered historical simulation (FHS).

Combines a volatility model (EWMA or simple GARCH(1,1)) with historical bootstrap of standardised residuals to produce volatility-adjusted scenario returns.

Parameters:

returns (ndarray) – 1-D array of historical returns.
vol_model (str, default: 'ewma') – Volatility model — "ewma" (default) or "garch".
decay (float, default: 0.94) – EWMA decay factor (lambda) or, for "garch", the persistence parameter beta.
n_sims (int, default: 1000) – Number of simulated next-period returns.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

simulated_returns: 1-D array of n_sims simulated next-period returns.
current_vol: Estimated current volatility used for rescaling.
standardised_residuals: Standardised residuals from the volatility model.

Return type:

Value at Risk¶

Value-at-Risk and Conditional VaR (Expected Shortfall) estimation.

Provides the two most important tail-risk measures in quantitative risk management. VaR is a regulatory standard (Basel II/III); CVaR (Expected Shortfall) is preferred by Basel IV and is mathematically superior because it is a coherent risk measure (satisfies sub-additivity).

Both measures can be estimated via historical simulation (non-parametric) or parametric (Gaussian) methods. For heavy-tailed distributions (equities, credit), historical simulation is generally more accurate; for smooth risk surfaces (rates), parametric is often sufficient.

value_at_risk(returns, confidence=0.95, method='historical')[source]¶

Estimate Value-at-Risk (VaR).

When to use:

Mathematical formulation:

Historical: VaR_alpha = -quantile(returns, 1 - alpha) Parametric: VaR_alpha = -(mu + sigma * Phi^{-1}(1 - alpha))

where alpha is the confidence level, mu and sigma are the sample mean and standard deviation, and Phi^{-1} is the standard normal inverse CDF.

How to interpret:

Parameters:

returns (Series) – Simple return series (e.g., daily percentage changes).
confidence (float, default: 0.95) – Confidence level (e.g., 0.95 for 95%, 0.99 for 99%). Basel III uses 0.99; internal risk management often uses 0.95.
method (str, default: 'historical') –
Estimation method: - "historical" – empirical quantile (non-parametric,

default). No distributional assumption; captures fat tails.
- "parametric" – Gaussian assumption. Smooth but underestimates tail risk for leptokurtic returns.

Return type:

Returns:

VaR as a positive float representing the loss threshold. For example, 0.025 means a 2.5% loss at the given confidence level.

Raises:

ValueError – If method is not recognized.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0, 0.01, 1000))
>>> var_95 = value_at_risk(returns, confidence=0.95)
>>> var_95 > 0
True

Caveats:

VaR is not sub-additive: the VaR of a portfolio can exceed the sum of individual VaRs. Use conditional_var for a coherent measure.
Historical VaR is sensitive to the sample window; recent crises dominate short windows.
Parametric VaR severely underestimates tail risk for fat- tailed distributions (equities, credit).

See also

References

Jorion (2006), “Value at Risk: The New Benchmark”
Basel Committee on Banking Supervision (2019), “Minimum capital requirements for market risk”

conditional_var(returns, confidence=0.95, method='historical')[source]¶

Estimate Conditional VaR (Expected Shortfall / CVaR).

When to use:

CVaR is preferred over VaR for: - Portfolio optimisation (mean-CVaR optimisation is convex). - Regulatory capital under Basel IV / FRTB. - Any situation where you care about tail severity, not just

tail frequency.

Use historical CVaR with long samples (>1000 obs) and parametric CVaR when you need smooth gradients or have short data.

Mathematical formulation:

Historical: CVaR_alpha = -mean(returns | returns <= VaR_quantile) Parametric: CVaR_alpha = -(mu - sigma * phi(z_alpha) / (1 - alpha))

where z_alpha = Phi^{-1}(1 - alpha), phi is the standard normal PDF, and Phi is the CDF.

How to interpret:

Parameters:

returns (Series) – Simple return series.
confidence (float, default: 0.95) – Confidence level (e.g., 0.95 for 95%).
method (str, default: 'historical') –
Estimation method: - "historical" – mean of returns in the tail

(default). Non-parametric; captures fat tails.
- "parametric" – Gaussian formula. Smooth but underestimates tail risk for heavy-tailed distributions.

Return type:

Returns:

CVaR as a positive float representing the expected tail loss. For example, 0.035 means an expected loss of 3.5% in the tail.

Raises:

ValueError – If method is not recognized.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0, 0.01, 1000))
>>> cvar = conditional_var(returns, confidence=0.95)
>>> var = value_at_risk(returns, confidence=0.95)
>>> cvar >= var  # CVaR is always >= VaR
True

See also

value_at_risk: The VaR threshold itself. wraquant.risk.monte_carlo.importance_sampling_var: Variance-

reduced tail estimation.

References

Artzner et al. (1999), “Coherent Measures of Risk”
Rockafellar & Uryasev (2000), “Optimization of Conditional Value-at-Risk”

garch_var(returns, alpha=0.05, vol_model='GARCH', p=1, q=1, dist='normal', horizon=1)[source]¶

Value at Risk using GARCH conditional volatility.

Combines GARCH volatility forecasting with parametric VaR, producing time-varying risk estimates that adapt to current market conditions. Superior to static VaR in volatile markets.

The conditional VaR at time t is:: VaR_t = -mu + sigma_t * z_alpha

where sigma_t is the GARCH-forecasted volatility and z_alpha is the quantile of the fitted error distribution.

Parameters:

returns (Series | ndarray) – Return series.
alpha (float, default: 0.05) – Significance level (0.05 = 95% VaR).
vol_model (str, default: 'GARCH') – GARCH variant (“GARCH”, “EGARCH”, “GJR”).
p (int, default: 1) – GARCH lag order.
q (int, default: 1) – ARCH lag order.
dist (str, default: 'normal') – Error distribution (“normal”, “t”, “skewt”).
horizon (int, default: 1) – Forecast horizon in periods.

Returns:

var (pd.Series) – Time-varying VaR series.
cvar (pd.Series) – Time-varying CVaR/ES series.
conditional_vol (pd.Series) – GARCH conditional volatility.
breaches (pd.Series) – Boolean where actual loss exceeded VaR.
breach_rate (float) – Fraction of breaches (should be ~alpha).
garch_params (dict) – Fitted GARCH parameters.

Return type:

Example

>>> from wraquant.risk.var import garch_var
>>> result = garch_var(returns, alpha=0.05, vol_model="GJR", dist="t")
>>> print(f"Breach rate: {result['breach_rate']:.3f} (target: 0.050)")

greeks_var(portfolio_greeks, spot, vol, rf=0.0, dt=0.003968253968253968, spot_shock=0.01, vol_shock=0.01, n_scenarios=10000, confidence=0.95, seed=None)[source]¶

VaR approximation using portfolio Greeks (delta-gamma-vega).

The P&L approximation is:

PnL approx delta * dS + 0.5 * gamma * dS^2 + vega * d_sigma + theta * dt

where dS and d_sigma are simulated from normal distributions with standard deviations spot_shock * spot and vol_shock respectively.

When to use:: Use delta-gamma-vega VaR for options portfolios where the P&L is nonlinear in the underlying. Standard (delta-only) VaR underestimates risk for portfolios with significant gamma or vega exposure. Full revaluation VaR is more accurate but much slower; this method is a fast approximation.

Parameters:

portfolio_greeks (dict[str, float]) – Dictionary with portfolio-level Greeks. Required keys: 'delta', 'gamma'. Optional keys: 'vega', 'theta'.
spot (float) – Current spot price of the underlying.
vol (float) – Current implied volatility (annualised, e.g. 0.20 for 20%).
rf (float, default: 0.0) – Risk-free rate (annualised). Used for drift in spot dynamics.
dt (float, default: 0.003968253968253968) – Time step as a fraction of a year (default 1/252 for one trading day).
spot_shock (float, default: 0.01) – Standard deviation of the spot return used for simulation (default 0.01 = 1%). Typically set to vol * sqrt(dt) for a realistic one-day shock.
vol_shock (float, default: 0.01) – Standard deviation of the volatility change (default 0.01 = 1 vol point).
n_scenarios (int, default: 10000) – Number of Monte Carlo scenarios (default 10,000).
confidence (float, default: 0.95) – VaR confidence level (default 0.95).
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

'var' (float) – Estimated VaR (positive = loss).
'cvar' (float) – Estimated CVaR / Expected Shortfall.
'mean_pnl' (float) – Mean P&L across scenarios.
'std_pnl' (float) – Standard deviation of P&L.
'delta_component' (float) – VaR contribution from delta.
'gamma_component' (float) – VaR contribution from gamma.
'vega_component' (float) – VaR contribution from vega.
'theta_component' (float) – Deterministic theta P&L.

Return type:

Example

>>> greeks = {'delta': 100, 'gamma': -50, 'vega': 200, 'theta': -10}
>>> result = greeks_var(greeks, spot=100, vol=0.20, seed=42)
>>> result['var'] > 0
True
>>> result['cvar'] >= result['var']
True

Notes

For a single option, compute Greeks with wraquant.price.greeks and pass them here. For a portfolio, sum the Greeks across positions first.

See also

value_at_risk: Standard return-based VaR. wraquant.price.greeks: Compute option Greeks.

Metrics¶

Core risk and performance metrics.

Provides the standard risk-adjusted return ratios used across portfolio management, fund evaluation, and strategy research. Each metric quantifies a different aspect of the risk-return trade-off.

sharpe_ratio(returns, risk_free=0.0, periods_per_year=252)[source]¶

Annualized Sharpe ratio.

The Sharpe ratio measures excess return per unit of total risk (standard deviation). It is the most widely cited risk-adjusted performance measure in finance.

When to use:

Mathematical formulation:

SR = (mean(r - r_f) / std(r - r_f)) * sqrt(N)

where r is the return series, r_f is the per-period risk-free rate, and N is periods_per_year.

How to interpret:

SR < 0: strategy loses money on a risk-adjusted basis.
0 < SR < 0.5: poor; barely compensating for risk.
0.5 < SR < 1.0: acceptable for long-only strategies.
1.0 < SR < 2.0: good; typical of well-designed quant strategies.
SR > 2.0: excellent; verify this is not overfitting.
SR > 3.0: suspicious; likely backtest artifact or very short sample.

Parameters:

returns (Series) – Simple return series (e.g., daily percentage changes divided by 100).
risk_free (float, default: 0.0) – Annual risk-free rate (e.g., 0.05 for 5%).
periods_per_year (int, default: 252) – Trading periods per year (252 for daily, 12 for monthly, 52 for weekly).

Return type:

Returns:

Annualized Sharpe ratio as a float.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> daily_returns = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> sr = sharpe_ratio(daily_returns, risk_free=0.04)
>>> isinstance(sr, float)
True

Caveats:

Assumes returns are IID and normally distributed; both assumptions are violated in practice.
Penalises upside volatility equally with downside volatility; use sortino_ratio if you only care about downside risk.
Annualisation via sqrt(N) is only exact for IID returns.

See also

References

Sharpe (1966), “Mutual Fund Performance”
Bailey & Lopez de Prado (2012), “The Sharpe Ratio Efficient Frontier”

sortino_ratio(returns, risk_free=0.0, periods_per_year=252)[source]¶

Annualized Sortino ratio (downside risk only).

When to use:

Prefer Sortino over Sharpe when returns are skewed or when you only care about downside risk. Particularly useful for options strategies, trend-following, and any strategy with convex payoffs.

Mathematical formulation:

Sortino = (mean(r - r_f) / DD) * sqrt(N)

where DD = sqrt(mean(min(r - r_f, 0)^2)) is the downside deviation (computed as the second lower partial moment).

How to interpret:

Parameters:

returns (Series) – Simple return series.
risk_free (float, default: 0.0) – Annual risk-free rate.
periods_per_year (int, default: 252) – Trading periods per year.

Return type:

Returns:

Annualized Sortino ratio. Returns inf when there are no negative excess returns and the mean is positive, and 0.0 when the mean is non-positive.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> daily_returns = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> sr = sortino_ratio(daily_returns, risk_free=0.04)
>>> sr > 0
True

See also

sharpe_ratio: Uses total standard deviation. max_drawdown: Peak-to-trough loss measure.

References

Sortino & van der Meer (1991), “Downside Risk”
Sortino & Satchell (2001), “Managing Downside Risk in Financial Markets”

information_ratio(returns, benchmark)[source]¶

Information ratio (active return / tracking error).

When to use:

Mathematical formulation:

IR = mean(r_p - r_b) / std(r_p - r_b)

where r_p is the portfolio return and r_b is the benchmark return. This version is not annualized; multiply by sqrt(periods_per_year) to annualize.

How to interpret:

IR < 0: underperforming the benchmark.
0 < IR < 0.3: modest skill; hard to distinguish from luck.
0.3 < IR < 0.5: good active management.
IR > 0.5: exceptional; sustained alpha generation.
IR > 1.0: very rare and likely indicates short sample bias.

Parameters:

returns (Series) – Portfolio return series.
benchmark (Series) – Benchmark return series (same frequency and index).

Return type:

Returns:

Information ratio (not annualized) as a float.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> portfolio = pd.Series(np.random.normal(0.0006, 0.01, 252))
>>> benchmark = pd.Series(np.random.normal(0.0004, 0.009, 252))
>>> ir = information_ratio(portfolio, benchmark)
>>> isinstance(ir, float)
True

See also

sharpe_ratio: Absolute risk-adjusted return. hit_ratio: Fraction of positive-return periods.

References

Grinold & Kahn (2000), “Active Portfolio Management”

max_drawdown(prices)[source]¶

Maximum drawdown from a price series.

Maximum drawdown is the largest peak-to-trough decline in the price series. It measures the worst historical loss an investor would have experienced if they bought at the peak and sold at the trough.

When to use:

Mathematical formulation:

MDD = min_t ( (P_t - max_{s<=t} P_s) / max_{s<=t} P_s )

The result is a negative number (or zero if the price series only goes up).

How to interpret:

MDD = -0.10 means the worst peak-to-trough loss was 10%.
MDD = -0.50 means the strategy lost half its value at worst.
For S&P 500 since 1928, the worst MDD was about -0.86 (1929-1932). Post-2000, the worst was about -0.57 (GFC).
A strategy with a Sharpe of 1.0 and max drawdown of -0.30 has a Calmar ratio of about 0.33.

Parameters:: prices (Series) – Price or equity curve series (not returns). Must be positive values.
Return type:: float
Returns:: Maximum drawdown as a negative float (e.g., -0.25 for a 25% drawdown). Returns 0.0 if the series is monotonically increasing.

Example

>>> import pandas as pd
>>> prices = pd.Series([100, 110, 105, 95, 108, 102])
>>> mdd = max_drawdown(prices)
>>> round(mdd, 4)
-0.1364

See also

References

Magdon-Ismail & Atiya (2004), “Maximum Drawdown”

hit_ratio(returns)[source]¶

Fraction of positive return periods (win rate).

When to use:

How to interpret:

0.50 = coin flip; no directional edge.
0.55 = statistically meaningful edge on daily data.
0.60+ = strong edge; verify you are not overfitting.
< 0.40 does not mean a bad strategy if the avg win >> avg loss.

Parameters:: returns (Series) – Simple return series.
Return type:: float
Returns:: Hit ratio as a float between 0 and 1.

Example

>>> import pandas as pd
>>> returns = pd.Series([0.01, -0.005, 0.008, -0.003, 0.012])
>>> hit_ratio(returns)
0.6

See also

sharpe_ratio: Risk-adjusted return measure. information_ratio: Alpha per unit of tracking error.

treynor_ratio(returns, benchmark, risk_free=0.0, periods_per_year=252)[source]¶

Treynor ratio: excess return per unit of systematic (beta) risk.

When to use:

Use Treynor for comparing portfolios that are part of a larger diversified portfolio (so only systematic risk matters). If the portfolio is the investor’s entire wealth, use Sharpe instead.

Mathematical formulation:

Treynor = (R_p - R_f) / beta_p

where R_p is annualized portfolio return, R_f is the risk-free rate, and beta_p is the portfolio’s beta vs the benchmark.

Parameters:

returns (Series) – Portfolio return series.
benchmark (Series) – Benchmark return series (same frequency and index).
risk_free (float, default: 0.0) – Annual risk-free rate.
periods_per_year (int, default: 252) – Trading periods per year.

Return type:

Returns:

Treynor ratio as a float. Higher is better. Negative indicates the portfolio underperformed the risk-free rate.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> portfolio = 1.2 * market + np.random.normal(0.0001, 0.005, 252)
>>> tr = treynor_ratio(portfolio, market, risk_free=0.04)
>>> isinstance(tr, float)
True

See also

sharpe_ratio: Total risk-adjusted return. jensens_alpha: Excess return above CAPM prediction.

References

Treynor (1965), “How to Rate Management of Investment Funds”

m_squared(returns, benchmark, risk_free=0.0, periods_per_year=252)[source]¶

M-squared (Modigliani-Modigliani) risk-adjusted performance.

When to use:: Use M-squared when you want to compare two portfolios with different risk levels on a common scale. M-squared answers: “if both portfolios had the same risk as the benchmark, which would earn more?”
Mathematical formulation:: M^2 = SR_p * sigma_b + R_f

Parameters:

returns (Series) – Portfolio return series.
benchmark (Series) – Benchmark return series.
risk_free (float, default: 0.0) – Annual risk-free rate.
periods_per_year (int, default: 252) – Trading periods per year.

Return type:

Returns:

M-squared as an annualized return (float). Positive means the portfolio outperforms the benchmark on a risk-adjusted basis.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> portfolio = pd.Series(np.random.normal(0.0006, 0.01, 252))
>>> benchmark = pd.Series(np.random.normal(0.0004, 0.009, 252))
>>> m2 = m_squared(portfolio, benchmark, risk_free=0.04)
>>> isinstance(m2, float)
True

See also

sharpe_ratio: Dimensionless risk-adjusted return.

References

Modigliani & Modigliani (1997), “Risk-Adjusted Performance”

jensens_alpha(returns, benchmark, risk_free=0.0, periods_per_year=252)[source]¶

Jensen’s alpha: excess return above CAPM-predicted return.

Mathematical formulation:

alpha = R_p - [R_f + beta * (R_m - R_f)]

where R_p is the portfolio return, R_m is the benchmark return, and beta is the portfolio’s beta to the benchmark.

Parameters:

returns (Series) – Portfolio return series.
benchmark (Series) – Benchmark return series.
risk_free (float, default: 0.0) – Annual risk-free rate.
periods_per_year (int, default: 252) – Trading periods per year.

Return type:

Returns:

Annualized Jensen’s alpha as a float.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> stock = 0.0002 + 1.0 * market + np.random.normal(0, 0.005, 252)
>>> alpha = jensens_alpha(stock, market, risk_free=0.04)
>>> isinstance(alpha, float)
True

See also

treynor_ratio: Risk-adjusted return using beta. appraisal_ratio: Alpha per unit of residual risk.

References

Jensen (1968), “The Performance of Mutual Funds in the Period 1945-1964”

appraisal_ratio(returns, benchmark, risk_free=0.0, periods_per_year=252)[source]¶

Appraisal ratio: Jensen’s alpha per unit of residual risk.

Mathematical formulation:

AR = alpha / sigma_epsilon

where alpha is Jensen’s alpha and sigma_epsilon is the standard deviation of the regression residuals (annualized).

Parameters:

returns (Series) – Portfolio return series.
benchmark (Series) – Benchmark return series.
risk_free (float, default: 0.0) – Annual risk-free rate.
periods_per_year (int, default: 252) – Trading periods per year.

Return type:

Returns:

Appraisal ratio as a float. Higher is better.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> stock = 0.0003 + 1.0 * market + np.random.normal(0, 0.005, 252)
>>> ar = appraisal_ratio(stock, market, risk_free=0.04)
>>> isinstance(ar, float)
True

See also

jensens_alpha: The numerator of the appraisal ratio. information_ratio: Alpha per unit of tracking error.

References

Treynor & Black (1973), “How to Use Security Analysis to Improve Portfolio Selection”

capture_ratios(returns, benchmark)[source]¶

Up-capture and down-capture ratios.

Capture ratios measure how much of the benchmark’s up and down movements the portfolio captures. An ideal portfolio has high up-capture (>100%) and low down-capture (<100%).

When to use:

Use capture ratios for: - Evaluating defensive vs aggressive positioning: a portfolio

with 90% up-capture and 70% down-capture is defensively positioned and will outperform in bear markets.

Manager selection: compare capture ratios across funds.
Style analysis: growth funds typically have high up-capture and high down-capture; value funds often have lower both.

Mathematical formulation:

Up-capture = (mean(r_p | r_b > 0) / mean(r_b | r_b > 0)) * 100 Down-capture = (mean(r_p | r_b < 0) / mean(r_b | r_b < 0)) * 100

Capture ratio = up-capture / down-capture

Parameters:

returns (Series) – Portfolio return series.
benchmark (Series) – Benchmark return series (same frequency and index).

Returns:

up_capture (float) – Percentage of benchmark’s up movements captured (>100 = amplified).
down_capture (float) – Percentage of benchmark’s down movements captured (<100 = dampened losses).
capture_ratio (float) – up_capture / down_capture. Values > 1 indicate the portfolio adds value through asymmetric participation.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> benchmark = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> portfolio = 0.8 * benchmark + np.random.normal(0.0001, 0.005, 252)
>>> caps = capture_ratios(portfolio, benchmark)
>>> caps["up_capture"] > 0
True

See also

treynor_ratio: Systematic risk-adjusted return.

Beta Estimation¶

Beta estimation models for systematic risk measurement.

Beta measures the sensitivity of an asset’s returns to a benchmark (typically the market). A beta of 1.0 means the asset moves in lockstep with the benchmark; >1.0 means amplified moves; <1.0 means dampened moves; <0 means the asset moves inversely.

This module provides six beta estimators, each suited to different situations:

Rolling OLS beta (rolling_beta) – the workhorse; shows how beta evolves over time. Use for regime analysis and dynamic hedging.
Blume adjustment (blume_adjusted_beta) – shrinks raw beta toward 1.0 using the empirical regression-to-mean relationship.
Vasicek adjustment (vasicek_adjusted_beta) – Bayesian shrinkage that incorporates estimation uncertainty.
Dimson beta (dimson_beta) – sums lagged betas to correct for non-synchronous trading in illiquid assets.
Conditional beta (conditional_beta) – separate up-market and down-market betas to capture asymmetric sensitivity.
EWMA beta (ewma_beta) – exponentially weighted beta that adapts quickly to recent market conditions.

References

Blume (1971), “On the Assessment of Risk”
Vasicek (1973), “A Note on Using Cross-Sectional Information in Bayesian Estimation of Security Betas”
Dimson (1979), “Risk Measurement When Shares are Subject to Infrequent Trading”

rolling_beta(returns, benchmark, window=60)[source]¶

Rolling OLS beta of asset returns against a benchmark.

When to use:

Use rolling beta to: - Monitor regime changes in market exposure (beta rising during

sell-offs indicates contagion).

Calibrate dynamic hedging ratios (e.g., beta-hedge a long position with index futures).
Detect structural breaks in a strategy’s factor exposure.

Mathematical formulation:

beta_t = Cov(r, b; t-w:t) / Var(b; t-w:t)

where r is the asset return, b is the benchmark return, and w is the rolling window size.

Parameters:

returns (Series) – Asset return series (e.g., daily simple returns).
benchmark (Series) – Benchmark return series (same frequency and aligned index). Typically a broad market index (S&P 500, MSCI World).
window (int, default: 60) – Rolling window size in periods. 60 trading days (~3 months) is standard for equity beta. Use 120-252 for more stable estimates; use 20-40 for faster-reacting estimates.

Return type:

Returns:

pd.Series of rolling beta values, indexed to match returns. The first window - 1 values are NaN (insufficient data).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> stock = 1.2 * market + np.random.normal(0, 0.005, 252)
>>> beta = rolling_beta(stock, market, window=60)
>>> abs(beta.iloc[-1] - 1.2) < 0.3
True

See also

ewma_beta: Exponentially weighted alternative (no fixed window). conditional_beta: Separate up/down market betas.

References

Ang & Chen (2007), “Asymmetric Correlations of Equity Portfolios”

blume_adjusted_beta(raw_beta)[source]¶

Blume-adjusted beta (mean-reversion adjustment).

Blume (1971) documented that betas regress toward 1.0 over time. The adjustment applies the empirical relationship:

adjusted_beta = 0.33 + 0.67 * raw_beta

This is the adjustment used by Bloomberg and most commercial risk systems when reporting “adjusted beta.”

When to use:: Use Blume adjustment for forward-looking beta estimates (e.g., cost of equity in CAPM). The raw historical beta is a biased predictor of future beta; the Blume adjustment reduces this bias.

Parameters:: raw_beta (float) – Historical OLS beta estimate.
Return type:: float
Returns:: Adjusted beta as a float, shrunk toward 1.0.

Example

>>> blume_adjusted_beta(1.5)
1.335
>>> blume_adjusted_beta(0.5)
0.665
>>> blume_adjusted_beta(1.0)
1.0

See also

vasicek_adjusted_beta: Bayesian shrinkage with uncertainty. rolling_beta: Source of the raw beta input.

References

Blume (1971), “On the Assessment of Risk”, Journal of Finance
Blume (1975), “Betas and Their Regression Tendencies”

vasicek_adjusted_beta(raw_beta, cross_sectional_mean=1.0, raw_se=0.2, prior_se=0.3)[source]¶

Vasicek Bayesian shrinkage beta adjustment.

When to use:

Mathematical formulation:

adjusted_beta = (prior_se^2 / (prior_se^2 + raw_se^2)) * raw_beta

(raw_se^2 / (prior_se^2 + raw_se^2)) * cross_sectional_mean

Parameters:

raw_beta (float) – Historical OLS beta estimate.
cross_sectional_mean (float, default: 1.0) – Prior mean beta (cross-sectional average). Typically 1.0 for market beta.
raw_se (float, default: 0.2) – Standard error of the raw beta estimate from OLS regression. Higher values cause more shrinkage.
prior_se (float, default: 0.3) – Standard deviation of the cross-sectional beta distribution. Represents uncertainty in the prior.

Return type:

Returns:

Vasicek-adjusted beta as a float.

Example

>>> vasicek_adjusted_beta(1.5, cross_sectional_mean=1.0, raw_se=0.2, prior_se=0.3)
1.3461538461538463
>>> # High SE -> more shrinkage toward prior
>>> vasicek_adjusted_beta(1.5, raw_se=0.5, prior_se=0.3)
1.1323529411764706

See also

blume_adjusted_beta: Simpler fixed-weight adjustment.

References

Vasicek (1973), “A Note on Using Cross-Sectional Information in Bayesian Estimation of Security Betas”

dimson_beta(returns, benchmark, lags=1)[source]¶

Dimson beta for illiquid or thinly traded assets.

When to use:

Mathematical formulation:

r_t = alpha + beta_0 * b_t + beta_1 * b_{t-1} + … + beta_k * b_{t-k} + eps

Dimson beta = sum(beta_0, beta_1, …, beta_k)

Parameters:

returns (Series) – Asset return series.
benchmark (Series) – Benchmark return series (same frequency and index).
lags (int, default: 1) – Number of lagged benchmark terms to include. 1 is standard for daily data; use 2-3 for very illiquid assets.

Returns:

total_beta (float) – Sum of all lag coefficients (the Dimson-adjusted beta).
lag_betas (list[float]) – Individual coefficients for each lag (index 0 = contemporaneous).
alpha (float) – Regression intercept.
r_squared (float) – R-squared of the multiple regression.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> benchmark = pd.Series(np.random.normal(0.0005, 0.01, 300))
>>> # Illiquid asset: reacts with a lag
>>> returns = 0.5 * benchmark + 0.4 * benchmark.shift(1).fillna(0) + \\
...     np.random.normal(0, 0.005, 300)
>>> result = dimson_beta(returns, benchmark, lags=1)
>>> result["total_beta"] > result["lag_betas"][0]
True

See also

rolling_beta: Standard OLS beta (assumes synchronous trading). ewma_beta: Exponentially weighted beta.

References

Dimson (1979), “Risk Measurement When Shares are Subject to Infrequent Trading”, Journal of Financial Economics

conditional_beta(returns, benchmark)[source]¶

Conditional (asymmetric) beta: separate up-market and down-market betas.

When to use:

Use conditional beta to: - Assess downside protection: an asset with low downside beta

and high upside beta is a desirable portfolio component.

Detect asymmetric risk exposure: if downside_beta >> upside_beta, the asset amplifies losses more than gains.
Evaluate hedge fund or options-like payoff profiles.

Mathematical formulation:

Up-market: r_t = alpha_up + beta_up * b_t + eps, for b_t > 0 Down-market: r_t = alpha_down + beta_down * b_t + eps, for b_t <= 0

Parameters:

returns (Series) – Asset return series.
benchmark (Series) – Benchmark return series (same frequency and index).

Returns:

upside_beta (float) – Beta in up-market periods.
downside_beta (float) – Beta in down-market periods.
beta_asymmetry (float) – downside_beta - upside_beta. Positive means the asset is more sensitive to down moves.
n_up (int) – Number of up-market observations.
n_down (int) – Number of down-market observations.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = pd.Series(np.random.normal(0.0005, 0.01, 500))
>>> stock = 1.0 * market + np.random.normal(0, 0.005, 500)
>>> result = conditional_beta(stock, market)
>>> isinstance(result["upside_beta"], float)
True

See also

rolling_beta: Time-varying beta (not conditional on direction). ewma_beta: Exponentially weighted beta.

References

Pettengill, Sundaram & Mathur (1995), “The Conditional Relation Between Beta and Returns”
Ang & Chen (2002), “Asymmetric Correlations of Equity Portfolios”

ewma_beta(returns, benchmark, halflife=60)[source]¶

Exponentially weighted moving average (EWMA) beta.

When to use:

Use EWMA beta when you need a smooth, responsive beta estimate that adapts quickly to regime changes. Compared to rolling beta: - EWMA has no “cliff effect” (old observations do not drop out

abruptly).

EWMA adapts faster to structural breaks (smaller halflife).
EWMA is smoother (no window-edge artifacts).

Mathematical formulation:

beta_t = EWCov(r, b; lambda) / EWVar(b; lambda)

where lambda = 1 - exp(-ln(2) / halflife) is the decay factor.

Parameters:

returns (Series) – Asset return series.
benchmark (Series) – Benchmark return series (same frequency and index).
halflife (int, default: 60) – Decay halflife in periods. 60 days is standard. Shorter halflife (20-30) reacts faster but is noisier. Longer halflife (90-120) is smoother but lags.

Return type:

Returns:

pd.Series of EWMA beta values.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> stock = 1.3 * market + np.random.normal(0, 0.005, 252)
>>> beta = ewma_beta(stock, market, halflife=60)
>>> abs(beta.iloc[-1] - 1.3) < 0.4
True

See also

rolling_beta: Fixed-window alternative. conditional_beta: Direction-dependent beta.

References

RiskMetrics Technical Document (1996), J.P. Morgan

Factor Risk¶

Factor risk models for return attribution and risk decomposition.

Factor models decompose portfolio risk into systematic (factor-driven) and idiosyncratic (asset-specific) components. This is fundamental to understanding where portfolio risk comes from and whether factor exposures are intentional or accidental.

This module provides four approaches:

Fundamental factor model (factor_risk_model) – regress returns on user-supplied factors (e.g., Fama-French, macro factors). Use when you know which factors matter.
Statistical factor model (statistical_factor_model) – extract latent factors via PCA. Use when you do not have a prior on which factors drive returns.
Fama-French regression (fama_french_regression) – specialised for the classic Fama-French framework with named factors (MKT, SMB, HML, etc.).
Factor contribution (factor_contribution) – given portfolio weights and factor exposures, decompose portfolio risk into factor contributions.

References

Fama & French (1993), “Common Risk Factors in the Returns on Stocks and Bonds”
Connor & Korajczyk (1986), “Performance Measurement with the Arbitrage Pricing Theory”
Menchero (2011), “The Barra Risk Model Handbook”

factor_risk_model(returns, factors)[source]¶

Regress asset returns on factors and decompose total risk.

When to use:

Mathematical formulation:

r_t = alpha + B * f_t + eps_t

Total variance = B’ * Sigma_f * B + sigma_eps^2 Factor risk share = B’ * Sigma_f * B / Total variance Specific risk share = sigma_eps^2 / Total variance

Parameters:

returns (Series | DataFrame) – Asset return series (pd.Series for one asset, pd.DataFrame for multiple assets – uses first column).
factors (DataFrame) – DataFrame of factor returns with columns as factor names and a compatible index.

Returns:

betas (dict[str, float]) – Factor loadings (regression coefficients). Positive beta = positive exposure.
alpha (float) – Regression intercept (excess return not explained by factors).
factor_risk (float) – Fraction of total variance explained by factors (0 to 1).
specific_risk (float) – Fraction of total variance from idiosyncratic sources (1 - factor_risk).
r_squared (float) – R-squared of the regression.
residual_vol (float) – Annualized volatility of residuals.
contributions (dict[str, float]) – Each factor’s individual contribution to systematic variance.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> mkt = np.random.normal(0.0005, 0.01, 252)
>>> smb = np.random.normal(0, 0.005, 252)
>>> stock = 1.1 * mkt + 0.3 * smb + np.random.normal(0, 0.005, 252)
>>> result = factor_risk_model(
...     pd.Series(stock),
...     pd.DataFrame({"MKT": mkt, "SMB": smb}),
... )
>>> result["factor_risk"] > 0.5
True

See also

statistical_factor_model: PCA-based (no prior on factors). fama_french_regression: Specialised Fama-French interface.

References

Menchero (2011), “The Barra Risk Model Handbook”

statistical_factor_model(returns, n_factors=3)[source]¶

PCA-based statistical factor model with risk decomposition.

When to use:: Use statistical factor models when you do not have a prior on which factors drive returns. PCA discovers the dominant sources of covariation. Useful for: - Constructing factor-mimicking portfolios. - Dimensionality reduction before portfolio optimisation. - Identifying hidden risk concentrations.

Parameters:

returns (DataFrame) – DataFrame of asset returns (columns = assets, rows = dates). Should have at least n_factors + 1 columns.
n_factors (int, default: 3) – Number of principal components to extract. 3-5 is typical for equity portfolios.

Returns:

factors (pd.DataFrame) – Extracted factor return series (columns: PC1, PC2, …).
loadings (np.ndarray) – Factor loadings matrix (n_assets x n_factors).
explained_variance (np.ndarray) – Variance explained by each factor.
explained_variance_ratio (np.ndarray) – Fraction of total variance explained by each factor.
cumulative_variance_ratio (np.ndarray) – Cumulative fraction of variance explained.
factor_risk (float) – Total fraction of variance explained by all extracted factors.
specific_risk (float) – Fraction of variance not explained.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> market = np.random.normal(0, 0.01, 252)
>>> returns = pd.DataFrame({
...     f"asset_{i}": market * (0.5 + i * 0.2) + np.random.normal(0, 0.005, 252)
...     for i in range(5)
... })
>>> result = statistical_factor_model(returns, n_factors=2)
>>> result["factor_risk"] > 0.3
True

See also

factor_risk_model: When you know which factors to use.

References

Connor & Korajczyk (1986), “Performance Measurement with the Arbitrage Pricing Theory”

fama_french_regression(returns, factors_df)[source]¶

Fama-French factor regression with full diagnostics.

When to use:: Use for performance attribution and alpha measurement. The classic 3-factor model (Mkt, SMB, HML) is the minimum; the 5-factor model adds RMW (profitability) and CMA (investment). Add Mom (momentum) for the 6-factor model.

Parameters:

returns (Series) – Asset or portfolio return series (excess of risk-free rate if the factors are excess returns).
factors_df (DataFrame) – DataFrame of factor returns. Column names should be descriptive (e.g., “Mkt-RF”, “SMB”, “HML”).

Returns:

alpha (float) – Jensen’s alpha (intercept).
betas (dict[str, float]) – Factor loadings.
t_stats (dict[str, float]) – t-statistics for each coefficient (including alpha under key “alpha”).
p_values (dict[str, float]) – p-values for each coefficient.
r_squared (float) – R-squared.
adj_r_squared (float) – Adjusted R-squared.
residual_vol (float) – Annualized residual volatility.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> mkt = np.random.normal(0.0005, 0.01, 252)
>>> smb = np.random.normal(0, 0.005, 252)
>>> hml = np.random.normal(0, 0.005, 252)
>>> stock = 0.0001 + 1.1 * mkt + 0.3 * smb - 0.2 * hml + \\
...     np.random.normal(0, 0.003, 252)
>>> factors = pd.DataFrame({"Mkt-RF": mkt, "SMB": smb, "HML": hml})
>>> result = fama_french_regression(pd.Series(stock), factors)
>>> abs(result["betas"]["Mkt-RF"] - 1.1) < 0.2
True

See also

factor_risk_model: General factor regression with risk decomposition.

References

Fama & French (1993), “Common Risk Factors in the Returns on Stocks and Bonds”
Fama & French (2015), “A Five-Factor Asset Pricing Model”

factor_contribution(weights, factor_betas, factor_cov)[source]¶

Decompose portfolio factor risk into per-factor contributions.

Given portfolio weights, a matrix of factor loadings, and the factor covariance matrix, computes how much each factor contributes to total portfolio factor risk (variance).

When to use:

Mathematical formulation:

Portfolio factor variance = w’ * B * Sigma_f * B’ * w

Factor i contribution = w’ * B_i * (Sigma_f * B’ * w)_i / total_var

Parameters:

weights (ndarray) – Portfolio weight vector (n_assets,).
factor_betas (ndarray) – Factor loading matrix (n_assets x n_factors). Each row is an asset’s factor exposures.
factor_cov (ndarray) – Factor covariance matrix (n_factors x n_factors).

Returns:

total_factor_var (float) – Total portfolio factor variance.
total_factor_vol (float) – Square root of factor variance.
factor_contributions (np.ndarray) – Each factor’s variance contribution (sums to total_factor_var).
factor_pct_contributions (np.ndarray) – Percentage contributions (sum to 1.0).

Return type:

Example

>>> import numpy as np
>>> weights = np.array([0.3, 0.3, 0.4])
>>> betas = np.array([[1.0, 0.5], [1.2, -0.3], [0.8, 0.1]])
>>> factor_cov = np.array([[0.0004, 0.00005], [0.00005, 0.0001]])
>>> result = factor_contribution(weights, betas, factor_cov)
>>> result["total_factor_var"] > 0
True

See also

factor_risk_model: Estimate factor betas from return data. statistical_factor_model: Extract latent factors via PCA.

Portfolio Analytics¶

Advanced portfolio risk analytics.

Extends the basic portfolio risk tools (volatility, risk contribution, diversification ratio) with VaR decomposition, risk budgeting, and benchmark-relative analytics. These functions are essential for institutional portfolio management, risk budgeting, and performance attribution.

Key concepts:

Component VaR: how much each asset contributes to portfolio VaR.
Marginal VaR: sensitivity of portfolio VaR to a small change in weight (used for position sizing).
Incremental VaR: change in portfolio VaR from adding/removing an asset entirely.
Risk budgeting: find weights that produce equal (or target) risk contributions.
Tracking error: active risk relative to a benchmark.
Active share: how different the portfolio is from the benchmark.

References

Litterman (1996), “Hot Spots and Hedges” (Euler decomposition)
Maillard, Roncalli & Teiletche (2010), “The Properties of Equally Weighted Risk Contribution Portfolios”
Cremers & Petajisto (2009), “How Active Is Your Fund Manager?”

component_var(weights, returns, alpha=0.05)[source]¶

Component Value-at-Risk: per-asset contribution to portfolio VaR.

When to use:

Use component VaR for: - Identifying which assets dominate portfolio tail risk. - Setting per-asset risk limits. - Reporting risk contributions to portfolio managers and risk

committees.

Mathematical formulation:

Component VaR_i = w_i * (partial VaR / partial w_i)

Under the delta-normal approximation: CVaR_i = w_i * (Sigma @ w)_i / sigma_p * VaR_p

Parameters:

weights (ndarray) – Portfolio weight vector (n_assets,).
returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
alpha (float, default: 0.05) – Significance level (0.05 = 95% VaR).

Return type:

Returns:

pd.Series of per-asset VaR contributions, indexed by asset names. Sum equals the portfolio VaR.

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "A": np.random.normal(0.0005, 0.01, 252),
...     "B": np.random.normal(0.0003, 0.015, 252),
... })
>>> weights = np.array([0.6, 0.4])
>>> cvar = component_var(weights, returns, alpha=0.05)
>>> cvar.sum() > 0  # total VaR is positive
True

See also

marginal_var: Sensitivity of VaR to weight changes. incremental_var: VaR change from adding/removing an asset.

marginal_var(weights, cov, alpha=0.05)[source]¶

Marginal VaR: sensitivity of portfolio VaR to weight changes.

Marginal VaR measures how much portfolio VaR changes for a small (infinitesimal) change in the weight of each asset. It is the gradient of portfolio VaR with respect to weights.

When to use:

Use marginal VaR for: - Position sizing: assets with high marginal VaR should have

smaller positions.

Optimisation: marginal VaR should be equal across assets at the optimal portfolio (risk parity condition).
Hedging: the hedge ratio is proportional to the marginal VaR.

Mathematical formulation:

Marginal VaR_i = dVaR/dw_i = z_alpha * (Sigma @ w)_i / sigma_p

Parameters:

weights (ndarray) – Portfolio weight vector (n_assets,).
cov (ndarray) – Covariance matrix (n_assets x n_assets).
alpha (float, default: 0.05) – Significance level (0.05 = 95% VaR).

Return type:

Returns:

np.ndarray of marginal VaR values per asset.

Example

>>> import numpy as np
>>> cov = np.array([[0.0004, 0.0001], [0.0001, 0.0009]])
>>> weights = np.array([0.6, 0.4])
>>> mvar = marginal_var(weights, cov, alpha=0.05)
>>> len(mvar) == 2
True

See also

component_var: Additive VaR decomposition (weight * marginal VaR).

incremental_var(weights, returns, alpha=0.05)[source]¶

Incremental VaR: change in portfolio VaR from adding each asset.

When to use:: Use incremental VaR when deciding whether to add or remove a position. Unlike marginal VaR (which is an infinitesimal measure), incremental VaR captures the full nonlinear impact including diversification effects.

Parameters:

weights (ndarray) – Portfolio weight vector (n_assets,).
returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
alpha (float, default: 0.05) – Significance level (0.05 = 95% VaR).

Return type:

Returns:

np.ndarray of incremental VaR values per asset. Positive means adding the asset increases portfolio VaR (adds risk); negative means it reduces VaR (diversification benefit).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "A": np.random.normal(0.001, 0.01, 252),
...     "B": np.random.normal(0.0005, 0.008, 252),
... })
>>> weights = np.array([0.6, 0.4])
>>> ivar = incremental_var(weights, returns, alpha=0.05)
>>> len(ivar) == 2
True

See also

component_var: Euler-based additive decomposition. marginal_var: Infinitesimal sensitivity.

risk_budgeting(cov, target_risk=None)[source]¶

Find portfolio weights that achieve target risk contributions.

Risk budgeting finds the weights such that each asset’s risk contribution (Euler decomposition) matches a target budget. With equal targets (default), this is the risk parity portfolio.

When to use:

Use risk budgeting for: - Risk parity portfolio construction (equal risk contribution). - Custom risk allocation (e.g., 60% risk from equities, 40%

from bonds, regardless of capital allocation).

Avoiding concentration: risk-budgeted portfolios avoid overweighting high-volatility assets.

Mathematical formulation:

Find w such that: w_i * (Sigma @ w)_i / sigma_p = b_i * sigma_p

where b_i is the target risk budget (sum to 1).

Parameters:

cov (ndarray) – Covariance matrix (n x n).
target_risk (ndarray | None, default: None) – Target risk contribution vector (sums to 1). If None, uses equal risk contributions (1/n for each asset).

Returns:

weights (np.ndarray) – Optimal portfolio weights.
risk_contributions (np.ndarray) – Achieved risk contributions (should match target).
portfolio_vol (float) – Portfolio volatility.
converged (bool) – Whether the optimiser converged.

Return type:

Example

>>> import numpy as np
>>> cov = np.array([[0.04, 0.006], [0.006, 0.01]])
>>> result = risk_budgeting(cov)
>>> np.allclose(result["risk_contributions"], 0.5, atol=0.05)
True

See also

wraquant.risk.portfolio.risk_contribution: Compute risk: contributions for given weights.

wraquant.opt.portfolio: Full portfolio optimisation suite.

References

Maillard, Roncalli & Teiletche (2010), “The Properties of Equally Weighted Risk Contribution Portfolios”

diversification_ratio(weights, cov)[source]¶

Diversification ratio of a portfolio.

The diversification ratio is the ratio of the weighted average of individual asset volatilities to the portfolio volatility. It measures the diversification benefit captured by the portfolio.

When to use:

Use as a portfolio quality metric. Higher is better: - DR = 1.0: no diversification benefit (perfectly correlated). - DR = 1.5: good diversification. - DR > 2.0: excellent diversification. The Maximum Diversification Portfolio (Choueifaty & Coignard) maximises this ratio.

Mathematical formulation:

DR = (w’ * sigma) / sqrt(w’ * Sigma * w)

where sigma is the vector of individual asset volatilities and Sigma is the covariance matrix.

Parameters:

weights (ndarray) – Portfolio weight vector (n_assets,).
cov (ndarray) – Covariance matrix (n_assets x n_assets).

Return type:

Returns:

Diversification ratio as a float (>= 1.0).

Example

>>> import numpy as np
>>> cov = np.array([[0.04, 0.006], [0.006, 0.01]])
>>> weights = np.array([0.5, 0.5])
>>> dr = diversification_ratio(weights, cov)
>>> dr >= 1.0
True

See also

concentration_ratio: Herfindahl-based concentration measure. risk_budgeting: Find weights for target risk contributions.

References

Choueifaty & Coignard (2008), “Toward Maximum Diversification”

concentration_ratio(weights, cov)[source]¶

Herfindahl concentration ratio of risk contributions.

When to use:

Use concentration ratio to: - Detect hidden risk concentrations even when capital weights

look diversified. A portfolio with equal weights can still have concentrated risk if one asset is much more volatile.

Monitor risk concentration over time.
Compare portfolios: lower concentration ratio = more diversified risk.

Mathematical formulation:

CR = sum(rc_i^2) where rc_i is asset i’s fractional risk contribution (sum to 1.0).

CR = 1/n for equal risk contribution; CR = 1.0 for single-asset.

Parameters:

weights (ndarray) – Portfolio weight vector (n_assets,).
cov (ndarray) – Covariance matrix (n_assets x n_assets).

Return type:

Returns:

Herfindahl concentration ratio between 1/n and 1.0.

Example

>>> import numpy as np
>>> cov = np.array([[0.04, 0.0], [0.0, 0.04]])
>>> weights = np.array([0.5, 0.5])
>>> cr = concentration_ratio(weights, cov)
>>> abs(cr - 0.5) < 0.01  # equal vol + equal weight -> equal risk
True

See also

diversification_ratio: Alternative diversification metric.

tracking_error(returns, benchmark)[source]¶

Active risk metrics relative to a benchmark.

Tracking error (TE) is the standard deviation of the active return (portfolio return minus benchmark return). It measures how much the portfolio’s performance deviates from the benchmark.

When to use:: Use tracking error for: - Index tracking: target TE < 50bp for passive strategies. - Active management: typical TE of 2-8% for active equity funds. - Risk budgeting: allocate TE budget across portfolio managers.

Parameters:

returns (Series) – Portfolio return series.
benchmark (Series) – Benchmark return series (same frequency and index).

Returns:

tracking_error (float) – Annualized tracking error.
information_ratio (float) – Annualized active return / tracking error.
active_return (float) – Annualized mean active return.
max_active_drawdown (float) – Worst cumulative active return drawdown.
active_return_std (float) – Daily active return standard deviation (non-annualized).

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> portfolio = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> benchmark = pd.Series(np.random.normal(0.0004, 0.009, 252))
>>> result = tracking_error(portfolio, benchmark)
>>> result["tracking_error"] > 0
True

See also

active_share: Weight-based difference from benchmark. wraquant.risk.metrics.information_ratio: Simpler IR calculation.

active_share(weights, benchmark_weights)[source]¶

Active share: weight-based deviation from benchmark.

Active share measures how different a portfolio’s holdings are from its benchmark. It is computed as half the sum of absolute weight differences.

When to use:

Use active share to classify portfolio management style: - Active share < 20%: closet indexer (charging active fees for

passive exposure).

20-60%: moderate active.
60-80%: genuinely active.
> 80%: concentrated active or different investment universe.

Mathematical formulation:

Active Share = (1/2) * sum_i |w_i - w_bench_i|

Parameters:

weights (ndarray) – Portfolio weight vector.
benchmark_weights (ndarray) – Benchmark weight vector (same length).

Return type:

Returns:

Active share as a float between 0 and 1.

Example

>>> import numpy as np
>>> portfolio = np.array([0.4, 0.3, 0.2, 0.1])
>>> benchmark = np.array([0.25, 0.25, 0.25, 0.25])
>>> as_ = active_share(portfolio, benchmark)
>>> 0 <= as_ <= 1
True

See also

tracking_error: Return-based deviation from benchmark.

References

Cremers & Petajisto (2009), “How Active Is Your Fund Manager?”

Portfolio Risk¶

Portfolio-level risk analytics.

portfolio_volatility(weights, cov_matrix)[source]¶

Compute portfolio volatility from weights and covariance matrix.

Parameters:

weights (ndarray) – Asset weight vector (n,).
cov_matrix (ndarray) – Covariance matrix (n, n).

Return type:

Returns:

Portfolio volatility (standard deviation).

risk_contribution(weights, cov_matrix)[source]¶

Compute each asset’s risk contribution to portfolio volatility.

Parameters:

weights (ndarray) – Asset weight vector (n,).
cov_matrix (ndarray) – Covariance matrix (n, n).

Return type:

Returns:

Array of fractional risk contributions that sum to 1.

diversification_ratio(weights, cov_matrix)[source]¶

Compute the diversification ratio.

The diversification ratio is the ratio of the weighted average of individual volatilities to the portfolio volatility. A value of 1 means no diversification benefit.

Parameters:

weights (ndarray) – Asset weight vector (n,).
cov_matrix (ndarray) – Covariance matrix (n, n).

Return type:

Returns:

Diversification ratio (>= 1).

Copulas¶

Copula models for dependency structure in multi-asset returns.

Copulas separate the marginal distributions of individual assets from their dependence structure. This is critical in finance because linear correlation (Pearson) understates co-movement in the tails – precisely where risk matters most. Copulas capture non-linear, asymmetric, and tail-specific dependence that correlation matrices cannot.

This module provides five copula families:

Gaussian copula (fit_gaussian_copula) – the simplest elliptical copula. Captures linear dependence but has zero tail dependence: extreme co-movements are asymptotically independent. Useful as a baseline but dangerous for tail-risk applications.
Student-t copula (fit_t_copula) – symmetric tail dependence controlled by degrees of freedom. Lower df => heavier tails and stronger tail dependence. The standard choice for equity portfolios where joint crashes are common.
Clayton copula (fit_clayton_copula) – Archimedean copula with lower tail dependence and no upper tail dependence. Use when you expect assets to crash together but rally independently (typical for equities and credit).
Gumbel copula (fit_gumbel_copula) – Archimedean copula with upper tail dependence and no lower tail dependence. Use for assets that rally together but crash independently (rare but possible for certain commodity pairs).
Frank copula (fit_frank_copula) – Archimedean copula with symmetric dependence and no tail dependence. Useful as a benchmark or when tail dependence is genuinely absent.

Utilities:

copula_simulate: Monte Carlo simulation from any fitted copula.
tail_dependence: empirical estimation of lower and upper tail dependence coefficients.
rank_correlation: Kendall’s tau and Spearman’s rho matrices.

How to choose:

Start with tail_dependence to check if the data has asymmetric tail dependence.
If lower > upper: use Clayton.
If upper > lower: use Gumbel.
If roughly symmetric tails: use Student-t with estimated df.
If no significant tail dependence: use Gaussian or Frank.

References

Nelsen (2006), “An Introduction to Copulas”
McNeil, Frey & Embrechts (2005), Ch. 5, “Copulas and Dependence”
Embrechts, Lindskog & McNeil (2003), “Modelling Dependence with Copulas and Applications to Risk Management”

fit_gaussian_copula(returns)[source]¶

Fit a Gaussian copula to multivariate returns.

Mathematical formulation:

C(u_1, …, u_d) = Phi_R(Phi^{-1}(u_1), …, Phi^{-1}(u_d))

where Phi is the standard normal CDF, Phi^{-1} is its inverse (quantile function), and Phi_R is the multivariate normal CDF with correlation matrix R.

The density is:

c(u) = |R|^{-1/2} exp(-1/2 z^T (R^{-1} - I) z)

where z_i = Phi^{-1}(u_i).

Tail dependence:

lambda_L = lambda_U = 0 (for rho < 1)

This means the Gaussian copula predicts that extreme co-movements vanish in the tails — a dangerous assumption for financial risk.

Interpretation:

The correlation matrix R describes the “normal-like” dependence structure of the data.
Key limitation: the Gaussian copula implies that joint extreme events (crashes, rallies) are asymptotically independent. The 2008 crisis demonstrated that assets crash together far more than the Gaussian copula predicts.
If tail_dependence() returns non-negligible lower tail dependence (> 0.1), the Gaussian copula is inappropriate — use fit_t_copula or fit_clayton_copula instead.
The fitted correlation matrix R is NOT the same as the Pearson correlation of raw returns. It is the correlation of the rank-transformed (copula) data.

When to use:

As a baseline for dependence modelling.
When tail dependence is genuinely absent (rare in equities).
For quick simulation of correlated returns.
Compare its AIC/BIC against Student-t to test for tail dependence.

Parameters:

returns (ndarray) – Array of shape (n_obs, n_assets) with raw returns. Each column is one asset’s return series.

Returns:

"correlation" (ndarray) – estimated copula correlation matrix R of shape (d, d). Off-diagonal values measure the normal-copula dependence. Higher values = stronger co-movement.
"copula_type" (str) – "gaussian".

Return type:

Example

>>> import numpy as np
>>> from wraquant.risk.copulas import fit_gaussian_copula, fit_t_copula
>>> rng = np.random.default_rng(42)
>>> returns = rng.multivariate_normal([0, 0], [[1, 0.6], [0.6, 1]], 500)
>>> gauss = fit_gaussian_copula(returns)
>>> print(f"Copula correlation: {gauss['correlation'][0, 1]:.3f}")
>>> # Compare with t-copula to test for tail dependence
>>> t_cop = fit_t_copula(returns)
>>> print(f"t-copula df: {t_cop['df']:.1f} (lower = heavier tails)")

See also

fit_t_copula: Student-t copula with symmetric tail dependence.: Use when df < 30 to capture crash co-movement.

fit_clayton_copula: Lower-tail dependence only (crash contagion). tail_dependence: Empirically check if tail dependence exists. copula_simulate: Generate correlated samples from the fitted copula.

fit_t_copula(returns, df=5.0)[source]¶

Fit a Student-t copula to multivariate returns.

Mathematical formulation:

C(u_1, …, u_d) = t_{R,df}(t_df^{-1}(u_1), …, t_df^{-1}(u_d))

where t_df is the univariate Student-t CDF with df degrees of freedom, and t_{R,df} is the multivariate t CDF with correlation matrix R.

Tail dependence coefficient (bivariate case):

lambda_L = lambda_U = 2 * t_{df+1}(-sqrt((df+1)(1-rho)/(1+rho)))

Unlike the Gaussian copula, this is strictly positive for finite df, meaning extreme co-movements DO occur.

Interpretation:

df (degrees of freedom) controls tail heaviness:
- df = 3-5: Very heavy tails, strong tail dependence. Typical for equity portfolios during turbulent markets. lambda ~ 0.3-0.5 for rho=0.5.
- df = 10-20: Moderate tails. Appropriate for investment- grade fixed income or diversified portfolios.
- df > 30: Nearly Gaussian (tail dependence < 0.05).
- df -> inf: Converges to Gaussian copula exactly.
Compare df across time periods: if df drops from 15 to 5, tail dependence has increased — crisis contagion is building.
Compare fitted df to Gaussian AIC: if t-copula AIC is much lower, tail dependence is statistically significant.

When to use:

Default choice for equity and credit risk modelling.
VaR/CVaR estimation for multi-asset portfolios.
When you expect both crashes AND rallies to be correlated (symmetric dependence).
If you need asymmetric tail dependence (crashes correlated but rallies independent), use fit_clayton_copula instead.

Parameters:

returns (ndarray) – Array of shape (n_obs, n_assets) with raw returns.
df (float, default: 5.0) – Degrees of freedom. Lower = heavier tails. Typical values: 3-5 for equities, 10-20 for bonds, 5-10 as a safe default. If unsure, start with 5 and compare AIC against df=10,15,20.

Returns:

"correlation" (ndarray) – Copula correlation matrix R.
"df" (float) – Degrees of freedom used.
"copula_type" (str) – "student_t".

Return type:

Example

>>> import numpy as np
>>> from wraquant.risk.copulas import fit_t_copula, tail_dependence
>>> rng = np.random.default_rng(42)
>>> # Simulate fat-tailed correlated returns
>>> returns = rng.multivariate_normal([0, 0], [[1, 0.5], [0.5, 1]], 1000)
>>> result = fit_t_copula(returns, df=5)
>>> print(f"Copula corr: {result['correlation'][0, 1]:.3f}")
>>> print(f"df={result['df']} → heavier tails than Gaussian")
>>> # Check tail dependence
>>> td = tail_dependence(returns)
>>> print(f"Lower tail dep: {td['lower']:.3f}")

Notes

Reference: Demarta, S. & McNeil, A.J. (2005). “The t Copula and Related Copulas.” International Statistical Review, 73, 111-129.

See also

fit_gaussian_copula: No tail dependence (Gaussian = t with df=inf). fit_clayton_copula: Asymmetric — lower tail dependence only. tail_dependence: Empirically estimate lambda_L and lambda_U.

fit_clayton_copula(u, v)[source]¶

Fit a Clayton copula (lower tail dependence) to bivariate data.

Interpretation:

theta > 0 controls dependence strength. Higher theta = stronger lower tail dependence. theta -> 0 = independence.
lower_tail_dependence = 2^{-1/theta} is the probability that both assets are in their lower tail simultaneously, given that one is. Values above 0.3 indicate meaningful crash co-movement.
The Clayton copula is the canonical choice for modelling “contagion” and “flight to quality” dynamics.

When to use:

Two equity assets or sectors where you expect crash contagion but independent rallies.
Credit portfolios where defaults cluster.
When tail_dependence shows lower > upper.

Parameters:

u (ndarray) – First marginal uniform observations. Pass raw data and the empirical CDF is applied automatically.
v (ndarray) – Second marginal uniform observations.

Returns:

"theta" – Clayton parameter (> 0).
"copula_type" – "clayton".
"lower_tail_dependence" – 2^{-1/theta}. > 0.3 is meaningful.

Return type:

See also

fit_gumbel_copula: Upper tail dependence (opposite asymmetry). fit_t_copula: Symmetric tail dependence.

fit_gumbel_copula(u, v)[source]¶

Fit a Gumbel copula (upper tail dependence) to bivariate data.

Interpretation:

theta >= 1 controls dependence strength. theta = 1 is independence; larger theta = stronger upper tail dependence.
upper_tail_dependence = 2 - 2^{1/theta} measures the probability of joint extreme upper co-movement.

When to use:

Commodity pairs with supply-shock co-movement.
When tail_dependence shows upper > lower.
When modelling “melt-up” contagion.

Parameters:

u (ndarray) – First marginal uniform observations.
v (ndarray) – Second marginal uniform observations.

Returns:

"theta" – Gumbel parameter (>= 1).
"copula_type" – "gumbel".
"upper_tail_dependence" – 2 - 2^{1/theta}.

Return type:

See also

fit_clayton_copula: Lower tail dependence (more common in equities). fit_t_copula: Symmetric tail dependence.

fit_frank_copula(u, v)[source]¶

Fit a Frank copula (symmetric dependence, no tail dependence).

Interpretation:

theta > 0: positive dependence; theta < 0: negative dependence; |theta| large = strong dependence.
The Frank copula allows negative dependence (unlike Clayton and Gumbel), making it useful for hedging pairs.
No tail dependence: extreme co-movements are modelled as asymptotically independent.

When to use:

As a benchmark against Clayton/Gumbel.
When tail_dependence shows negligible tail dependence.
For pairs with negative dependence (hedging relationships).

Parameters:

u (ndarray) – First marginal uniform observations.
v (ndarray) – Second marginal uniform observations.

Returns:

"theta" – Frank parameter. Positive = positive dependence, negative = negative dependence.
"copula_type" – "frank".

Return type:

See also

fit_gaussian_copula: Multivariate with no tail dependence. fit_clayton_copula: Lower tail dependence. fit_gumbel_copula: Upper tail dependence.

copula_simulate(copula_params, n_sims=10000, copula_type=None, seed=None)[source]¶

Simulate from a fitted copula.

Interpretation:

Output columns are uniform on (0, 1) – they represent the dependence structure only, not the marginal distributions.
Correlated low values in all columns = a joint crash scenario.
To convert to returns: returns[:, j] = norm.ppf(U[:, j], loc=mu_j, scale=sigma_j) or use empirical quantile functions for non-parametric marginals.
Use for Monte Carlo VaR/CVaR, portfolio simulation, or stress testing.

Parameters:

copula_params (dict[str, Any]) – Output from one of the fit_*_copula functions.
n_sims (int, default: 10000) – Number of samples to draw.
copula_type (str | None, default: None) – Override copula type (defaults to copula_params["copula_type"]).
seed (int | None, default: None) – Random seed for reproducibility.

Return type:

Returns:

Array of shape (n_sims, d) with uniform marginals in (0, 1).

Raises:

ValueError – If the copula type is not recognized.

Example

>>> from scipy.stats import norm
>>> cop = fit_gaussian_copula(returns)
>>> U = copula_simulate(cop, n_sims=10000, seed=42)
>>> # Transform to normal marginals:
>>> sim_returns = norm.ppf(U, loc=mu, scale=sigma)

tail_dependence(u, v, method='empirical', threshold=0.05)[source]¶

Estimate lower and upper tail dependence coefficients.

Tail dependence measures the probability that one variable is in its extreme tail given that the other is. This is the key quantity that copula selection hinges on.

Interpretation:

lower ~ P(V <= q | U <= q): the probability of a joint crash. Values above 0.2-0.3 are economically significant and indicate that the Gaussian copula is inappropriate.
upper ~ P(V > 1-q | U > 1-q): the probability of a joint rally.
If lower >> upper: use Clayton copula (crash contagion).
If upper >> lower: use Gumbel copula (rally contagion).
If both are similar: use Student-t copula (symmetric tails).
If both are near zero: Gaussian or Frank copula is adequate.

Caveat:

Empirical tail dependence estimates are noisy with small samples. Use threshold = 0.10 for more observations per tail (less extreme, more precise) or 0.05 for fewer observations (more extreme, noisier).

Parameters:

u (ndarray) – First marginal uniform observations.
v (ndarray) – Second marginal uniform observations.
method (str, default: 'empirical') – "empirical" (default) for non-parametric estimation.
threshold (float, default: 0.05) – Quantile threshold for tail estimation. 0.05 = 5th/95th percentile (default). 0.10 = more data in the tail estimate but less extreme.

Return type:

Returns:

Dict with "lower" and "upper" tail dependence estimates.

Example

>>> import numpy as np
>>> from scipy.stats import norm
>>> rng = np.random.default_rng(0)
>>> # Simulate from a t-copula (symmetric tail dependence)
>>> z = rng.multivariate_normal([0,0], [[1,0.5],[0.5,1]], 5000)
>>> u = norm.cdf(z[:, 0])
>>> v = norm.cdf(z[:, 1])
>>> td = tail_dependence(u, v)
>>> print(f"Lower: {td['lower']:.2f}, Upper: {td['upper']:.2f}")

rank_correlation(returns, method='both')[source]¶

Compute Kendall’s tau and/or Spearman’s rho rank correlation.

Interpretation:

Kendall’s tau: Probability of concordance minus probability of discordance. tau = 0.5 means ~75% of pairs move in the same direction.
Spearman’s rho: Pearson correlation of the ranks. Generally |rho| >= |tau| for the same data.
Both are robust to outliers (unlike Pearson).
Copula relationships: for the Gaussian copula, rho_pearson = 2*sin(pi*tau/6). For Clayton, theta = 2*tau/(1-tau).

When to use:

Always prefer rank correlation over Pearson for copula modelling.
Use Kendall’s tau for small samples (more robust).
Use Spearman’s rho for comparison with Pearson.

Parameters:

returns (ndarray) – Array of shape (n_obs, n_assets).
method (str, default: 'both') – "kendall", "spearman", or "both" (default).

Returns:

"kendall_tau" – Kendall’s tau correlation matrix.
"spearman_rho" – Spearman’s rho correlation matrix.

Return type:

Raises:

ValueError – If method is not recognized.

Tail Risk & EVT¶

Tail risk analytics for non-normal return distributions.

Standard risk measures (VaR, volatility) assume or approximate normality. In practice, financial returns are fat-tailed (excess kurtosis) and left-skewed. This module provides tail-aware risk measures that account for higher moments and drawdown-based risk.

Functions:

Cornish-Fisher VaR – adjusts the normal VaR quantile for skewness and kurtosis using the Cornish-Fisher expansion.
ES decomposition – per-asset contribution to Expected Shortfall.
Conditional Drawdown at Risk (CDaR) – the average of worst-alpha% drawdowns (analogous to CVaR but for drawdowns).
Tail ratio analysis – 95th/5th percentile ratio with diagnostics.
Drawdown at Risk (DaR) – worst alpha-quantile drawdown.

References

Cornish & Fisher (1937), “Moments and Cumulants in the Specification of Distributions”
Chekhlov, Uryasev & Zabarankin (2005), “Drawdown Measure in Portfolio Optimization”

cornish_fisher_var(returns, alpha=0.05)[source]¶

Cornish-Fisher expansion VaR (skewness and kurtosis adjusted).

When to use:

Use Cornish-Fisher VaR when: - Returns are detectably non-normal (skewness != 0 or kurtosis != 3). - You want a quick analytical adjustment without fitting a full

distribution (e.g., Student-t or EVT).

The sample is too short for reliable historical VaR but long enough to estimate skewness/kurtosis (>100 observations).

Mathematical formulation:

z_cf = z + (z^2 - 1) * S/6 + (z^3 - 3z) * K/24 - (2z^3 - 5z) * S^2/36

CF-VaR = -(mu + sigma * z_cf)

where z = Phi^{-1}(alpha), S = skewness, K = excess kurtosis.

How to interpret:

Parameters:

returns (Series) – Simple return series.
alpha (float, default: 0.05) – Significance level (0.05 = 95% VaR).

Returns:

cf_var (float) – Cornish-Fisher adjusted VaR (positive number = loss).
normal_var (float) – Standard parametric (Gaussian) VaR.
z_cf (float) – Adjusted quantile.
z_normal (float) – Standard normal quantile.
skewness (float) – Sample skewness.
excess_kurtosis (float) – Sample excess kurtosis.
adjustment_factor (float) – cf_var / normal_var.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0, 0.01, 1000))
>>> result = cornish_fisher_var(returns, alpha=0.05)
>>> result["cf_var"] > 0
True

See also

wraquant.risk.var.value_at_risk: Historical and parametric VaR. tail_ratio_analysis: Tail shape diagnostics.

References

Cornish & Fisher (1937), “Moments and Cumulants in the Specification of Distributions”
Maillard (2012), “A User’s Guide to the Cornish Fisher Expansion”

expected_shortfall_decomposition(weights, returns, alpha=0.05)[source]¶

Decompose Expected Shortfall (CVaR) into per-asset contributions.

When to use:

Use ES decomposition for: - Identifying which assets drive tail losses. - Setting per-asset ES limits. - Comparing tail-risk concentration to normal-market risk

concentration.

Mathematical formulation:

ES_i = w_i * E[r_i | r_p <= VaR_alpha(r_p)]

where r_p = w’ @ r is the portfolio return.

Parameters:

weights (ndarray) – Portfolio weight vector (n_assets,).
returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
alpha (float, default: 0.05) – Significance level (0.05 = worst 5% of days).

Return type:

Returns:

pd.Series of per-asset ES contributions. Sum equals portfolio ES (as a positive number).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "A": np.random.normal(0.0005, 0.01, 500),
...     "B": np.random.normal(0.0003, 0.015, 500),
... })
>>> weights = np.array([0.6, 0.4])
>>> es = expected_shortfall_decomposition(weights, returns, alpha=0.05)
>>> es.sum() > 0  # positive = loss
True

See also

component_var: Euler decomposition of VaR. cornish_fisher_var: Skewness-adjusted VaR.

conditional_drawdown_at_risk(returns, alpha=0.05)[source]¶

Conditional Drawdown at Risk (CDaR).

When to use:: Use CDaR when drawdown is a primary risk constraint (e.g., hedge funds with max drawdown mandates). CDaR penalises sustained drawdowns, not just point-in-time losses. A portfolio optimised to minimise CDaR will have better drawdown recovery properties than one optimised for VaR.

Parameters:

returns (Series) – Simple return series.
alpha (float, default: 0.05) – Fraction of worst drawdowns to average (0.05 = worst 5%).

Return type:

Returns:

CDaR as a positive float (e.g., 0.15 = average worst-5% drawdown is 15%).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0.0005, 0.01, 500))
>>> cdar = conditional_drawdown_at_risk(returns, alpha=0.05)
>>> cdar >= 0
True

See also

drawdown_at_risk: Quantile-based drawdown measure (DaR). wraquant.risk.metrics.max_drawdown: Single worst drawdown.

References

Chekhlov, Uryasev & Zabarankin (2005), “Drawdown Measure in Portfolio Optimization”

tail_ratio_analysis(returns)[source]¶

Tail ratio analysis with interpretation.

When to use:

Use tail ratio analysis to: - Assess payoff asymmetry: trend-following should have tail ratio > 1

(large gains, small frequent losses).

Detect negative skew: mean-reversion and short vol strategies typically have tail ratio < 1.
Compare strategies beyond Sharpe ratio.

Parameters:

returns (Series) – Simple return series.

Returns:

tail_ratio (float) – 95th percentile / abs(5th percentile).
right_tail (float) – 95th percentile return.
left_tail (float) – 5th percentile return.
tail_ratio_99 (float) – 99th/1st percentile ratio.
skewness (float) – Sample skewness.
excess_kurtosis (float) – Sample excess kurtosis.
interpretation (str) – Human-readable assessment.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0, 0.01, 1000))
>>> result = tail_ratio_analysis(returns)
>>> result["tail_ratio"] > 0
True

See also

cornish_fisher_var: Skewness-adjusted VaR.

drawdown_at_risk(returns, alpha=0.05)[source]¶

Drawdown at Risk (DaR): worst alpha-quantile drawdown.

DaR is to drawdowns what VaR is to returns. It is the alpha-percentile of the drawdown distribution – i.e., the drawdown that is exceeded only alpha% of the time.

When to use:: Use DaR when setting drawdown limits: - “With 95% confidence, the drawdown will not exceed DaR.” - Useful for fund prospectuses and investor communications. - More intuitive than VaR for many stakeholders because

drawdowns are easier to understand than daily P&L.

Parameters:

returns (Series) – Simple return series.
alpha (float, default: 0.05) – Significance level (0.05 = 95th percentile drawdown).

Return type:

Returns:

DaR as a positive float (e.g., 0.12 = 12% drawdown).

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0.0005, 0.01, 500))
>>> dar = drawdown_at_risk(returns, alpha=0.05)
>>> dar >= 0
True

See also

conditional_drawdown_at_risk: Average of worst drawdowns (CDaR). wraquant.risk.metrics.max_drawdown: Single worst drawdown.

DCC Multivariate¶

Dynamic Conditional Correlation (DCC-GARCH) models.

Provides univariate GARCH(1,1) fitting, DCC parameter estimation via MLE, rolling DCC-based correlations, correlation forecasting, and time-varying covariance matrix computation. All implementations use pure numpy/scipy.

dcc_garch(returns, p=1, q=1)[source]¶

Fit a DCC-GARCH(p, q) model.

Currently supports p=1, q=1 (DCC(1,1)).

Procedure:

Fit univariate GARCH(1,1) to each return series.
Compute standardized residuals.
Estimate DCC parameters (a, b) via MLE.

Parameters:

returns (ndarray) – Array of shape (T, k) with k asset return series.
p (int, default: 1) – DCC lag order for innovation (currently must be 1).
q (int, default: 1) – DCC lag order for conditional (currently must be 1).

Returns:

"a" – DCC innovation parameter.
"b" – DCC persistence parameter.
"garch_params" – list of per-asset GARCH(1,1) parameter dicts.
"qbar" – unconditional correlation matrix of standardized residuals.
"conditional_vols" – array of per-asset conditional volatilities (T, k).
"std_residuals" – standardized residuals (T, k).

Return type:

rolling_correlation_dcc(returns, window=None)[source]¶

Compute DCC-based time-varying correlation matrices.

Fits DCC-GARCH to the full sample, then extracts the time-varying correlation matrix at each time step.

Parameters:

returns (ndarray) – Array of shape (T, k).
window (int | None, default: None) – Ignored (present for API compatibility with simple rolling-window methods). DCC uses the full sample.

Returns:

"correlations" – array of shape (T, k, k) with the time-varying correlation matrices.
"dcc_model" – the fitted DCC model dict.

Return type:

forecast_correlation(dcc_model, horizon=1)[source]¶

Forecast future correlation matrices from a fitted DCC model.

Parameters:

dcc_model (dict[str, Any]) – Output of dcc_garch().
horizon (int, default: 1) – Number of steps ahead to forecast.

Returns:

"forecasted_correlations" – array of shape (horizon, k, k).
"forecasted_covariances" – array of shape (horizon, k, k).

Return type:

conditional_covariance(returns, dcc_params=None)[source]¶

Compute time-varying covariance matrices from DCC-GARCH.

If dcc_params is not supplied, the model is fitted to returns automatically.

Parameters:

returns (ndarray) – Array of shape (T, k).
dcc_params (dict[str, Any] | None, default: None) – Pre-fitted DCC model dict (optional).

Returns:

"covariances" – array of shape (T, k, k) with conditional covariance matrices.
"correlations" – array of shape (T, k, k) with conditional correlation matrices.
"volatilities" – array of shape (T, k) with conditional volatilities.

Return type:

Monte Carlo¶

Advanced Monte Carlo methods for risk measurement.

Variance reduction techniques, bootstrap methods, and filtered historical simulation for improved VaR/ES estimation.

antithetic_variates(mu, sigma, n_sims, n_assets=1, seed=None)[source]¶

Generate antithetic variate samples for variance reduction.

Produces n_sims paired samples (original + antithetic) from a normal distribution. The antithetic counterpart mirrors each draw about the mean, reducing variance for monotone functions.

Parameters:

mu (float | ndarray) – Mean(s) of the distribution. Scalar or array of length n_assets.
sigma (float | ndarray) – Standard deviation(s). Scalar or array of length n_assets.
n_sims (int) – Number of pairs to generate (total output = 2 * n_sims rows).
n_assets (int, default: 1) – Number of assets / dimensions.
seed (int | None, default: None) – Random seed for reproducibility.

Return type:

Returns:

Array of shape (2 * n_sims, n_assets) containing the original and antithetic draws interleaved.

block_bootstrap(returns, block_size, n_sims=1000, seed=None)[source]¶

Block bootstrap for autocorrelated time series.

Resamples contiguous blocks of the input series to preserve serial dependence (Kunsch 1989, Liu & Singh 1992).

Parameters:

returns (ndarray) – 1-D array of return observations.
block_size (int) – Length of each block.
n_sims (int, default: 1000) – Number of bootstrap replications.
seed (int | None, default: None) – Random seed for reproducibility.

Return type:

Returns:

2-D array of shape (n_sims, len(returns)) where each row is one bootstrap replicate.

filtered_historical_simulation(returns, vol_model='ewma', decay=0.94, n_sims=1000, seed=None)[source]¶

Filtered historical simulation (FHS).

Combines a volatility model (EWMA or simple GARCH(1,1)) with historical bootstrap of standardised residuals to produce volatility-adjusted scenario returns.

Parameters:

returns (ndarray) – 1-D array of historical returns.
vol_model (str, default: 'ewma') – Volatility model — "ewma" (default) or "garch".
decay (float, default: 0.94) – EWMA decay factor (lambda) or, for "garch", the persistence parameter beta.
n_sims (int, default: 1000) – Number of simulated next-period returns.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

simulated_returns: 1-D array of n_sims simulated next-period returns.
current_vol: Estimated current volatility used for rescaling.
standardised_residuals: Standardised residuals from the volatility model.

Return type:

importance_sampling_var(returns, n_sims=10000, target_quantile=0.01, shift=None, seed=None)[source]¶

Estimate VaR via importance sampling.

Shifts the sampling distribution toward the tail to obtain more accurate estimates of extreme quantiles with fewer simulations.

Parameters:

returns (ndarray) – 1-D array of historical returns.
n_sims (int, default: 10000) – Number of Monte Carlo draws.
target_quantile (float, default: 0.01) – Quantile level for VaR (e.g., 0.01 for 1%).
shift (float | None, default: None) – Mean shift for the importance sampling distribution. If None, automatically set to the empirical quantile.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

var: Estimated VaR (positive number = loss).
effective_sample_size: Effective sample size after reweighting, as a fraction of n_sims.

Return type:

stationary_bootstrap(returns, avg_block_size=10.0, n_sims=1000, seed=None)[source]¶

Stationary bootstrap with random block sizes (Politis & Romano 1994).

Block lengths follow a geometric distribution with mean avg_block_size, producing a strictly stationary resampled series.

Parameters:

returns (ndarray) – 1-D array of return observations.
avg_block_size (float, default: 10.0) – Expected block size (must be > 1).
n_sims (int, default: 1000) – Number of bootstrap replications.
seed (int | None, default: None) – Random seed for reproducibility.

Return type:

Returns:

2-D array of shape (n_sims, len(returns)) where each row is one bootstrap replicate.

stratified_sampling(returns, n_strata=10, n_sims=10000, seed=None)[source]¶

Stratified sampling for VaR estimation.

Divides the probability space into equal strata and draws uniformly within each stratum, ensuring better coverage of the tails.

Parameters:

returns (ndarray) – 1-D array of historical returns (used to fit a normal distribution).
n_strata (int, default: 10) – Number of strata.
n_sims (int, default: 10000) – Total number of draws (distributed evenly across strata).
seed (int | None, default: None) – Random seed for reproducibility.

Return type:

Returns:

1-D array of n_sims stratified samples drawn from the fitted normal distribution.

Stress Testing¶

Comprehensive stress testing for portfolio risk analysis.

Stress testing answers “what if?” questions that historical VaR and CVaR cannot address. While VaR extrapolates from the empirical distribution, stress tests evaluate specific adverse scenarios – including scenarios that have never occurred in the sample.

This module provides seven complementary stress testing approaches:

Scenario-based (stress_test_returns) – apply user-defined additive shocks (e.g., “what if every day’s return is 10% worse?”).
Historical replay (historical_stress_test) – measure portfolio performance during known crises (GFC, COVID, dot-com).
Volatility scaling (vol_stress_test) – scale return dispersion by multipliers (1.5x, 2x, 3x) while preserving the mean.
Spot stress (spot_stress_test) – shift price levels by percentage amounts (-30% to +10%).
Sensitivity ladder (sensitivity_ladder) – P&L sensitivity to a single factor across a range of shock levels.
Reverse stress test (reverse_stress_test) – find the scenarios that produce a target loss (regulatory requirement).
Joint stress test (joint_stress_test) – simultaneous volatility, spot, and correlation shocks.
Marginal contribution (marginal_stress_contribution) – identify the asset contributing most to stress loss.

How to interpret stress test results:: Stress tests produce point estimates, not probability distributions. The output tells you “if X happens, the P&L impact is Y.” They are valuable for: - Setting position limits and stop-losses. - Capital adequacy assessment (CCAR, DFAST). - Identifying concentrated risk exposures. - Communicating tail risk to stakeholders.

References

Berkowitz (2000), “A Coherent Framework for Stress-Testing”
McNeil, Frey & Embrechts (2005), “Quantitative Risk Management”

stress_test_returns(returns, scenarios)[source]¶

Apply user-defined additive shock scenarios to a return series.

When to use:: Use when you want to evaluate the impact of a uniform adverse shift in returns. For example, “what if every daily return is 10bp worse due to a funding cost shock?” For more realistic scenario analysis, use historical_stress_test (replays real crises) or joint_stress_test (simultaneous multi-factor shocks).
How to interpret:: Compare stressed_var_95 and stressed_cvar_95 across scenarios to identify which shock level pushes your portfolio into unacceptable loss territory. If a moderate shock (-5%) already produces a severe stressed CVaR, the portfolio has insufficient risk budget.

Parameters:

returns (Series | DataFrame) – Historical return series (Series) or multi-asset returns (DataFrame). For DataFrames, the cross-asset mean is used.
scenarios (dict[str, float]) – Mapping of scenario name to additive return shock (e.g. {"crash": -0.10, "boom": 0.05}). A shock of -0.10 subtracts 10% from every observation.

Returns:

"scenario_results" – dict mapping scenario name to a dict with "stressed_mean", "stressed_var_95" (5th percentile of stressed returns), "stressed_cvar_95" (mean of returns below the 5th percentile).
"base_mean" – mean of the original (unstressed) returns.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0.0005, 0.01, 252))
>>> result = stress_test_returns(returns, {"mild": -0.005, "severe": -0.02})
>>> result["scenario_results"]["severe"]["stressed_mean"] < result["base_mean"]
True

See also

historical_stress_test: Replay known crisis periods. vol_stress_test: Scale volatility by multipliers. joint_stress_test: Simultaneous multi-factor shocks.

historical_stress_test(returns, crisis_periods=None)[source]¶

Test portfolio returns against known historical crisis periods.

When to use:

Use historical stress testing for: - Board and regulator presentations (“how would we have

performed in the GFC?”).

Identifying whether the portfolio’s risk profile has improved or deteriorated relative to past crises.
Calibrating position limits against known worst cases.

How to interpret:

Built-in crisis periods (used when crisis_periods is None):

GFC 2008: 2008-09-01 to 2009-03-31
COVID 2020: 2020-02-19 to 2020-03-23
Dot-Com 2000: 2000-03-10 to 2002-10-09
Euro Debt 2011: 2011-07-01 to 2011-11-30
Taper Tantrum 2013: 2013-05-22 to 2013-09-05
Volmageddon 2018: 2018-02-02 to 2018-02-08
Flash Crash 2010: 2010-05-06

Parameters:

returns (Series | DataFrame) – Return series with a DatetimeIndex.
crisis_periods (dict[str, tuple[str, str]] | None, default: None) – Mapping of crisis name to (start, end) date strings. Periods not covered by the data are skipped.

Returns:

"crisis_results" – dict mapping crisis name to a dict with "cumulative_return" (compounded return over the crisis), "max_drawdown" (worst peak-to-trough within the crisis), "mean_daily_return", "n_days".
"periods_found" – list of crisis names that overlap with the data.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> idx = pd.bdate_range("2008-01-01", "2009-12-31")
>>> returns = pd.Series(np.random.normal(-0.001, 0.02, len(idx)), index=idx)
>>> result = historical_stress_test(returns)
>>> "gfc_2008" in result["periods_found"]
True

See also

stress_test_returns: User-defined additive shocks. reverse_stress_test: Find scenarios that produce a target loss.

vol_stress_test(returns, vol_shocks=None)[source]¶

Stress test by scaling return volatility with multipliers.

When to use:: Use volatility stress tests to: - Assess margin adequacy under elevated vol regimes. - Calibrate dynamic position sizing rules. - Compare the portfolio’s sensitivity to vol scaling

(a diversified portfolio should be less sensitive than a concentrated one).
How to interpret:: The stressed_vol should scale linearly with the multiplier (by construction). The key outputs are stressed_var_95 and stressed_cvar_95: if doubling vol (multiplier 2.0) causes CVaR to more than double, the portfolio has convex (nonlinear) tail exposure.

Parameters:

returns (Series | DataFrame) – Historical return series.
vol_shocks (list[float] | None, default: None) – List of volatility multipliers (e.g. [1.5, 2.0, 3.0]). Defaults to [1.5, 2.0, 2.5, 3.0].

Returns:

"vol_results" – dict mapping multiplier (as string) to "stressed_vol", "stressed_var_95", "stressed_cvar_95", "stressed_mean".
"base_vol" – volatility of the original returns.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0.0005, 0.01, 500))
>>> result = vol_stress_test(returns, vol_shocks=[1.5, 2.0])
>>> result["vol_results"]["2.0"]["stressed_vol"] > result["base_vol"]
True

See also

stress_test_returns: Additive shock scenarios. joint_stress_test: Combined vol, spot, and correlation shocks.

spot_stress_test(prices, spot_shocks=None)[source]¶

Shift spot (price) levels by specified percentage amounts.

Each shock is applied as a multiplicative factor to the final price (e.g., -0.10 means a 10% drop from the last price). This is useful for mark-to-market stress testing of current positions.

When to use:: Use spot stress tests for: - Options portfolio Greeks analysis (delta P&L under spot move). - Margin calculation under adverse spot scenarios. - Reporting to counterparties (“what is my exposure if the

underlying drops 20%?”).
How to interpret:: The shocked_price shows the resulting price level after the shock. For a DataFrame (multi-asset), the same percentage shock is applied to each asset’s last price. price_change is the absolute dollar change.

Parameters:

prices (Series | DataFrame) – Price series or DataFrame of asset prices.
spot_shocks (list[float] | None, default: None) – List of percentage shocks (e.g. [-0.30, -0.20, -0.10, 0.10, 0.20]). Defaults to [-0.30, -0.20, -0.10, -0.05, 0.05, 0.10].

Returns:

"spot_results" – dict mapping shock (as string) to "shocked_price", "price_change", "pct_change".
"base_price" – the last observed price.

Return type:

Example

>>> import pandas as pd
>>> prices = pd.Series([100.0, 102.0, 101.0, 103.0])
>>> result = spot_stress_test(prices, spot_shocks=[-0.10, 0.10])
>>> result["spot_results"]["-0.1"]["shocked_price"]
92.7

See also

vol_stress_test: Scale volatility by multipliers. sensitivity_ladder: P&L sensitivity to a single factor.

sensitivity_ladder(portfolio_returns, factor_returns, shock_range=None)[source]¶

Compute portfolio P&L across a range of factor shocks.

When to use:

Use sensitivity ladders to: - Understand how exposed the portfolio is to a single risk

factor (e.g., S&P 500, 10Y yield, oil price).

Construct hedging ratios (the beta tells you how much factor exposure to neutralise).
Present risk to traders and PMs in an intuitive format.

Mathematical formulation:

Step 1: Fit r_p = alpha + beta * r_f + epsilon via OLS. Step 2: For each shock s, estimate P&L = alpha + beta * s.

How to interpret:

Parameters:

portfolio_returns (Series) – Portfolio return series.
factor_returns (Series) – Single-factor return series (same index).
shock_range (ndarray | list[float] | None, default: None) – Array of factor shock values (e.g. np.linspace(-0.10, 0.10, 21)). Defaults to np.linspace(-0.10, 0.10, 21).

Returns:

"ladder" – dict mapping shock level (float) to estimated portfolio P&L.
"beta" – regression beta (sensitivity).
"alpha" – regression intercept (return when factor = 0).
"r_squared" – R-squared of the regression.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> factor = pd.Series(np.random.normal(0, 0.01, 252))
>>> portfolio = 0.8 * factor + np.random.normal(0, 0.005, 252)
>>> result = sensitivity_ladder(portfolio, factor)
>>> abs(result["beta"] - 0.8) < 0.2  # beta close to true value
True

See also

spot_stress_test: Direct price-level shocks. joint_stress_test: Multi-factor simultaneous shocks.

reverse_stress_test(returns, target_loss, n_sims=10000, seed=None)[source]¶

Find scenarios that produce at least the specified target loss.

When to use:

Use reverse stress tests when you need to: - Identify the conditions under which the portfolio breaches

a risk limit (e.g., -20% annual loss).

Satisfy regulatory requirements for reverse stress testing.
Understand how “unlikely” a catastrophic loss really is.

How to interpret:

Caveats:

Parameters:

returns (Series | DataFrame) – Historical return series.
target_loss (float) – Target cumulative loss as a negative number (e.g. -0.20 for a 20% loss).
n_sims (int, default: 10000) – Number of Monte Carlo paths to simulate.
seed (int | None, default: None) – Random seed for reproducibility.

Returns:

"scenarios_found" – number of simulated paths that hit the target.
"probability" – estimated probability of hitting the target.
"avg_loss" – mean loss across qualifying scenarios.
"worst_loss" – worst loss observed in qualifying scenarios.
"threshold_percentile" – percentile at which the target loss sits in the simulated distribution.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.Series(np.random.normal(0.0003, 0.01, 252))
>>> result = reverse_stress_test(returns, target_loss=-0.30, n_sims=5000, seed=42)
>>> result["probability"] >= 0
True

See also

historical_stress_test: Replay known crisis periods. stress_test_returns: User-defined additive shocks.

joint_stress_test(returns, vol_shock=2.0, spot_shock=-0.1, correlation_shock=0.0)[source]¶

Apply combined volatility, spot, and correlation shocks.

When to use:

Use joint stress tests for: - Portfolio optimisation stress testing: feed the stressed

covariance matrix into a mean-variance optimiser.

Capital adequacy under combined adverse conditions.
Comparing diversification benefits under normal vs. stressed conditions (correlation shock toward 1.0 eliminates diversification).

Procedure:

Scale demeaned returns by vol_shock (volatility multiplier).
Shift the mean by spot_shock (additive level shift).
Blend the correlation matrix toward uniform correlation: stressed_corr = (1 - c) * corr + c * ones_matrix, where c = correlation_shock.

How to interpret:

Parameters:

returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
vol_shock (float, default: 2.0) – Volatility multiplier (e.g. 2.0 = double vol).
spot_shock (float, default: -0.1) – Additive shift to mean returns (e.g. -0.10 = subtract 10% from each asset’s mean return).
correlation_shock (float, default: 0.0) – Blend factor toward perfect correlation. 0 = unchanged, 0.5 = halfway to perfect correlation, 1 = all pairwise correlations set to 1.0.

Returns:

"stressed_mean" – stressed mean return per asset (dict).
"stressed_vol" – stressed volatility per asset (dict).
"stressed_corr" – stressed correlation matrix (ndarray).
"stressed_cov" – stressed covariance matrix (ndarray).
"base_mean" – original mean returns (dict).
"base_vol" – original volatilities (dict).

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "SPY": np.random.normal(0.0005, 0.01, 252),
...     "TLT": np.random.normal(0.0002, 0.005, 252),
... })
>>> result = joint_stress_test(returns, vol_shock=2.0, correlation_shock=0.5)
>>> result["stressed_vol"]["SPY"] > result["base_vol"]["SPY"]
True

See also

vol_stress_test: Volatility scaling only. marginal_stress_contribution: Identify worst-contributing asset.

marginal_stress_contribution(portfolio_weights, returns, scenario)[source]¶

Identify which asset contributes most to portfolio stress loss.

When to use:: Use marginal stress contribution after running a stress test to answer: “which position is killing the portfolio under this scenario?” This guides targeted hedging decisions (e.g., buy puts on the worst-contributing asset).
How to interpret:: asset_contributions shows each asset’s dollar P&L under the scenario (weight_i * scenario_return_i). pct_contributions normalises these to sum to 1.0, showing the fraction of total loss attributable to each asset. worst_asset is the asset with the most negative contribution.

Parameters:

portfolio_weights (ndarray) – Weight vector aligned with returns.columns. Must have the same length as the number of columns.
returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
scenario (dict[str, float]) – Mapping of asset name to shocked return value. Assets not in the scenario use their historical mean return.

Returns:

"total_stress_loss" – portfolio return under the scenario.
"asset_contributions" – dict mapping asset name to its P&L contribution (weight * return).
"pct_contributions" – dict mapping asset name to its percentage contribution to total loss.
"worst_asset" – name of the asset contributing the most loss.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "AAPL": np.random.normal(0.001, 0.02, 100),
...     "MSFT": np.random.normal(0.001, 0.015, 100),
... })
>>> weights = np.array([0.6, 0.4])
>>> scenario = {"AAPL": -0.15, "MSFT": -0.05}
>>> result = marginal_stress_contribution(weights, returns, scenario)
>>> result["worst_asset"]
'AAPL'

See also

joint_stress_test: Generate stressed parameters for scenarios. wraquant.risk.portfolio.risk_contribution: Euler risk

decomposition (non-scenario-based).

correlation_stress(returns, shock_levels=None)[source]¶

Stress test portfolio by increasing pairwise correlations.

When to use:

Use correlation stress for: - Evaluating diversification fragility: how much does portfolio

risk increase if diversification breaks down?

Regulatory stress testing: correlation breakdown is a standard CCAR scenario.
Risk committee presentations: “if all correlations jump to 0.8, our portfolio vol goes from X% to Y%.”

Parameters:

returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
shock_levels (list[float] | None, default: None) – Blend factors toward perfect correlation. 0 = unchanged, 1 = all correlations set to 1.0. Defaults to [0.0, 0.25, 0.5, 0.75, 1.0].

Returns:

"results" – dict mapping shock level to a dict with "portfolio_vol" (equal-weighted portfolio volatility), "avg_correlation" (mean off-diagonal correlation), "stressed_corr" (the stressed correlation matrix).
"base_vol" – equal-weighted portfolio vol with no shock.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "A": np.random.normal(0, 0.01, 252),
...     "B": np.random.normal(0, 0.012, 252),
...     "C": np.random.normal(0, 0.008, 252),
... })
>>> result = correlation_stress(returns, shock_levels=[0.0, 0.5, 1.0])
>>> result["results"][1.0]["portfolio_vol"] >= result["base_vol"]
True

See also

joint_stress_test: Combined vol, spot, and correlation shocks. vol_stress_test: Volatility-only scaling.

liquidity_stress(returns, volumes=None, liquidity_haircuts=None, portfolio_value=1000000.0)[source]¶

Estimate liquidation cost under adverse market conditions.

Models the cost of unwinding a portfolio under stressed liquidity conditions. If volume data is provided, uses a market-impact model; otherwise, applies user-defined haircuts to each asset.

When to use:: Use liquidity stress for: - Estimating portfolio liquidation costs during crises. - Measuring liquidity-adjusted VaR (LVaR). - Satisfying regulatory requirements for liquidity stress testing

(e.g., SEC Rule 22e-4 for mutual funds).

Parameters:

returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
volumes (DataFrame | None, default: None) – Optional DataFrame of trading volumes (same shape and index as returns). If provided, liquidity cost is estimated as spread * sqrt(position / ADV).
liquidity_haircuts (dict[str, float] | None, default: None) – Optional dict mapping asset name to liquidation cost (e.g., {"AAPL": 0.001, "ILLIQ": 0.05}). If neither volumes nor haircuts are provided, uses volatility as a proxy.
portfolio_value (float, default: 1000000.0) – Total portfolio value for position sizing.

Returns:

"total_cost" – Estimated total liquidation cost ($).
"total_cost_pct" – Cost as a fraction of portfolio value.
"asset_costs" – dict mapping asset to its liquidation cost.
"days_to_liquidate" – estimated days to liquidate if limited to 10% of ADV per day (only if volumes provided).

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "A": np.random.normal(0, 0.01, 252),
...     "B": np.random.normal(0, 0.02, 252),
... })
>>> result = liquidity_stress(returns, portfolio_value=1_000_000)
>>> result["total_cost"] > 0
True

See also

vol_stress_test: Volatility scaling stress test. wraquant.execution.cost: Transaction cost modeling.

scenario_library(returns, scenarios=None)[source]¶

Apply pre-defined crisis scenarios from the built-in library.

When to use:

Use the scenario library for: - Quick stress testing without designing custom scenarios. - Regulatory reporting: standard scenarios that regulators

expect to see.

Benchmarking: compare your portfolio’s sensitivity to well-known crises.

Available scenarios:

"gfc_2008" – Global Financial Crisis
"covid_2020" – COVID-19 crash
"dot_com_2000" – Dot-com bubble burst
"rate_hike_2022" – 2022 rate hiking cycle
"stagflation" – Stagflation scenario
"flash_crash" – Flash crash (intraday)
"em_crisis" – Emerging markets crisis

Parameters:

returns (DataFrame) – Multi-asset return DataFrame (columns = assets).
scenarios (list[str] | None, default: None) – List of scenario names from the library. If None, all scenarios are applied.

Returns:

"scenario_results" – dict mapping scenario name to a dict with "stressed_portfolio_return" (equity shock applied), "stressed_vol" (vol-scaled portfolio volatility), "scenario_params" (the raw scenario parameters).
"available_scenarios" – list of all available scenario names.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> returns = pd.DataFrame({
...     "SPY": np.random.normal(0.0005, 0.01, 252),
...     "TLT": np.random.normal(0.0002, 0.005, 252),
... })
>>> result = scenario_library(returns, scenarios=["gfc_2008", "covid_2020"])
>>> "gfc_2008" in result["scenario_results"]
True

See also

historical_stress_test: Replay actual crisis returns. joint_stress_test: Custom multi-factor stress test.

Scenarios¶

Scenario analysis and Monte Carlo simulation for portfolio risk.

monte_carlo_var(returns, weights, n_sims=10000, confidence=0.95)[source]¶

Estimate portfolio VaR via Monte Carlo simulation.

Draws from a multivariate normal distribution fitted to historical returns.

Parameters:

returns (DataFrame) – DataFrame of asset returns (columns = assets).
weights (ndarray) – Portfolio weight vector.
n_sims (int, default: 10000) – Number of simulation draws.
confidence (float, default: 0.95) – Confidence level.

Return type:

Returns:

VaR as a positive float.

stress_test(returns, weights, shocks)[source]¶

Compute the portfolio loss under a stress scenario.

Each entry in shocks maps an asset name to a shocked return. Assets not in shocks use their historical mean.

Parameters:

returns (DataFrame) – DataFrame of asset returns (columns = assets).
weights (ndarray) – Portfolio weight vector.
shocks (dict[str, float]) – Mapping of asset name to shocked return value.

Return type:

Returns:

Portfolio return under the stress scenario (typically negative).

Historical Events¶

Historical crisis analysis and drawdown attribution.

This module provides tools for analysing portfolio behaviour during historical crises, measuring event impacts, quantifying contagion, and attributing drawdowns to individual assets.

These functions complement the stress testing module (risk.stress) by focusing on what actually happened rather than hypothetical scenarios. Use them for:

Post-mortem analysis: understand what drove past losses.
Regime-aware portfolio construction: identify assets that provide protection in crises.
Contagion monitoring: detect when correlations spike during stress.
Investor reporting: show drawdown history with recovery timelines.

References

Forbes & Rigobon (2002), “No Contagion, Only Interdependence”
Bacon (2008), “Practical Portfolio Performance Measurement and Attribution”

crisis_drawdowns(returns, top_n=5)[source]¶

Identify the top N drawdowns with full lifecycle metrics.

Scans the return series for the largest peak-to-trough drawdowns and reports start date, trough date, recovery date, duration, and magnitude for each.

When to use:

Use crisis drawdowns for: - Investor reporting: show the worst historical losses and

recovery times.

Strategy evaluation: compare drawdown profiles across strategies.
Risk limit calibration: set max drawdown limits based on historical experience.

Parameters:

returns (Series) – Simple return series with a DatetimeIndex (or integer index).
top_n (int, default: 5) – Number of largest drawdowns to return.

Returns:

start – Date the drawdown began (peak).
trough – Date of maximum drawdown.
end – Date the drawdown recovered (or last date if still in drawdown).
drawdown – Magnitude of drawdown (negative number).
days_to_trough – Trading days from start to trough.
days_to_recovery – Trading days from trough to recovery (NaN if not recovered).
total_days – Total drawdown duration.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> idx = pd.bdate_range("2020-01-01", periods=500)
>>> returns = pd.Series(np.random.normal(0.0003, 0.01, 500), index=idx)
>>> dd = crisis_drawdowns(returns, top_n=3)
>>> len(dd) <= 3
True

See also

drawdown_attribution: Which assets caused the drawdowns. wraquant.risk.metrics.max_drawdown: Single worst drawdown.

event_impact(returns, event_dates, window=10)[source]¶

Measure portfolio returns around specific events.

For each event date, extracts the returns in a window before and after the event and computes cumulative return, max drawdown, and volatility within each window.

When to use:

Use event impact analysis for: - Post-mortem: “how did the portfolio react to the Fed rate hike?” - Event studies: systematic analysis of recurring events

(earnings, FOMC, NFP).

Scenario planning: calibrate stress scenarios based on actual event impacts.

Parameters:

returns (Series) – Return series with a DatetimeIndex.
event_dates (list[str]) – List of event date strings (ISO format, e.g., “2020-03-16”). Dates not in the index are matched to the nearest available date.
window (int, default: 10) – Number of trading days before and after the event to analyse.

Returns:

pre_cumulative (float) – Cumulative return in the window before the event.
post_cumulative (float) – Cumulative return in the window after the event.
event_day_return (float) – Return on the event day itself.
pre_vol (float) – Volatility in the pre-event window.
post_vol (float) – Volatility in the post-event window.
total_impact (float) – Cumulative return over the full window (pre + event + post).

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> idx = pd.bdate_range("2020-01-01", periods=252)
>>> returns = pd.Series(np.random.normal(0.0005, 0.01, 252), index=idx)
>>> result = event_impact(returns, ["2020-03-16", "2020-06-15"], window=5)
>>> len(result) >= 1
True

See also

wraquant.risk.stress.historical_stress_test: Replay known crises. crisis_drawdowns: Top drawdown periods.

contagion_analysis(returns_df, crisis_dates)[source]¶

Compare normal vs. crisis-period correlations to detect contagion.

When to use:

Use contagion analysis for: - Evaluating diversification reliability: do correlations spike

when you need diversification most?

Stress testing: adjust portfolio correlations based on empirically observed crisis behaviour.
Regime-aware portfolio construction: allocate less to assets that become highly correlated during crises.

Parameters:

returns_df (DataFrame) – Multi-asset return DataFrame with DatetimeIndex.
crisis_dates (tuple[str, str]) – Tuple of (start_date, end_date) strings defining the crisis period.

Returns:

normal_corr (pd.DataFrame) – Correlation matrix during non-crisis period.
crisis_corr (pd.DataFrame) – Correlation matrix during the crisis period.
corr_change (pd.DataFrame) – Change in correlation (crisis - normal).
avg_normal_corr (float) – Average off-diagonal correlation in normal period.
avg_crisis_corr (float) – Average off-diagonal correlation in crisis period.
contagion_detected (bool) – True if average crisis correlation significantly exceeds normal.
n_normal (int) – Number of normal-period observations.
n_crisis (int) – Number of crisis-period observations.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> idx = pd.bdate_range("2019-01-01", periods=500)
>>> returns = pd.DataFrame({
...     "A": np.random.normal(0.0005, 0.01, 500),
...     "B": np.random.normal(0.0003, 0.012, 500),
... }, index=idx)
>>> result = contagion_analysis(returns, ("2020-02-01", "2020-06-01"))
>>> "contagion_detected" in result
True

See also

wraquant.risk.stress.joint_stress_test: Apply correlation shocks.

References

Forbes & Rigobon (2002), “No Contagion, Only Interdependence: Measuring Stock Market Comovements”

drawdown_attribution(returns_df, weights)[source]¶

Attribute portfolio drawdowns to individual asset contributions.

For each point in the drawdown, decomposes the portfolio’s loss from peak into per-asset contributions. This shows which assets are responsible for the drawdown at each point in time.

When to use:

Use drawdown attribution for: - Post-mortem analysis: “which position caused the 2020 drawdown?” - Risk monitoring: track per-asset drawdown contributions in

real time.

Portfolio construction: identify assets that consistently contribute to drawdowns and consider hedging or removing them.

Parameters:

returns_df (DataFrame) – Multi-asset return DataFrame (columns = assets).
weights (ndarray) – Portfolio weight vector aligned with columns.

Returns:

portfolio_dd – Total portfolio drawdown at each point.
One column per asset showing that asset’s contribution to the drawdown.

Return type:

Example

>>> import pandas as pd, numpy as np
>>> np.random.seed(42)
>>> idx = pd.bdate_range("2020-01-01", periods=252)
>>> returns = pd.DataFrame({
...     "A": np.random.normal(0.0005, 0.01, 252),
...     "B": np.random.normal(0.0003, 0.015, 252),
... }, index=idx)
>>> weights = np.array([0.6, 0.4])
>>> attr = drawdown_attribution(returns, weights)
>>> "portfolio_dd" in attr.columns
True

See also

crisis_drawdowns: Identify top drawdown periods. wraquant.risk.stress.marginal_stress_contribution: Stress-based

attribution.

Credit Risk¶

Credit risk models and default probability estimation.

Credit risk is the risk that a borrower fails to meet its obligations. This module provides tools spanning three major approaches:

Structural models – model the firm’s equity as a contingent claim on its assets. Default occurs when asset value falls below the debt barrier.
- merton_model: the foundational structural model (Merton 1974). Treats equity as a European call option on total firm assets. Iteratively solves for implied asset value and volatility, then computes distance-to-default and default probability.
Credit scoring – statistical models that predict default from accounting ratios or market data.
- altman_z_score: the original 1968 Altman Z-Score for publicly traded manufacturing firms. Combines five accounting ratios into a single score that classifies firms as “safe” (Z > 2.99), “grey zone” (1.81-2.99), or “distress” (Z < 1.81).
Reduced-form / intensity models – model default as a random event driven by a hazard rate (default intensity).
- default_probability: cumulative default probability from a rating transition matrix raised to the power of the horizon.
- credit_spread: implied spread from PD and recovery rate.
- cds_spread: fair CDS premium from a constant hazard rate, integrating protection and premium legs.
- loss_given_default: LGD = exposure * (1 - recovery rate).
- expected_loss: EL = PD * LGD * EAD – the central formula of regulatory capital calculation (Basel II IRB).

How to choose:

For public equities with observable stock prices: merton_model gives market-implied default probabilities that update daily.
For quick screening of financial health: altman_z_score using balance sheet data.
For pricing CDS or credit-linked instruments: cds_spread with calibrated hazard rates.
For portfolio credit risk (e.g., a loan book): default_probability from rating agency transition matrices + expected_loss.

References

Merton (1974), “On the Pricing of Corporate Debt”
Altman (1968), “Financial Ratios, Discriminant Analysis and the Prediction of Corporate Bankruptcy”
Lando (2004), “Credit Risk Modeling: Theory and Applications”

altman_z_score(working_capital, total_assets, retained_earnings, ebit, market_cap, total_liabilities, sales)[source]¶

Altman Z-Score for bankruptcy prediction.

Interpretation:

Z > 2.99 (“safe”): Firm is financially healthy. Default probability is very low (< 1% over 2 years).
1.81 <= Z <= 2.99 (“grey zone”): Ambiguous. Firm could go either way. Warrants deeper analysis.
Z < 1.81 (“distress”): High bankruptcy risk. Historically, ~95% of firms that defaulted had Z < 1.81 one year prior.

Component interpretation:

x1 (WC/TA): Liquidity. Negative = current liabilities exceed current assets.
x2 (RE/TA): Cumulative profitability and firm age. Young firms have low retained earnings.
x3 (EBIT/TA): Operating profitability.
x4 (Market Cap/TL): Market leverage.
x5 (Sales/TA): Asset turnover efficiency.

When to use:

Quick screening of a large universe of firms.
As a factor in multi-factor credit models.
For early warning systems.

Limitations:

Designed for publicly traded manufacturing firms.
Does not capture market-implied information (use merton_model for that).
Accounting data can be manipulated.

Parameters:

working_capital (float) – Current assets minus current liabilities.
total_assets (float) – Total assets.
retained_earnings (float) – Cumulative retained earnings.
ebit (float) – Earnings before interest and taxes.
market_cap (float) – Market capitalisation of equity.
total_liabilities (float) – Total liabilities.
sales (float) – Net sales / revenue.

Returns:

z_score: The computed Z-Score.
zone: One of "safe" (Z > 2.99), "grey" (1.81 <= Z <= 2.99), or "distress" (Z < 1.81).
x1 .. x5: Individual component ratios.

Return type:

dict[str, float | str]

cds_spread(default_intensity, recovery_rate, maturity)[source]¶

Fair CDS spread from a constant hazard rate (default intensity).

Computes the breakeven CDS premium by equating the expected protection leg (what the protection seller pays at default) with the expected premium leg (what the protection buyer pays over time).

Under a constant hazard rate model, the fair spread is approximately lambda * (1 - R), but this function uses the exact continuous-time formula with quarterly premium payments for greater accuracy.

Interpretation:

The output is the annualized spread (decimal). Multiply by 10,000 for basis points.
Compare to market CDS spreads to detect relative value.
If model spread > market spread: protection is cheap (market underestimates default risk).
If model spread < market spread: protection is expensive or there is a risk premium.
CDS spreads are approximately equal to bond spreads over swaps (CDS-bond basis ~ 0 in normal markets).

When to use:

Pricing CDS contracts given a calibrated hazard rate.
Calibrating hazard rates from market CDS spreads (invert numerically).
Converting between PD and spread for credit analysis.

Parameters:

default_intensity (float) – Constant hazard rate (lambda), annualized. E.g., 0.02 means a 2% probability of default per year.
recovery_rate (float) – Recovery rate in [0, 1]. Standard assumption is 0.40 for senior unsecured corporate debt.
maturity (float) – CDS maturity in years (standard: 1, 3, 5, 7, 10).

Return type:

Returns:

Annualized CDS spread as a decimal fraction. Multiply by 10,000 for basis points.

Example

>>> spread = cds_spread(0.02, 0.40, 5.0)
>>> print(f"5Y CDS: {spread*10000:.0f} bps")  # ~120 bps

credit_spread(default_prob, recovery_rate, rf_rate=0.0)[source]¶

Implied credit spread from a default probability.

The formula: spread = -ln(1 - PD * LGD), where LGD = 1 - R. For small PD, this simplifies to spread ~ PD * LGD.

Interpretation:

Output is annualized as a decimal (0.01 = 100 bps).
Multiply by 10,000 for basis points.
If market spread > model spread: bond is cheap (excess compensation for credit risk).
If market spread < model spread: bond is expensive or the model PD is too high.

Parameters:

default_prob (float) – Annualized probability of default (e.g., 0.02 for 2% annual PD).
recovery_rate (float) – Recovery rate in [0, 1]. Investment grade typically 0.40-0.50; high yield 0.25-0.40.
rf_rate (float, default: 0.0) – Risk-free rate (unused in simple model but accepted for API consistency).

Return type:

Returns:

Annualized credit spread as a decimal fraction. Multiply by 10,000 for basis points.

Example

>>> spread = credit_spread(0.02, 0.40)
>>> print(f"Spread: {spread*10000:.0f} bps")  # ~120 bps

default_probability(rating_transitions, horizon)[source]¶

Cumulative default probability from a rating transition matrix.

Interpretation:

The output is a vector where each element is the cumulative default probability for a given starting rating.
AAA will have the smallest PD; CCC the largest.
Compare to historical default rates published by Moody’s or S&P to calibrate.
Use with expected_loss for portfolio credit risk.

When to use:

Converting a rating agency transition matrix into PDs for capital calculations.
Stress testing: modify the transition matrix (increase downgrade probabilities) and re-compute PDs.

The last row/column of the transition matrix is assumed to represent the default (absorbing) state.

Parameters:

rating_transitions (ndarray) – Square transition matrix of shape (n, n) where element [i, j] is the one-period probability of migrating from rating i to rating j. Rows must sum to 1.
horizon (int) – Number of periods to compound over (e.g., 5 for 5-year cumulative PD).

Return type:

Returns:

1-D array of cumulative default probabilities for each non-default rating, length n - 1.

Example

>>> import numpy as np
>>> # Simple 3-state matrix: AAA, BBB, Default
>>> T = np.array([[0.95, 0.04, 0.01],
...               [0.02, 0.90, 0.08],
...               [0.00, 0.00, 1.00]])
>>> pd_5yr = default_probability(T, horizon=5)
>>> print(f"5yr PD from AAA: {pd_5yr[0]:.4f}")
>>> print(f"5yr PD from BBB: {pd_5yr[1]:.4f}")

expected_loss(pd_val, lgd, ead)[source]¶

Expected loss (EL = PD x LGD x EAD).

The expected loss is the central formula of credit risk management and Basel II/III regulatory capital calculation. It represents the average loss you expect from a credit exposure.

Interpretation:

EL is the mean of the loss distribution. It should be covered by pricing (loan margins, bond spreads) rather than capital reserves.
Capital reserves cover the unexpected loss (UL), which is the tail beyond EL.
For a portfolio, EL is additive: sum over all exposures.

Parameters:

pd_val (float) – Probability of default (annualized, e.g., 0.02 for 2%).
lgd (float) – Loss given default as a fraction of EAD (e.g., 0.45 for 45% loss rate).
ead (float) – Exposure at default (dollar amount at risk).

Return type:

Returns:

Expected loss in the same units as ead.

Example

>>> el = expected_loss(pd_val=0.02, lgd=0.45, ead=1_000_000)
>>> print(f"Expected loss: ${el:,.0f}")  # $9,000

loss_given_default(exposure, recovery_rate)[source]¶

Expected loss given default.

LGD = EAD * (1 - Recovery Rate). This is the dollar amount you expect to lose if the borrower defaults.

Interpretation:

A recovery rate of 0.40 means you recover 40 cents on the dollar; LGD is 60% of exposure.
Recovery rates vary by seniority: secured senior ~65%, unsecured senior ~45%, subordinated ~25%.
Use with expected_loss for Basel II/III capital calculations.

Parameters:

exposure (float) – Exposure at default (EAD) – the amount at risk.
recovery_rate (float) – Expected recovery rate in [0, 1].

Return type:

Returns:

Loss given default = exposure * (1 - recovery_rate).

merton_model(equity, debt, vol, rf_rate, maturity)[source]¶

Merton structural credit risk model.

The model iteratively solves for the unobservable asset value and asset volatility from the observable equity value and equity volatility, using the Black-Scholes option pricing relationship.

Interpretation:

distance_to_default (DD): How many standard deviations the firm’s asset value is above the default barrier. DD > 4: very safe. DD 2-4: investment grade. DD 1-2: high yield. DD < 1: distress. Moody’s KMV uses this as the primary input to their EDF (Expected Default Frequency) model.
default_probability: N(-DD), the probability that asset value drifts below debt by maturity. This is a risk-neutral probability – real-world default rates are typically lower.
asset_vol: Implied asset volatility. Higher = more default risk. Asset vol is always lower than equity vol because equity is a levered claim.
credit_spread: The yield premium investors should demand for holding risky debt. Compare to market CDS spreads to detect mispricing.

When to use:

Market-implied default probabilities from daily equity data.
Screening for distressed firms (DD < 2).
As a factor in credit scoring models.
For relative value: compare Merton-implied spreads to market CDS spreads.

Red flags:

DD < 1: firm is in acute distress.
Asset vol > 50%: inputs may be unreliable.
Equity vol is stale or missing: model won’t converge.

Parameters:

equity (float) – Current market value of equity (market cap).
debt (float) – Face value of outstanding debt (the “strike”).
vol (float) – Equity volatility (annualized, e.g., 0.30 for 30%).
rf_rate (float) – Continuous risk-free rate (annualized).
maturity (float) – Time to maturity of debt in years (typically 1).

Returns:

asset_value (float) – Implied total asset value.
asset_vol (float) – Implied asset volatility.
d1, d2 (float) – Black-Scholes d1 and d2.
distance_to_default (float) – d2, number of std devs above the default barrier.
default_probability (float) – N(-d2), risk-neutral probability of default.
credit_spread (float) – Implied credit spread over the risk-free rate (annualized).

Return type:

Example

>>> result = merton_model(equity=50e6, debt=40e6, vol=0.35,
...                       rf_rate=0.04, maturity=1.0)
>>> print(f"DD: {result['distance_to_default']:.2f}")
>>> print(f"PD: {result['default_probability']:.4f}")
>>> print(f"Spread: {result['credit_spread']*10000:.0f} bps")

See also

altman_z_score: Accounting-based bankruptcy predictor. cds_spread: Reduced-form CDS pricing.

Notes

Reference: Merton, R.C. (1974). “On the Pricing of Corporate Debt: The Risk Structure of Interest Rates.” Journal of Finance, 29(2), 449-470.

Survival Analysis¶

Survival analysis estimators for financial applications.

Survival analysis models the time until an event occurs – in finance, this could be time-to-default, time until a drawdown ends, time between large losses, or fund lifetime. The key challenge is censoring: not all subjects experience the event during the observation period.

This module provides pure numpy/scipy implementations of the standard survival analysis toolkit:

Non-parametric estimators:

kaplan_meier: the Kaplan-Meier product-limit estimator of the survival function S(t) = P(T > t). The most common survival curve estimator. Handles right-censored data.
nelson_aalen: cumulative hazard estimator H(t). Related to Kaplan-Meier via S(t) = exp(-H(t)) but more natural for hazard rate estimation.
hazard_rate: kernel-smoothed instantaneous hazard rate from Nelson-Aalen increments. Useful for visualising how default risk changes over time.

Semi-parametric model:

cox_partial_likelihood: Cox proportional hazards model. Estimates the effect of covariates on the hazard rate without specifying the baseline hazard. The workhorse of survival regression: “does leverage, size, or profitability affect time-to-default?”

Parametric models:

exponential_survival: S(t) = exp(-lambda * t). Assumes constant hazard (memoryless). Simple but often too restrictive.
weibull_survival: S(t) = exp(-(t/lambda)^k). Generalises exponential (k=1). k < 1 = decreasing hazard (burn-in), k > 1 = increasing hazard (aging/wear-out), which corresponds to increasing default risk with time for distressed firms.

Hypothesis testing:

log_rank_test: compares two survival curves. Use to test whether two groups (e.g., investment-grade vs. high-yield) have significantly different survival distributions.

Utility:

median_survival_time: smallest t where S(t) <= 0.5.

Financial applications:

Credit risk: model time-to-default with Cox PH, using leverage, profitability, and market indicators as covariates.
Drawdown analysis: model time to recover from a drawdown; Weibull shape > 1 suggests recovery becomes less likely over time.
Fund closure: Kaplan-Meier curves for hedge fund lifetimes, stratified by strategy type.
Trade duration: model how long a position is held before hitting a stop-loss or take-profit.

References

Cox (1972), “Regression Models and Life-Tables”
Kaplan & Meier (1958), “Nonparametric Estimation from Incomplete Observations”
Lando (2004), “Credit Risk Modeling: Theory and Applications”

cox_partial_likelihood(durations, event_observed, covariates, max_iter=100, tol=1e-09)[source]¶

Cox proportional hazards model via Newton-Raphson.

In finance: “Does leverage, profitability, or market beta affect the hazard of default, controlling for other factors?”

Interpretation:

beta[j] is the log hazard ratio for covariate j. exp(beta[j]) > 1 means the covariate increases the hazard (bad for survival). exp(beta[j]) < 1 means it decreases the hazard (protective).
se[j]: Standard error. beta[j] / se[j] gives a z-statistic. |z| > 1.96 is significant at the 5% level.
log_partial_likelihood: Higher (less negative) = better fit. Use for comparing nested models via likelihood ratio tests.

Red flags:

Very large beta (|beta| > 5): possible separation/convergence issues.
n_iter = max_iter: did not converge, results unreliable.
se contains NaN: Hessian is singular, model is degenerate.

Parameters:

durations (ndarray) – 1-D array of observed durations.
event_observed (ndarray) – 1-D boolean/int array indicating event occurrence.
covariates (ndarray) – 2-D array of shape (n_subjects, n_covariates).
max_iter (int, default: 100) – Maximum Newton-Raphson iterations.
tol (float, default: 1e-09) – Convergence tolerance on the gradient norm.

Returns:

beta (ndarray) – Regression coefficients. exp(beta) gives hazard ratios.
se (ndarray) – Standard errors of coefficients.
log_partial_likelihood (float) – Maximised log partial likelihood.
n_iter (int) – Number of iterations to convergence.

Return type:

dict[str, ndarray | float]

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> n = 200
>>> leverage = rng.uniform(0.2, 0.8, n)
>>> durations = rng.exponential(5 / (1 + leverage), n)
>>> events = np.ones(n, dtype=bool)
>>> result = cox_partial_likelihood(durations, events, leverage.reshape(-1, 1))
>>> print(f"Leverage HR: {np.exp(result['beta'][0]):.2f}")

See also

kaplan_meier: Non-parametric survival curve (no covariates). weibull_survival: Parametric survival model.

exponential_survival(lambda_param, t)[source]¶

Exponential survival function S(t) = exp(-lambda * t).

Interpretation:

lambda = 0.1 means roughly a 10% chance of the event per unit of time.
Mean survival time = 1 / lambda.
If the hazard is actually increasing or decreasing over time, this model is too simplistic. Use weibull_survival instead.

Parameters:

lambda_param (float) – Hazard rate (constant, > 0). Higher = faster time to event.
t (float | ndarray) – Time point(s) at which to evaluate.

Return type:

Returns:

Survival probability at each t.

See also

weibull_survival: Generalises exponential with time-varying hazard.

hazard_rate(durations, event_observed, bandwidth=None)[source]¶

Kernel-smoothed hazard rate estimate.

Applies Epanechnikov kernel smoothing to the Nelson-Aalen increments to produce a smooth instantaneous hazard rate function h(t).

Interpretation:

A flat hazard rate means constant risk (exponential model).
An increasing hazard means risk accelerates with time (typical for credit deterioration, infrastructure aging).
A decreasing hazard means early failures dominate and survivors become stronger (“infant mortality”).
A bathtub-shaped hazard (decreasing then increasing) is common in reliability engineering.

Parameters:

durations (ndarray) – 1-D array of observed durations.
event_observed (ndarray) – 1-D boolean/int array indicating event occurrence.
bandwidth (float | None, default: None) – Kernel bandwidth. If None, uses Silverman’s rule. Larger bandwidth = smoother curve. Smaller = more detail but noisier.

Returns:

timeline (ndarray) – Evaluation grid.
hazard (ndarray) – Smoothed hazard rate at each point.

Return type:

See also

nelson_aalen: The cumulative hazard from which this is derived. kaplan_meier: Survival function estimator.

kaplan_meier(durations, event_observed)[source]¶

Kaplan-Meier survival curve estimator.

In finance, this answers questions like:

“What fraction of bonds survive to year 5 without defaulting?”
“How long do hedge funds typically survive before closing?”
“What is the probability a drawdown lasts longer than 60 days?”

Interpretation:

survival: S(t) is the probability of surviving beyond time t. A steep drop indicates a period of high hazard.
variance (Greenwood’s formula): Use to construct 95% confidence bands as S(t) +/- 1.96 * sqrt(variance(t)).
The median survival time is where S(t) first drops below 0.5.
A flat survival curve = low hazard rate (few events).
A curve that drops quickly early = high initial hazard (e.g., new funds failing in the first year).

Parameters:

durations (ndarray) – 1-D array of observed durations (time-to-event or time-to-censoring).
event_observed (ndarray) – 1-D boolean/int array where 1 (True) indicates the event occurred and 0 (False) indicates right-censoring.

Returns:

timeline (ndarray) – Sorted unique event times.
survival (ndarray) – Survival probability S(t) at each event time.
variance (ndarray) – Greenwood’s variance estimate for constructing confidence bands.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> durations = rng.exponential(5.0, 200)  # avg survival 5 years
>>> events = rng.binomial(1, 0.7, 200)     # 70% observed, 30% censored
>>> km = kaplan_meier(durations, events)
>>> print(f"5-year survival: {km['survival'][km['timeline'] <= 5][-1]:.2f}")

See also

nelson_aalen: Cumulative hazard estimator. log_rank_test: Compare two survival curves.

log_rank_test(durations1, event1, durations2, event2)[source]¶

Log-rank test comparing two survival curves.

Tests the null hypothesis that the survival functions of two groups are identical. This is the standard test for comparing survival experiences between groups.

In finance: “Do investment-grade bonds have significantly different time-to-default than high-yield bonds?” or “Do value stocks have different drawdown durations than growth stocks?”

Interpretation:

p_value < 0.05: reject H0 – the two groups have significantly different survival experiences.
observed1 >> expected1: Group 1 has more events than expected (worse survival).
observed1 << expected1: Group 1 has fewer events than expected (better survival).
The test is most powerful when the hazard ratio is constant (proportional hazards). For crossing survival curves, the Wilcoxon (Breslow) test may be more appropriate.

Parameters:

durations1 (ndarray) – Durations for group 1.
event1 (ndarray) – Event indicators for group 1.
durations2 (ndarray) – Durations for group 2.
event2 (ndarray) – Event indicators for group 2.

Returns:

test_statistic (float) – Chi-squared statistic (1 df).
p_value (float) – P-value. < 0.05 rejects equality.
observed1 (float) – Total observed events in group 1.
expected1 (float) – Expected events under H0.

Return type:

Example

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> d1 = rng.exponential(5.0, 100)  # group 1: avg survival 5y
>>> d2 = rng.exponential(3.0, 100)  # group 2: avg survival 3y
>>> e1 = np.ones(100, dtype=bool)
>>> e2 = np.ones(100, dtype=bool)
>>> result = log_rank_test(d1, e1, d2, e2)
>>> print(f"p-value: {result['p_value']:.4f}")  # should be small

median_survival_time(durations, event_observed)[source]¶

Median survival time from the Kaplan-Meier estimator.

The median survival time is the smallest time t at which the estimated survival function drops to or below 0.5 – i.e., the time by which half the subjects have experienced the event.

Interpretation:

This is the “half-life” of the population.
More robust than mean survival (which is heavily influenced by censoring and long survivors).
Returns np.inf if the survival curve never reaches 0.5, which happens with heavy censoring or if more than half the subjects never experience the event.

In finance:

“The median time-to-default for CCC-rated firms is 2.3 years.”
“The median drawdown recovery time for equity portfolios is 45 trading days.”

Parameters:

durations (ndarray) – 1-D array of observed durations.
event_observed (ndarray) – 1-D boolean/int array indicating event occurrence.

Return type:

Returns:

Median survival time, or np.inf if the survival curve never drops to 0.5 (too much censoring).

nelson_aalen(durations, event_observed)[source]¶

Nelson-Aalen cumulative hazard estimator.

Interpretation:

H(t) represents the accumulated risk up to time t.
The slope of H(t) is the instantaneous hazard rate: steep segments indicate periods of high risk.
A linear H(t) suggests a constant hazard rate (exponential survival).
A concave H(t) suggests a decreasing hazard (survival gets easier over time).
A convex H(t) suggests an increasing hazard (risk accelerates – typical for aging/wear-out or credit deterioration).

Parameters:

durations (ndarray) – 1-D array of observed durations.
event_observed (ndarray) – 1-D boolean/int array indicating event occurrence.

Returns:

timeline (ndarray) – Sorted unique event times.
cumulative_hazard (ndarray) – H(t) at each time.
variance (ndarray) – Variance estimate at each time.

Return type:

See also

kaplan_meier: Direct survival function estimator. hazard_rate: Smoothed instantaneous hazard from Nelson-Aalen.

weibull_survival(lambda_param, k, t)[source]¶

Weibull survival function S(t) = exp(-(t / lambda)^k).

The Weibull distribution generalises the exponential by allowing the hazard rate to increase or decrease over time. It is the most commonly used parametric survival model in practice.

Interpretation of the shape parameter k:

k = 1: Constant hazard (reduces to exponential). The event is equally likely at any time.
k < 1: Decreasing hazard (“burn-in”). Early failures are most common; survivors become stronger. Typical for infant mortality in manufactured goods.
k > 1: Increasing hazard (“aging”). Risk increases over time. Typical for credit deterioration in distressed firms or aging infrastructure.

In finance:

k > 1 for time-to-default: firms that have survived a long time in distress become more likely to default (debt maturity approaches, liquidity dries up).
k < 1 for drawdown recovery: if a drawdown has already lasted a long time, recovery becomes more likely (mean reversion kicks in).

Parameters:

lambda_param (float) – Scale parameter (> 0). Larger = longer survival.
k (float) – Shape parameter (> 0). k=1 is exponential. k>1 is increasing hazard. k<1 is decreasing hazard.
t (float | ndarray) – Time point(s) at which to evaluate.

Return type:

Returns:

Survival probability at each t.

See also

exponential_survival: Simplest case (k=1). kaplan_meier: Non-parametric alternative.

Integrations¶

Advanced risk and portfolio optimisation integrations.

Provides wrappers around PyPortfolioOpt, Riskfolio-Lib, skfolio, copulas, pyvinecopulib, and pyextremes for portfolio construction, copula modelling, and extreme value analysis.

pypfopt_efficient_frontier(expected_returns, cov_matrix)[source]¶

Compute the efficient frontier using PyPortfolioOpt.

Solves for the maximum-Sharpe-ratio portfolio and returns the optimal weights along with performance metrics.

Parameters:

expected_returns (Series) – Expected annual returns per asset.
cov_matrix (DataFrame) – Annualised covariance matrix of asset returns.

Returns:

Dictionary containing:

weights – dict mapping asset names to optimal weights.
expected_return – portfolio expected annual return.
volatility – portfolio expected annual volatility.
sharpe_ratio – portfolio Sharpe ratio.

Return type:

riskfolio_portfolio(returns, method='MV')[source]¶

Optimise a portfolio using Riskfolio-Lib.

Parameters:

returns (DataFrame) – Historical return data. Each column is an asset.
method (str, default: 'MV') –
Risk measure for optimisation:
- 'MV' – minimum variance
- 'CVaR' – conditional value-at-risk
- 'MAD' – mean absolute deviation

Returns:

Dictionary containing:

weights – dict mapping asset names to optimal weights.
method – risk measure used.

Return type:

skfolio_optimize(returns, objective='min_variance')[source]¶

Optimise a portfolio using skfolio.

Parameters:

returns (DataFrame) – Historical return data. Each column is an asset.
objective (str, default: 'min_variance') –
Optimisation objective:
- 'min_variance' – minimum variance portfolio
- 'max_sharpe' – maximum Sharpe ratio portfolio

Returns:

Dictionary containing:

weights – dict mapping asset names to optimal weights.
objective – objective used.

Return type:

copulas_fit(data, copula_type='gaussian')[source]¶

Fit a copula model to multivariate data.

Uses the copulas library to fit a copula and provides methods for sampling and density evaluation.

Parameters:

data (DataFrame) – Multivariate data. Each column is a variable.
copula_type (str, default: 'gaussian') –
Type of copula to fit:
- 'gaussian' – Gaussian copula
- 'vine' – vine copula (via copulas library)

Returns:

Dictionary containing:

copula – fitted copula object.
copula_type – type of copula used.
columns – list of column names from the input data.
n_samples – number of observations used for fitting.

Return type:

copulae_fit(data, family='gaussian')[source]¶

Fit a copula using the copulae library.

Parameters:

data (ndarray | DataFrame) – Data of shape (n, d). Can be raw observations (marginals are automatically converted to pseudo-observations) or pre-transformed uniform marginals on [0, 1].
family (str, default: 'gaussian') –
Copula family to fit:
- 'gaussian' – Gaussian copula (no tail dependence).
- 'student' – Student-t copula (symmetric tail dependence).
- 'clayton' – Clayton copula (lower tail dependence).
- 'gumbel' – Gumbel copula (upper tail dependence).
- 'frank' – Frank copula (symmetric, no tail dependence).

Returns:

Dictionary containing:

params – fitted copula parameters (structure depends on family).
log_likelihood – log-likelihood of the fitted model.
aic – Akaike information criterion.
bic – Bayesian information criterion (approximate).
fitted_copula – the fitted copula object for further use (sampling, CDF evaluation, etc.).

Return type:

Example

>>> import numpy as np
>>> from wraquant.risk.integrations import copulae_fit
>>> rng = np.random.default_rng(42)
>>> data = rng.normal(0, 1, (200, 3))
>>> result = copulae_fit(data, family="gaussian")
>>> result["log_likelihood"]

Notes

Reference: Joe (2014). Dependence Modeling with Copulas. Chapman & Hall/CRC.

See also

copulas_fit: Alternative using the copulas library.
fit_gaussian_copula: Built-in Gaussian copula implementation.
vine_copula: Vine copula fitting via pyvinecopulib.

vine_copula(data, structure='regular')[source]¶

Fit a vine copula using pyvinecopulib.

Parameters:

data (ndarray | DataFrame) – Multivariate data. Each column is a variable. Values should ideally be on the unit interval (pseudo-observations); if not, the data is rank-transformed automatically.
structure (str, default: 'regular') – Vine structure type. Currently only 'regular' is supported.

Returns:

Dictionary containing:

vinecop – fitted pyvinecopulib.Vinecop object.
structure – vine structure used.
n_vars – number of variables.
loglik – log-likelihood of the fitted model.

Return type:

extreme_value_analysis(data)[source]¶

Perform extreme value analysis using pyextremes.

Fits a Generalized Extreme Value (GEV) distribution to block maxima extracted from the data.

Parameters:

data (Series | ndarray) – Univariate time series of observations (e.g. losses or negative returns).

Returns:

Dictionary containing:

shape – GEV shape parameter (xi).
loc – GEV location parameter (mu).
scale – GEV scale parameter (sigma).
return_levels – dict of return levels for common return periods (10, 50, 100 years).

Return type: