Applications¶

Reference implementations for analysis, not production pricing or capital libraries. Each module states its modelling assumptions.

rcopula.risk ¶

Portfolio and systemic risk measures under copula dependence.

The reason copulas matter for risk is narrow and specific: the dependence model changes the tail of the aggregate loss even when every margin and every pairwise correlation is held fixed. Two portfolios calibrated to the same Kendall's tau can differ by tens of percent in 99.9% VaR, because a Gaussian copula says extreme losses arrive independently and a t or Clayton copula says they arrive together. This module makes that difference measurable.

What is here:

=========================== ================================================= :func:value_at_risk The quantile of the loss distribution. :func:expected_shortfall Mean loss beyond VaR -- coherent, unlike VaR. :func:simulate_losses Aggregate portfolio loss under a copula model. :func:risk_contributions Euler allocation of ES to positions. :func:diversification_benefit How much dependence costs you. :func:covar / :func:delta_covar Systemic risk: system VaR given a firm. :func:marginal_expected_shortfall A firm's loss in a system-wide crisis. :func:stress_scenario Conditional simulation given a stressed factor. :func:rank_reorder Impose a copula on fixed marginal samples. =========================== =================================================

Sign convention. Everything here works in losses: positive numbers are bad, and VaR at level alpha is the alpha-quantile of the loss. If you have returns, negate them first.

.. warning::

That convention flips which tail matters, and it is easy to get backwards. Aggregate loss risk is driven by upper tail dependence -- the tendency of large losses to arrive together. Clayton is lower-tail dependent, so applying it directly to losses clusters the small ones and understates risk. Measured on an equally weighted five-name portfolio, all calibrated to the same Kendall's tau of 0.5:

========================== ============= ========== copula 99% ES lambda_U ========================== ============= ========== Clayton (lower tail) 4.91 0.00 Gaussian (no tail dep) 7.16 0.00 Student t, df=4 (both) 7.61 0.40 Gumbel (upper tail) 8.41 0.59 ========================== ============= ==========

The ordering follows lambda_U exactly, and Clayton comes out below Gaussian. "Crashes cluster, so use Clayton" is right about returns and wrong about losses: negating flips the tails. For losses, reach for Gumbel, a Student-t, or a survival (180-degree rotated) Clayton.

References

McNeil, A. J., Frey, R. and Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques and Tools, 2nd ed. Princeton. Chapters 2 and 8 for risk measures and their aggregation. Acerbi, C. and Tasche, D. (2002). On the coherence of expected shortfall. Journal of Banking & Finance 26(7), 1487-1503. Adrian, T. and Brunnermeier, M. K. (2016). CoVaR. American Economic Review 106(7), 1705-1741. Acharya, V. V., Pedersen, L. H., Philippon, T. and Richardson, M. (2017). Measuring systemic risk. Review of Financial Studies 30(1), 2-47. Marginal expected shortfall. Tasche, D. (2007). Euler allocation: theory and practice. arXiv:0708.2542.

RiskSummary ¶

Bases: NamedTuple

A VaR and expected-shortfall pair, both measured at the same confidence level.

A small named tuple, so it unpacks like a plain tuple (alpha, var, es = summary) and also allows attribute access.

Attributes:

Name	Type	Description
`alpha`	`float`	Confidence level in `(0, 1)`, e.g. `0.99` for 99%.
`var`	`float`	Value at Risk at `alpha`, in the same units as the losses: the loss exceeded on only `1 - alpha` of occasions.
`expected_shortfall`	`float`	Expected shortfall at `alpha`, in the same units as the losses: the average loss on those occasions when VaR is reached or exceeded.

Examples:

>>> from rcopula.risk import RiskSummary
>>> s = RiskSummary(0.99, 2.5, 3.1)
>>> s.var
2.5

value_at_risk ¶

value_at_risk(losses: ArrayLike, alpha: float = 0.99) -> float

The loss you should exceed on only 1 - alpha of occasions (Value at Risk).

With alpha=0.99 and daily losses, this is the loss that is beaten on roughly one day in a hundred. Technically it is the alpha-quantile of the loss distribution (the "inverted CDF" quantile, so the result is always one of the input values).

Parameters:

Name	Type	Description	Default
`losses`	`array_like of float, shape (n,)`	Simulated or realised losses. Positive numbers are losses, negative numbers are gains; if you have returns, pass `-returns`. Any shape is accepted and flattened.	required
`alpha`	`float`	Confidence level, strictly between 0 and 1, e.g. `0.99` for 99% VaR.	`0.99`

Returns:

Type	Description
`float`	The VaR, in the same units as `losses`.

Raises:

Type	Description
`ValueError`	If `alpha` is not strictly inside `(0, 1)`, or `losses` is empty.

Notes

VaR is not coherent: it can penalise diversification, because it says nothing about how bad the tail is beyond the quantile. Regulatory frameworks have been moving to expected shortfall for exactly that reason. It is reported here because it is still ubiquitous, not because it is the better measure.

Examples:

>>> import numpy as np
>>> from rcopula.risk import value_at_risk
>>> losses = np.arange(1, 1001, dtype=float)
>>> float(value_at_risk(losses, 0.99))
990.0

expected_shortfall ¶

expected_shortfall(losses: ArrayLike, alpha: float = 0.99) -> float

The average loss on the days when VaR is reached or exceeded (expected shortfall, CVaR).

Where VaR tells you where the tail starts, expected shortfall tells you how bad the tail is on average. Use it as the headline tail-risk number; it is the measure Basel III's market-risk rules (FRTB) use, at alpha=0.975.

Parameters:

Name	Type	Description	Default
`losses`	`array_like of float, shape (n,)`	Simulated or realised losses. Positive numbers are losses, negative numbers are gains; if you have returns, pass `-returns`. Any shape is accepted and flattened.	required
`alpha`	`float`	Confidence level, strictly between 0 and 1.	`0.99`

Returns:

Type	Description
`float`	The expected shortfall, in the same units as `losses`. It is always at least as large as :func:`value_at_risk` at the same `alpha`.

Raises:

Type	Description
`ValueError`	If `alpha` is not strictly inside `(0, 1)`, or `losses` is empty.

Notes

.. math:: \mathrm{ES}\alpha = \mathbb{E}[L \mid L \ge \mathrm{VaR}\alpha].

Unlike VaR this is coherent -- in particular subadditive, so it never penalises diversification -- which is why Basel III replaced 99% VaR with 97.5% ES.

Examples:

>>> import numpy as np
>>> from rcopula.risk import expected_shortfall, value_at_risk
>>> losses = np.arange(1, 1001, dtype=float)
>>> float(expected_shortfall(losses, 0.99))
995.0

ES always dominates VaR at the same level:

>>> rng = np.random.default_rng(0)
>>> x = rng.standard_t(3, size=100_000)
>>> bool(expected_shortfall(x, 0.99) > value_at_risk(x, 0.99))
True

rank_reorder ¶

rank_reorder(samples: ArrayLike, copula: Copula, random_state: Generator | int | None = None) -> NDArray[np.float64]

Make separately simulated risks co-move as a copula says, without changing their values.

Use this when you already have simulated losses for each business line or asset (for example from approved standalone models) and only need to decide how they co-move. Each column keeps exactly the same values; only the order of rows changes, so that the columns' ranks follow a draw from copula. This is the Iman-Conover style "rank reordering" used in capital aggregation.

Parameters:

Name	Type	Description	Default
`samples`	`array_like of float, shape (n, d)`	Marginal samples, one column per risk and one row per scenario. The columns are treated independently; their current row order is discarded.	required
`copula`	`Copula`	The dependence structure to impose. Its `dim` must equal `d`.	required
`random_state`	`(int, Generator or None)`	Seed or generator for the copula draw that supplies the ranks. Pass an int for reproducible results.	`None`

Returns:

Type	Description
`numpy.ndarray of float, shape (n, d)`	The reordered samples. Sorting any column gives back the sorted input column exactly.

Raises:

Type	Description
`ValueError`	If the number of columns differs from `copula.dim`.

Notes

Each column is sorted and then re-ordered to follow the ranks of a draw from copula. The marginal distributions survive exactly -- the same values come back, only rearranged -- while the dependence becomes the copula's.

This is the standard move in risk aggregation when the marginals are already fixed: capital models arrive with an approved loss distribution per business line and the question is only how to combine them. Refitting the margins to make a joint model tractable would change numbers that have been signed off; reordering does not touch them.

Examples:

The margins are preserved to the last value:

>>> import numpy as np
>>> from rcopula import ClaytonCopula
>>> from rcopula.risk import rank_reorder
>>> rng = np.random.default_rng(0)
>>> x = np.column_stack([rng.lognormal(size=5000), rng.exponential(size=5000)])
>>> y = rank_reorder(x, ClaytonCopula(4.0), random_state=0)
>>> bool(np.array_equal(np.sort(x, axis=0), np.sort(y, axis=0)))
True

...while the dependence becomes the copula's:

>>> from scipy import stats
>>> bool(abs(stats.kendalltau(y[:, 0], y[:, 1]).statistic - 2 / 3) < 0.03)
True

simulate_losses ¶

simulate_losses(copula: Copula, margins: Margin | list[Margin], weights: ArrayLike | None = None, n: int = 100000, random_state: Generator | int | None = None) -> NDArray[np.float64]

Simulate total portfolio losses when the positions' losses are tied by a copula.

This is the Monte Carlo engine behind the other functions here: feed its output to :func:value_at_risk or :func:expected_shortfall. Each draw picks a joint scenario from the copula (how the positions co-move), turns each coordinate into a loss using that position's own loss distribution, and adds them up with weights.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Dependence between the individual loss drivers. Its `dim` is the number of positions `d`.	required
`margins`	`scipy.stats frozen distribution or list of them, length d`	Marginal loss distributions, one per position (positive values are losses). A single distribution is reused for every position.	required
`weights`	`array_like of float, shape (d,)`	Position weights or exposures. Defaults to equal weights `1/d` summing to one. Pass exposures directly (e.g. notionals) for an unnormalised total in currency units.	`None`
`n`	`int`	Number of simulated scenarios.	`100_000`
`random_state`	`(int, Generator or None)`	Seed or generator for reproducibility.	`None`

Returns:

Type	Description
`numpy.ndarray of float, shape (n,)`	Simulated aggregate losses, one per scenario.

Raises:

Type	Description
`ValueError`	If `weights` does not have exactly `d` entries, or a list of `margins` has the wrong length.

Examples:

>>> import numpy as np
>>> from scipy import stats
>>> from rcopula import ClaytonCopula
>>> from rcopula.risk import simulate_losses, value_at_risk
>>> losses = simulate_losses(
...     ClaytonCopula(2.0, dim=3), stats.lognorm(0.8), n=50_000, random_state=0
... )
>>> losses.shape
(50000,)
>>> bool(value_at_risk(losses, 0.99) > np.mean(losses))
True

diversification_benefit ¶

diversification_benefit(copula: Copula, margins: Margin | list[Margin], alpha: float = 0.99, weights: ArrayLike | None = None, n: int = 100000, random_state: Generator | int | None = None) -> dict[str, float]

How much expected shortfall you save because positions do not all crash together.

Compares the portfolio's expected shortfall under your copula with the "no diversification" case in which all positions move in lockstep (the comonotone copula, where each position is at the same percentile of its own loss distribution). The gap is the diversification benefit.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	The dependence model to evaluate. Its `dim` is the number of positions `d`.	required
`margins`	`scipy.stats frozen distribution or list of them, length d`	Marginal loss distributions, one per position. A single distribution is reused for every position.	required
`alpha`	`float`	Confidence level for expected shortfall, strictly between 0 and 1.	`0.99`
`weights`	`array_like of float, shape (d,)`	Position weights or exposures. Defaults to equal weights `1/d`.	`None`
`n`	`int`	Number of simulated scenarios for each of the two portfolios.	`100_000`
`random_state`	`(int, Generator or None)`	Seed or generator. An int seeds both simulations identically.	`None`

Returns:

Type	Description
`dict of str to float`	`"es"` Expected shortfall under `copula`, in loss units. `"es_comonotone"` Expected shortfall with no diversification (comonotone copula). `"benefit"` `es_comonotone - es`: the absolute reduction, in loss units. `"benefit_pct"` The reduction as a percentage of `es_comonotone` (0 to 100 in normal use; `0.0` if `es_comonotone` is zero).

Raises:

Type	Description
`ValueError`	If `alpha` is outside `(0, 1)` or `weights` has the wrong length.

Notes

Under the comonotone copula, risk measures are simply additive across positions -- that is the worst case, and the benchmark regulators use for a "no diversification" capital charge. The benefit is the shortfall against it.

Both figures are Monte Carlo estimates, so benefit carries simulation noise; for heavy-tailed margins it can come out slightly negative.

Examples:

Independence diversifies; strong dependence does not:

>>> from scipy import stats
>>> from rcopula import ClaytonCopula, IndependenceCopula
>>> from rcopula.risk import diversification_benefit
>>> margins = stats.lognorm(0.6)
>>> free = diversification_benefit(
...     IndependenceCopula(4), margins, n=40_000, random_state=0
... )
>>> tied = diversification_benefit(
...     ClaytonCopula(8.0, dim=4), margins, n=40_000, random_state=0
... )
>>> bool(free["benefit_pct"] > tied["benefit_pct"])
True

risk_contributions ¶

risk_contributions(copula: Copula, margins: Margin | list[Margin], weights: ArrayLike | None = None, alpha: float = 0.99, n: int = 100000, random_state: Generator | int | None = None) -> NDArray[np.float64]

Split the portfolio's expected shortfall into one share per position (Euler rule).

Use this to charge risk capital to desks or assets: the shares add up to the portfolio's expected shortfall, and a position's share reflects how much it loses in the portfolio's worst scenarios, not its standalone risk.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Dependence between the positions' losses. Its `dim` is the number of positions `d`.	required
`margins`	`scipy.stats frozen distribution or list of them, length d`	Marginal loss distributions, one per position. A single distribution is reused for every position.	required
`weights`	`array_like of float, shape (d,)`	Position weights or exposures. Defaults to equal weights `1/d`.	`None`
`alpha`	`float`	Confidence level for expected shortfall, strictly between 0 and 1.	`0.99`
`n`	`int`	Number of simulated scenarios.	`100_000`
`random_state`	`(int, Generator or None)`	Seed or generator. Using the same seed and `n` as :func:`simulate_losses` reproduces the same scenarios.	`None`

Returns:

Type	Description
`numpy.ndarray of float, shape (d,)`	Each position's contribution to expected shortfall, in loss units. The entries sum to the portfolio's expected shortfall.

Raises:

Type	Description
`ValueError`	If `alpha` is outside `(0, 1)`, or `weights` does not have exactly `d` entries.

Notes

Each position's contribution is its average loss conditional on the portfolio being in its own tail, weighted by exposure:

.. math:: \mathrm{ES}i = w_i\,\mathbb{E}[L_i \mid L \ge \mathrm{VaR}_\alpha].}

Contributions sum exactly to portfolio ES, which is what makes the Euler rule the standard allocation: capital charged to desks adds up to the capital held. Note that a position can carry a large share of the risk while being a small share of the portfolio, if it loads on the same tail events.

Examples:

>>> import numpy as np
>>> from scipy import stats
>>> from rcopula import ClaytonCopula
>>> from rcopula.risk import expected_shortfall, risk_contributions, simulate_losses
>>> cop, margins = ClaytonCopula(3.0, dim=3), stats.lognorm(0.7)
>>> parts = risk_contributions(cop, margins, n=60_000, random_state=0)
>>> total = expected_shortfall(
...     simulate_losses(cop, margins, n=60_000, random_state=0), 0.99
... )
>>> bool(np.isclose(parts.sum(), total))
True

covar ¶

covar(system: ArrayLike, firm: ArrayLike, alpha: float = 0.95, beta: float = 0.95, band: float = 0.05) -> float

The system's VaR on days when one firm is in distress (CoVaR).

A systemic-risk measure: how bad do losses across the whole market (or a banking system, or a portfolio) get when this particular firm is having a bad day? Compare it with the system's ordinary VaR, or use :func:delta_covar for the difference.

Parameters:

Name	Type	Description	Default
`system`	`array_like of float, shape (n,)`	Loss series for the system (index, sector, portfolio). Positive numbers are losses.	required
`firm`	`array_like of float, shape (n,)`	Loss series for the firm, paired row-by-row with `system` (same dates or same simulated scenarios).	required
`alpha`	`float`	Confidence level for the system's VaR, strictly between 0 and 1.	`0.95`
`beta`	`float`	The firm's distress level, as a quantile of its own losses: `0.95` means "the firm is at its 95th-percentile loss".	`0.95`
`band`	`float`	Half-width of the conditioning window, as a quantile fraction: rows where the firm's loss lies between its `beta - band` and `beta + band` quantiles are kept. Conditioning on an exact quantile has probability zero, so a window is unavoidable; `0.05` keeps roughly 10% of the sample.	`0.05`

Returns:

Type	Description
`float`	CoVaR, in the units of `system`.

Raises:

Type	Description
`ValueError`	If `system` and `firm` have different lengths, fewer than 10 rows fall in the conditioning window, or `alpha` is outside `(0, 1)`.

Notes

:math:\mathrm{CoVaR}_{\alpha|\beta} is the alpha-quantile of the system loss, conditional on the firm's own loss sitting at its beta-quantile (Adrian & Brunnermeier 2016).

Examples:

A firm whose losses are tied to the system raises CoVaR above the unconditional VaR; an unrelated firm does not:

>>> import numpy as np
>>> from rcopula import ClaytonCopula, IndependenceCopula
>>> from rcopula.risk import covar, value_at_risk
>>> tied = ClaytonCopula(6.0).rvs(60_000, random_state=0)
>>> free = IndependenceCopula(2).rvs(60_000, random_state=0)
>>> var = value_at_risk(tied[:, 0], 0.95)          # unconditional system VaR
>>> bool(covar(tied[:, 0], tied[:, 1], 0.95) > var + 0.02)
True
>>> bool(abs(covar(free[:, 0], free[:, 1], 0.95) - var) < 0.02)
True

delta_covar ¶

delta_covar(system: ArrayLike, firm: ArrayLike, alpha: float = 0.95, beta: float = 0.95, band: float = 0.05) -> float

How much the system's VaR rises when a firm goes from normal to distressed (ΔCoVaR).

It is :func:covar at the firm's distress level beta minus :func:covar at the firm's median (beta=0.5), so it measures the firm's marginal contribution to system risk.

Parameters:

Name	Type	Description	Default
`system`	`array_like of float, shape (n,)`	Loss series for the system. Positive numbers are losses.	required
`firm`	`array_like of float, shape (n,)`	Loss series for the firm, paired row-by-row with `system`.	required
`alpha`	`float`	Confidence level for the system's VaR, strictly between 0 and 1.	`0.95`
`beta`	`float`	The firm's distress level, as a quantile of its own losses.	`0.95`
`band`	`float`	Half-width of the conditioning window in quantile units; see :func:`covar`.	`0.05`

Returns:

Type	Description
`float`	ΔCoVaR, in the units of `system`. Near zero when the firm's losses are unrelated to the system's.

Raises:

Type	Description
`ValueError`	For the same reasons as :func:`covar`.

Notes

The difference between the system's VaR when the firm is in distress and when it is at its median. This is the quantity Adrian & Brunnermeier argue should drive systemic capital surcharges -- a firm can be small and safe on its own yet have a large ΔCoVaR, which is precisely the case unconditional measures miss.

Examples:

>>> from rcopula import ClaytonCopula, IndependenceCopula
>>> from rcopula.risk import delta_covar
>>> tied = ClaytonCopula(6.0).rvs(60_000, random_state=0)
>>> free = IndependenceCopula(2).rvs(60_000, random_state=0)
>>> bool(delta_covar(tied[:, 0], tied[:, 1]) > delta_covar(free[:, 0], free[:, 1]))
True

marginal_expected_shortfall ¶

marginal_expected_shortfall(firm: ArrayLike, system: ArrayLike, alpha: float = 0.95) -> float

A firm's average loss on the system's worst days (marginal expected shortfall, MES).

Use it to ask "how much would this firm lose in a market-wide crisis?": take the days when the system's loss is at or beyond its VaR, and average the firm's loss over those days.

Parameters:

Name	Type	Description	Default
`firm`	`array_like of float, shape (n,)`	Loss series for the firm. Positive numbers are losses.	required
`system`	`array_like of float, shape (n,)`	Loss series for the system, paired row-by-row with `firm`.	required
`alpha`	`float`	Confidence level that defines the system's "worst days", strictly between 0 and 1.	`0.95`

Returns:

Type	Description
`float`	MES, in the units of `firm`.

Raises:

Type	Description
`ValueError`	If `firm` and `system` have different lengths, or `alpha` is outside `(0, 1)`.

Notes

.. math:: \mathrm{MES}\alpha = \mathbb{E}[L \mid L_{\text{system}} \ge \mathrm{VaR}_\alpha].}

The conditioning runs the opposite way to CoVaR -- MES asks what a firm loses in a crisis, CoVaR asks what a firm's distress does to the system -- and the two can rank firms differently.

Examples:

>>> from rcopula import ClaytonCopula, IndependenceCopula
>>> from rcopula.risk import marginal_expected_shortfall as mes
>>> tied = ClaytonCopula(6.0).rvs(60_000, random_state=0)
>>> free = IndependenceCopula(2).rvs(60_000, random_state=0)
>>> bool(mes(tied[:, 0], tied[:, 1]) > mes(free[:, 0], free[:, 1]))
True

stress_scenario ¶

stress_scenario(copula: Copula, margins: Margin | list[Margin], stressed: dict[int, float], n: int = 200000, band: float = 0.02, random_state: Generator | int | None = None) -> NDArray[np.float64]

Simulate every factor given that some of them are pushed to stressed levels.

Answers "if factor 0 hits its 99th percentile, what happens to the rest?" -- with the dependence structure supplying the answer rather than an assumed correlation.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Dependence between the factors. Its `dim` is the number of factors `d`.	required
`margins`	`scipy.stats frozen distribution or list of them, length d`	Marginal distribution of each factor, used to convert the retained draws to the factors' own units. A single distribution is reused for every factor.	required
`stressed`	`dict of int to float`	Maps a factor's column index (0-based) to the quantile to condition on, e.g. `{0: 0.99}` for a 99th-percentile shock to the first factor. Must name at least one factor.	required
`n`	`int`	Number of scenarios simulated before filtering.	`200_000`
`band`	`float`	Half-width of the conditioning window in quantile units: a draw is kept when each stressed factor lies within `level +/- band`.	`0.02`
`random_state`	`(int, Generator or None)`	Seed or generator for reproducibility.	`None`

Returns:

Type	Description
`numpy.ndarray of float, shape (m, d)`	The retained scenarios, in the factors' own units (via `margins`). `m` is smaller than `n`: only the draws satisfying the conditioning survive (about `2 * band * n` for one stressed factor).

Raises:

Type	Description
`ValueError`	If `stressed` is empty, names an index outside `0..d-1`, or fewer than 10 draws survive the conditioning.

Notes

Conditioning is done by rejection on a window, not analytically. That keeps it valid for every family in the package, at the cost of discarding most draws -- roughly 2 * band of them survive per conditioned factor, so conditioning on several at once needs a large n. An exact conditional sampler via the Rosenblatt transform is the natural upgrade.

Examples:

Under Clayton, stressing one factor drags the others down with it:

>>> import numpy as np
>>> from scipy import stats
>>> from rcopula import ClaytonCopula
>>> from rcopula.risk import stress_scenario
>>> draws = stress_scenario(
...     ClaytonCopula(6.0, dim=3), stats.norm(), {0: 0.99}, n=100_000, random_state=0
... )
>>> bool(draws[:, 1].mean() > 0.5)      # unconditionally it would be 0
True

rcopula.credit ¶

Credit portfolio and structured-product utilities.

Copulas entered mainstream finance through credit: Li (2000) proposed joining individual default times with a Gaussian copula, and the whole CDO market was built on that idea. Its failure in 2007-08 is the most consequential demonstration of the point this package keeps returning to -- the Gaussian copula has no tail dependence, so it prices simultaneous defaults as essentially impossible, and senior tranches were sold on that assumption.

Everything here is deliberately model-agnostic: pass any copula and see what changes. Swapping :class:~rcopula.GaussianCopula for :class:~rcopula.StudentCopula or :class:~rcopula.ClaytonCopula at matched Kendall's tau reprices senior tranches by multiples, and reproducing that is the most useful thing this module does.

============================== ============================================== :func:default_indicators Which names default by the horizon. :func:default_times Li (2000) copula-linked default times. :func:portfolio_loss Simulated fractional portfolio loss. :func:tranche_loss Loss allocated to one tranche. :func:tranche_expected_loss Expected tranche loss, as a fraction. :func:tranche_spread Approximate fair spread. :func:nth_to_default_probability Basket default swap leg. :func:vasicek_loss_cdf Large-pool closed form -- the validation anchor. :func:implied_correlation Correlation implying a given tranche loss. ============================== ==============================================

.. warning::

These are reference implementations for analysis and teaching, not a production pricing library. Spreads use a flat-hazard, single-period approximation with a flat discount rate, no accrual on default and no counterparty adjustment.

References

Li, D. X. (2000). On default correlation: a copula function approach. Journal of Fixed Income 9(4), 43-54. Vasicek, O. (2002). The distribution of loan portfolio value. Risk 15(12), 160-162. The large-homogeneous-pool limit. Gordy, M. B. (2003). A risk-factor model foundation for ratings-based bank capital rules. Journal of Financial Intermediation 12(3), 199-232. Hull, J. and White, A. (2004). Valuation of a CDO and an nth-to-default CDS without Monte Carlo simulation. Journal of Derivatives 12(2), 8-23. MacKenzie, D. and Spears, T. (2014). "The formula that killed Wall Street": the Gaussian copula and modelling practices in investment banking. Social Studies of Science 44(3), 393-417.

default_indicators ¶

default_indicators(copula: Copula, default_prob: ArrayLike, n: int = 100000, random_state: Generator | int | None = None) -> NDArray[np.bool_]

Simulate which borrowers ("names") default before the horizon.

Each simulated scenario is one possible future: a row of True/False flags saying which names in the portfolio defaulted. Each name defaults with the probability you give it; the copula decides how often names default together. Use it as the raw building block for loss distributions, basket default swaps or your own payoff logic.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Any `rcopula` copula whose dimension `d` equals the number of names. It controls the dependence between defaults only, not their marginal probabilities.	required
`default_prob`	`float or array_like of float, shape (d,)`	Probability that each name defaults by the horizon, each in `[0, 1]` (e.g. `0.02` for 2%). A single number applies to every name.	required
`n`	`int`	Number of simulated scenarios.	`100_000`
`random_state`	`(int, Generator or None)`	Seed or generator for reproducible draws.	`None`

Returns:

Type	Description
`numpy.ndarray of bool, shape (n, d)`	`True` where the name defaulted in that scenario.

Raises:

Type	Description
`ValueError`	If any default probability lies outside `[0, 1]`, or `default_prob` has a length other than 1 or `d`.

Notes

A name defaults when its copula coordinate falls below its default probability: :math:U_i \le p_i. That is exactly the one-factor threshold model when copula is Gaussian with exchangeable correlation, but works for any family.

Examples:

The marginal default rate is whatever you asked for, regardless of copula:

>>> import numpy as np
>>> from rcopula import ClaytonCopula, GaussianCopula
>>> from rcopula.credit import default_indicators
>>> d = default_indicators(GaussianCopula(0.3, dim=50), 0.02, 40_000, 0)
>>> bool(abs(d.mean() - 0.02) < 0.002)
True

but the joint behaviour is not. Clayton produces far more all-or-nothing outcomes at the same marginal probability:

>>> gauss = default_indicators(GaussianCopula(0.3, dim=50), 0.02, 40_000, 0)
>>> clayton = default_indicators(ClaytonCopula(2.0, dim=50), 0.02, 40_000, 0)
>>> bool(clayton.sum(axis=1).max() > gauss.sum(axis=1).max())
True

default_times ¶

default_times(copula: Copula, hazard_rate: ArrayLike, n: int = 100000, random_state: Generator | int | None = None) -> NDArray[np.float64]

Simulate when each name defaults, in years, with defaults linked by a copula.

This is the Li (2000) default-time model behind classic CDO pricing. Each name keeps its own credit curve (a constant hazard rate, i.e. a constant annual default intensity); the copula decides whether early defaults tend to cluster. Use it when timing matters, e.g. first-to-default within five years: (times <= 5.0).

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Any `rcopula` copula whose dimension `d` equals the number of names.	required
`hazard_rate`	`float or array_like of float, shape (d,)`	Constant annual hazard rate of each name, strictly positive (`0.02` means roughly a 2% chance of default per year; a common rule of thumb is `spread / (1 - recovery)`). A single number applies to every name.	required
`n`	`int`	Number of simulated scenarios.	`100_000`
`random_state`	`(int, Generator or None)`	Seed or generator for reproducible draws.	`None`

Returns:

Type	Description
`numpy.ndarray of float, shape (n, d)`	Default time of each name in each scenario, in years (the same time unit as `1 / hazard_rate`). Values are positive and unbounded.

Raises:

Type	Description
`ValueError`	If any hazard rate is zero or negative, or `hazard_rate` has a length other than 1 or `d`.

Notes

With a flat hazard :math:\lambda_i, the survival function is :math:e^{-\lambda_i t}, so :math:\tau_i = -\log(1 - U_i)/\lambda_i turns a copula draw into a default time while preserving each name's own credit curve.

Examples:

Marginal default times are exponential whatever the copula:

>>> import numpy as np
>>> from scipy import stats
>>> from rcopula import ClaytonCopula
>>> from rcopula.credit import default_times
>>> t = default_times(ClaytonCopula(3.0, dim=4), 0.02, 20_000, random_state=0)
>>> bool(stats.kstest(t[:, 0], stats.expon(scale=1 / 0.02).cdf).pvalue > 0.01)
True

portfolio_loss ¶

portfolio_loss(copula: Copula, default_prob: ArrayLike, lgd: ArrayLike = 0.6, exposure: ArrayLike | None = None, n: int = 100000, random_state: Generator | int | None = None) -> NDArray[np.float64]

Simulate the credit portfolio's loss, as a fraction of total exposure.

Each scenario's loss is the sum, over the names that defaulted, of exposure times loss-given-default, divided by total exposure. A result of 0.04 means the portfolio lost 4% of its notional. Feed the output to :func:tranche_loss, :func:tranche_spread or a VaR/ES function.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Any `rcopula` copula whose dimension `d` equals the number of names.	required
`default_prob`	`float or array_like of float, shape (d,)`	Per-name default probability to the horizon, each in `[0, 1]`. A single number applies to every name.	required
`lgd`	`float or array_like of float, shape (d,)`	Loss given default, as a fraction of exposure (`0.6` means 40% recovery), each in `[0, 1]`. A single number applies to every name.	`0.6`
`exposure`	`float or array_like of float, shape (d,)`	Per-name exposure (notional), in any currency unit, each finite and `>= 0` with a positive total. Defaults to equal. Only relative sizes matter: losses are returned as a fraction of total exposure, so the result always lies in `[0, 1]`.	`None`
`n`	`int`	Number of simulated scenarios.	`100_000`
`random_state`	`(int, Generator or None)`	Seed or generator for reproducible draws.	`None`

Returns:

Type	Description
`numpy.ndarray of float, shape (n,)`	Fractional portfolio loss in each scenario.

Raises:

Type	Description
`ValueError`	If any default probability or LGD lies outside `[0, 1]`, any exposure is negative or not finite, the exposures sum to zero, or a per-name input has a length other than 1 or `d`.

Examples:

>>> import numpy as np
>>> from rcopula import GaussianCopula
>>> from rcopula.credit import portfolio_loss
>>> loss = portfolio_loss(GaussianCopula(0.2, dim=100), 0.03, 0.6, n=20_000,
...                       random_state=0)
>>> bool(0 <= loss.min() and loss.max() <= 1)
True
>>> bool(abs(loss.mean() - 0.03 * 0.6) < 0.003)      # mean loss = PD x LGD
True

tranche_loss ¶

tranche_loss(loss: ArrayLike, attachment: float, detachment: float) -> NDArray[np.float64]

Share of a tranche's notional wiped out by each portfolio loss.

A tranche is a slice of the portfolio's losses between an attachment and a detachment point. For a 3-7% tranche, a 2% portfolio loss leaves it untouched, a 5% loss wipes out half of it, and anything above 7% wipes it out completely.

Parameters:

Name	Type	Description	Default
`loss`	`float or array_like of float, shape (n,)`	Portfolio loss(es) as a fraction of total notional, typically the output of :func:`portfolio_loss`.	required
`attachment`	`float`	Lower edge of the tranche, as a fraction of portfolio notional (e.g. `0.03` for 3%).	required
`detachment`	`float`	Upper edge of the tranche, as a fraction of portfolio notional. Must satisfy `0 <= attachment < detachment <= 1`.	required

Returns:

Type	Description
numpy.ndarray of float, same shape as ``loss``	Tranche loss as a fraction of tranche notional, each in `[0, 1]`.

Raises:

Type	Description
`ValueError`	If `0 <= attachment < detachment <= 1` does not hold.

Notes

A tranche spanning :math:[a, d] absorbs :math:\min(\max(L - a, 0),\, d - a), rescaled by its width. Equity tranches take the first losses; senior tranches are untouched until the subordination below them is exhausted, which is precisely why their pricing is so sensitive to the probability of many simultaneous defaults.

Examples:

>>> import numpy as np
>>> from rcopula.credit import tranche_loss
>>> losses = np.array([0.0, 0.02, 0.05, 0.10, 0.30])
>>> tranche_loss(losses, 0.03, 0.07)     # a 3-7% mezzanine tranche
array([0. , 0. , 0.5, 1. , 1. ])

tranche_expected_loss ¶

tranche_expected_loss(loss: ArrayLike, attachment: float, detachment: float) -> float

Average share of a tranche's notional lost, across simulated scenarios.

This is the number that drives a tranche's price: the mean of :func:tranche_loss over all scenarios. 0.25 means the tranche expects to lose a quarter of its notional by the horizon.

Parameters:

Name	Type	Description	Default
`loss`	`array_like of float, shape (n,)`	Simulated portfolio losses as fractions of total notional, typically the output of :func:`portfolio_loss`.	required
`attachment`	`float`	Lower edge of the tranche, as a fraction of portfolio notional.	required
`detachment`	`float`	Upper edge of the tranche, as a fraction of portfolio notional. Must satisfy `0 <= attachment < detachment <= 1`.	required

Returns:

Type	Description
`float`	Expected tranche loss as a fraction of tranche notional, in `[0, 1]`.

Raises:

Type	Description
`ValueError`	If `0 <= attachment < detachment <= 1` does not hold.

Examples:

>>> from rcopula import GaussianCopula
>>> from rcopula.credit import portfolio_loss, tranche_expected_loss
>>> loss = portfolio_loss(GaussianCopula(0.2, dim=100), 0.05, n=40_000,
...                       random_state=0)
>>> equity = tranche_expected_loss(loss, 0.0, 0.03)
>>> senior = tranche_expected_loss(loss, 0.15, 0.30)
>>> bool(equity > senior)          # equity absorbs losses first
True

tranche_spread ¶

tranche_spread(loss: ArrayLike, attachment: float, detachment: float, maturity: float = 5.0, discount_rate: float = 0.0) -> float

Rough fair running spread on a tranche, in basis points per year.

Turns a tranche's expected loss into the annual premium a protection buyer would pay. Use it to compare how much different copulas reprice the same tranche, not to quote a trade.

Parameters:

Name	Type	Description	Default
`loss`	`array_like of float, shape (n,)`	Simulated portfolio losses as fractions of total notional, typically the output of :func:`portfolio_loss`. They should be losses to `maturity`.	required
`attachment`	`float`	Lower edge of the tranche, as a fraction of portfolio notional.	required
`detachment`	`float`	Upper edge of the tranche, as a fraction of portfolio notional. Must satisfy `0 <= attachment < detachment <= 1`.	required
`maturity`	`float`	Tranche maturity in years. Must be positive.	`5.0`
`discount_rate`	`float`	Continuously compounded annual rate (`0.03` for 3%), applied to both legs. Zero means no discounting. Because both legs are discounted on the same schedule the rate only has a second-order effect: a higher rate slightly tightens the spread (it down-weights the late, amortised premium payments).	`0.0`

Returns:

Type	Description
`float`	Approximate fair spread in basis points per year (`250.0` means 2.5% of tranche notional per year).

Raises:

Type	Description
`ValueError`	If `maturity` is not positive, or the attachment/detachment points are invalid.

Notes

Equates the protection leg to the premium leg assuming the tranche loss accrues linearly over [0, T] on a flat discount curve :math:e^{-rt}. The protection leg pays :math:\mathrm{EL}/T per year; the premium leg pays :math:s on the outstanding notional :math:1 - \mathrm{EL}\,t/T. Both are discounted with the same curve:

.. math::

s = \frac{(\mathrm{EL}/T)\,A_0}{A_0 - (\mathrm{EL}/T)\,A_1},\qquad
A_0 = \int_0^T e^{-rt}\,dt,\quad A_1 = \int_0^T t\,e^{-rt}\,dt,

which reduces to :math:s = \mathrm{EL} / (T (1 - \mathrm{EL}/2)) at :math:r = 0 (the second factor is the adjustment for notional amortising as losses accrue). Real pricing integrates over a default-time distribution with a discount curve and premium accruals; this is for comparing models, not for quoting.

Examples:

Equity trades far wider than senior:

>>> from rcopula import GaussianCopula
>>> from rcopula.credit import portfolio_loss, tranche_spread
>>> loss = portfolio_loss(GaussianCopula(0.2, dim=100), 0.05, n=40_000,
...                       random_state=0)
>>> bool(tranche_spread(loss, 0.0, 0.03) > 10 * tranche_spread(loss, 0.15, 0.30))
True

nth_to_default_probability ¶

nth_to_default_probability(copula: Copula, default_prob: ArrayLike, n_th: int = 1, n: int = 100000, random_state: Generator | int | None = None) -> float

Probability that at least n_th names in a basket default by the horizon.

This is the trigger probability of an n_th-to-default basket swap: protection pays out once the n_th default occurs. n_th=1 is first-to-default.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Any `rcopula` copula whose dimension `d` equals the basket size.	required
`default_prob`	`float or array_like of float, shape (d,)`	Per-name default probability to the horizon, each in `[0, 1]`. A single number applies to every name.	required
`n_th`	`int`	How many defaults trigger the payout, between 1 and `d`.	`1`
`n`	`int`	Number of simulated scenarios.	`100_000`
`random_state`	`(int, Generator or None)`	Seed or generator for reproducible draws.	`None`

Returns:

Type	Description
`float`	Simulated probability, in `[0, 1]`, that `n_th` or more names default.

Raises:

Type	Description
`ValueError`	If `n_th` is outside `1..d`, or any default probability lies outside `[0, 1]`.

Notes

First-to-default is worth most when defaults are independent (many chances for a first one); n-th-to-default is worth most when they are dependent (defaults arrive together). That inversion is the whole economics of correlation trading, and it is visible directly here.

Examples:

>>> from rcopula import ClaytonCopula, IndependenceCopula
>>> from rcopula.credit import nth_to_default_probability as ntd
>>> free = IndependenceCopula(10)
>>> tied = ClaytonCopula(5.0, dim=10)
>>> bool(ntd(free, 0.05, 1, 40_000, 0) > ntd(tied, 0.05, 1, 40_000, 0))
True
>>> bool(ntd(tied, 0.05, 5, 40_000, 0) > ntd(free, 0.05, 5, 40_000, 0))
True

vasicek_loss_cdf ¶

vasicek_loss_cdf(x: ArrayLike, default_prob: float, correlation: float) -> NDArray[np.float64]

Probability that a very large, uniform loan pool loses at most x (Vasicek).

A closed-form answer, with no simulation, to "how likely is my portfolio loss to stay below x?" for a pool of many small, identical loans whose defaults are linked by a one-factor Gaussian copula. It is the formula behind Basel IRB capital, and a check on simulated results.

Parameters:

Name	Type	Description	Default
`x`	`float or array_like of float`	Loss level(s) as a fraction of the pool (losses here assume 100% loss given default; to apply an LGD, evaluate at `x / lgd`).	required
`default_prob`	`float`	Default probability of each loan to the horizon, strictly in `(0, 1)`.	required
`correlation`	`float`	Asset correlation of the one-factor Gaussian model, strictly in `(0, 1)` (Basel uses roughly `0.12` to `0.24` for corporates).	required

Returns:

Type	Description
numpy.ndarray of float, same shape as ``x``	:math:`P(L \le x)`, each in `[0, 1]`. A scalar `x` gives a 0-d array; wrap it in `float()` if needed.

Raises:

Type	Description
`ValueError`	If `default_prob` or `correlation` is not strictly inside `(0, 1)`.

Notes

.. math::

P(L \le x) = \Phi\!\left(
    \frac{\sqrt{1-\rho}\,\Phi^{-1}(x) - \Phi^{-1}(p)}{\sqrt{\rho}}\right).

The closed-form limit of the one-factor Gaussian model as the number of names grows, and the foundation of the Basel IRB capital formula. It is the validation anchor for this module: a simulated Gaussian-copula portfolio with many names must converge to it.

Examples:

Simulation converges to the closed form:

>>> import numpy as np
>>> from rcopula import GaussianCopula
>>> from rcopula.credit import portfolio_loss, vasicek_loss_cdf
>>> pd_, rho = 0.05, 0.2
>>> loss = portfolio_loss(GaussianCopula(rho, dim=800), pd_, lgd=1.0,
...                       n=40_000, random_state=0)
>>> empirical = np.mean(loss <= 0.10)
>>> closed_form = float(vasicek_loss_cdf(0.10, pd_, rho))
>>> bool(abs(empirical - closed_form) < 0.02)
True

implied_correlation ¶

implied_correlation(target_expected_loss: float, default_prob: float, attachment: float, detachment: float, lgd: float = 1.0, n_names: int = 200, n: int = 40000, random_state: Generator | int | None = None) -> float

Find the Gaussian-copula correlation that reproduces a tranche's expected loss.

The credit market's "implied correlation", analogous to implied volatility: the single correlation number that makes the Gaussian copula model match an observed tranche price, here expressed as an expected loss.

Parameters:

Name	Type	Description	Default
`target_expected_loss`	`float`	Expected tranche loss to match, as a fraction of tranche notional (as returned by :func:`tranche_expected_loss`), in `[0, 1]`.	required
`default_prob`	`float`	Default probability of every name to the horizon, in `[0, 1]`.	required
`attachment`	`float`	Lower edge of the tranche, as a fraction of portfolio notional.	required
`detachment`	`float`	Upper edge of the tranche, as a fraction of portfolio notional. Must satisfy `0 <= attachment < detachment <= 1`.	required
`lgd`	`float`	Loss given default, as a fraction of exposure, for every name.	`1.0`
`n_names`	`int`	Number of equally weighted names in the simulated pool.	`200`
`n`	`int`	Number of simulated scenarios per trial correlation.	`40_000`
`random_state`	`(int, Generator or None)`	Seed or generator for reproducible results. Whatever is passed, every trial correlation reuses the same random numbers (common random numbers): a `Generator` or `None` is used once to draw a seed, which is then reused for every trial, so the Monte Carlo mismatch is a smooth function of the correlation and the root search is stable.	`None`

Returns:

Type	Description
`float`	Implied correlation, found to within about `1e-4`.

Raises:

Type	Description
`ValueError`	If the target is unattainable for any correlation in the search range `[0, 0.999]`, or the tranche bounds are invalid.

Notes

The correlation is found by Brent root search over [0, 0.999] on a Monte Carlo estimate, so it carries simulation error, and each call runs dozens of full simulations. (Exactly 1 is excluded because the equicorrelation matrix is singular there.)

Inverting the one-factor Gaussian model until it reproduces an observed tranche price, tranche by tranche on the same pool, produces the correlation skew -- different tranches implying different correlations, which is a contradiction if the model were right, and is the standard evidence that it is not.

Examples:

Round-trip: price a tranche at a known correlation, then recover it.

>>> from rcopula import GaussianCopula
>>> from rcopula.credit import (
...     implied_correlation, portfolio_loss, tranche_expected_loss)
>>> loss = portfolio_loss(GaussianCopula(0.25, dim=200), 0.05, lgd=1.0,
...                       n=40_000, random_state=1)
>>> el = tranche_expected_loss(loss, 0.03, 0.07)
>>> rho = implied_correlation(el, 0.05, 0.03, 0.07, n_names=200, random_state=1)
>>> bool(abs(rho - 0.25) < 0.06)
True

rcopula.derivatives ¶

Multi-asset option pricing under copula dependence.

A single correlation number is enough for a multi-asset payoff only when the joint law is multivariate lognormal. It is not: each underlying has its own volatility smile, and the smiles say the marginals are not lognormal. Copulas resolve the tension exactly -- take each marginal risk-neutral distribution from its own smile, and supply the dependence separately -- which is what :class:SmileMargin and :func:basket_implied_vol are for.

============================ ================================================== :func:lognormal_terminal Black-76 terminal distribution, as a margin. :class:SmileMargin Risk-neutral marginal implied by a volatility smile. :func:basket_option Option on a weighted basket. :func:rainbow_option Best-of / worst-of on several assets. :func:spread_option Option on the difference of two assets. :func:cms_spread_option Option on a spread of two CMS rates. :func:cms_convexity_adjustment The CMS convexity correction. :func:cms_margin A convexity-adjusted CMS rate, as a margin. :func:black76 Single-asset closed form. :func:margrabe Exchange option -- exact, the validation anchor. :func:kirk_spread Kirk's spread-option approximation. :func:implied_volatility Invert Black-76 for a quoted price. :func:basket_implied_vol The basket's own smile, implied by the copula. ============================ ==================================================

Prices are Monte-Carlo unless stated otherwise, and a standard_error is returned alongside so the noise is visible rather than assumed away.

.. warning::

Reference implementations for analysis, not a production pricing library. Single flat rate, no dividends beyond what the forward already embeds, no early exercise, no calibration to a term structure.

References

Margrabe, W. (1978). The value of an option to exchange one asset for another. Journal of Finance 33(1), 177-186. Kirk, E. (1995). Correlation in the energy markets. In Managing Energy Price Risk. Risk Publications. Breeden, D. T. and Litzenberger, R. H. (1978). Prices of state-contingent claims implicit in option prices. Journal of Business 51(4), 621-651. Cherubini, U., Luciano, E. and Vecchiato, W. (2004). Copula Methods in Finance. Wiley. Black, F. (1976). The pricing of commodity contracts. Journal of Financial Economics 3(1-2), 167-179. Hull, J. C. (2018). Options, Futures, and Other Derivatives, 10th ed., ch. 34. The convexity adjustment used by :func:cms_convexity_adjustment. Hagan, P. S. (2003). Convexity conundrums: pricing CMS swaps, caps, and floors. Wilmott Magazine, March, 38-44. The replication view of the same adjustment.

MonteCarloPrice ¶

Bases: NamedTuple

A simulated option price together with how noisy that estimate is.

Every Monte-Carlo pricer in this module returns one of these. Treat price +/- 2 * standard_error as a rough 95% band: if two prices differ by less than a couple of standard errors, the simulation cannot tell them apart. Increase n in the pricer to shrink the error (it falls like 1 / sqrt(n)).

Being a NamedTuple, it unpacks as price, se, n = result.

Attributes:

Name	Type	Description
`price`	`float`	Discounted present value of the option, in the same currency units as the underlying prices (times the notional, where there is one).
`standard_error`	`float`	Monte-Carlo standard error of `price`, same units. Measures simulation noise only, not model error.
`n`	`int`	Number of simulated scenarios the estimate is based on.

SmileMargin ¶

SmileMargin(strikes: ArrayLike, vols: ArrayLike, forward: float, maturity: float, rate: float = 0.0, tails: str = 'lognormal')

The distribution of an asset's price at expiry, read off its volatility smile.

Option prices across strikes tell you the market's (risk-neutral) probability of the asset finishing at each level. This class turns a quoted smile into that distribution, so a basket or rainbow option can be priced with each underlying keeping the skew its own options show, instead of assuming it is lognormal. Pass instances as margins to :func:basket_option, :func:rainbow_option, :func:basket_implied_vol and friends.

Parameters:

Name	Type	Description	Default
`strikes`	`array_like of float, shape (k,)`	Strikes at which the smile is quoted, strictly increasing, in the same currency units as `forward`. At least four are needed. Outside the quoted range the tails are extrapolated (see `tails` and Notes).	required
`vols`	`array_like of float, shape (k,)`	Black-76 implied volatilities at those strikes, as decimals (`0.25` for 25%).	required
`forward`	`float`	Forward price of the asset for delivery at `maturity`.	required
`maturity`	`float`	Time to expiry in years.	required
`rate`	`float`	Continuously compounded discount rate, as a decimal.	`0.0`
`tails`	`('lognormal', 'flat')`	How to extend the distribution beyond the quoted strikes. `"lognormal"` attaches lognormal tails using the implied vol at the nearest end strike, so quantiles can fall outside the quoted range. `"flat"` holds the CDF constant outside the range (the behaviour of earlier versions): all tail mass sits on the end strikes and :meth:`ppf` never leaves `[strikes[0], strikes[-1]]`.	`"lognormal"`

Attributes:

Name	Type	Description
`strikes`	`numpy.ndarray of float, shape (k,)`	The quoted strikes.
`vols`	`numpy.ndarray of float, shape (k,)`	The quoted implied vols.
`forward`	`float`	Forward price.
`maturity`	`float`	Time to expiry in years.
`rate`	`float`	Discount rate.

Raises:

Type	Description
`ValueError`	If `strikes` and `vols` differ in length, fewer than four strikes are given, the strikes are not strictly increasing, `forward` or `maturity` is not positive, or `tails` is unknown.

Notes

Breeden & Litzenberger (1978): the undiscounted call price determines the whole risk-neutral law, with

.. math:: F(K) = 1 + e^{rT}\,\frac{\partial C}{\partial K}.

This is what lets a basket be priced without assuming its components are lognormal. Each underlying keeps the distribution its own smile implies, and the copula supplies the dependence -- which a single correlation number cannot do once the marginals are non-lognormal.

The derivative is taken numerically on the quoted grid and the result is clipped to [0, 1] and forced to be non-decreasing. Between strikes the CDF is linearly interpolated.

Beyond the quoted range (tails="lognormal"), the left tail is a lognormal with the forward and the implied vol of the lowest strike, rescaled so the CDF is continuous there: :math:F(x) = F(K_1)\,G_L(x)/G_L(K_1) for :math:x < K_1. The right tail is the mirror image using the highest strike's vol: :math:1 - F(x) = (1 - F(K_k))\,\bar G_R(x)/\bar G_R(K_k). Holding the end-strike vol flat in the wings is the usual simple extrapolation; it keeps the CDF continuous and monotone and puts no mass on the end strikes themselves. With tails="flat" the CDF is held constant outside the range instead.

Examples:

A flat smile must reproduce the lognormal margin it came from:

>>> import numpy as np
>>> from rcopula.derivatives import SmileMargin, lognormal_terminal
>>> strikes = np.linspace(50, 200, 60)
>>> flat = SmileMargin(strikes, np.full(60, 0.25), forward=100.0, maturity=1.0)
>>> exact = lognormal_terminal(100.0, 0.25, 1.0)
>>> bool(abs(flat.cdf(110.0) - exact.cdf(110.0)) < 0.01)
True

A downward-sloping smile puts more mass in the left tail, as it should:

>>> skewed = SmileMargin(strikes, 0.25 + 0.0007 * (100 - strikes),
...                      forward=100.0, maturity=1.0)
>>> bool(skewed.cdf(70.0) > flat.cdf(70.0))
True

cdf ¶

cdf(x: ArrayLike) -> NDArray[np.float64]

Probability that the asset finishes at or below x (risk-neutral CDF).

Parameters:

Name	Type	Description	Default
`x`	`array_like of float`	Price level(s) at expiry, any shape.	required

Returns:

Type	Description
numpy.ndarray of float, same shape as ``x``	Probabilities in [0, 1].

ppf ¶

ppf(q: ArrayLike) -> NDArray[np.float64]

The price level the asset finishes below with probability q (quantile).

This is the inverse of :meth:cdf, computed by inverting the interpolated CDF. It is what the copula machinery calls to turn uniform draws into prices.

Parameters:

Name	Type	Description	Default
`q`	`array_like of float`	Probabilities in [0, 1], any shape.	required

Returns:

Type	Description
numpy.ndarray of float, same shape as ``q``	Price levels. With `tails="lognormal"` these extend below `strikes[0]` (towards 0 as `q -> 0`) and above `strikes[-1]` (to `inf` at `q = 1`); with `tails="flat"` they stay within `[strikes[0], strikes[-1]]`.

pdf ¶

pdf(x: ArrayLike) -> NDArray[np.float64]

How likely each price level at expiry is (risk-neutral density).

Equivalently, the (undiscounted) second derivative of the call price with respect to strike. Computed numerically from the CDF grid and floored at zero.

Parameters:

Name	Type	Description	Default
`x`	`array_like of float`	Price level(s) at expiry, any shape.	required

Returns:

Type	Description
numpy.ndarray of float, same shape as ``x``	Non-negative density values (probability per unit of price).

CmsLeg ¶

Bases: NamedTuple

One leg of a CMS spread option: which swap rate, and how volatile it is.

A small record passed in a list of two to :func:cms_spread_option, e.g. CmsLeg(0.045, 0.22, 10.0) for a 10-year rate at 4.5% with 22% lognormal vol.

Attributes:

Name	Type	Description
`forward`	`float`	Forward swap rate, as a decimal.
`vol`	`float`	Volatility of that rate, on the scale implied by `model` (lognormal decimal, or absolute rate units for `"normal"`).
`tenor`	`float`	Swap tenor in years -- the "constant maturity".
`frequency`	`int, default 2`	Fixed-leg payments per year.
`model`	`{"lognormal", "normal"}, default "lognormal"`	Marginal dynamics.

lognormal_terminal ¶

lognormal_terminal(forward: float, vol: float, maturity: float) -> Margin

The distribution of an asset's price at expiry, assuming a flat Black-76 vol.

Use this as a margin (one per underlying) for the multi-asset pricers in this module when you have a single at-the-money volatility per asset and no smile. If you do have a smile, use :class:SmileMargin instead.

Parameters:

Name	Type	Description	Default
`forward`	`float`	Forward price of the asset for delivery at `maturity` (currency units). The distribution's mean equals this.	required
`vol`	`float`	Annualised Black-76 (lognormal) volatility as a decimal, e.g. `0.25` for 25%. Must be positive.	required
`maturity`	`float`	Time to expiry in years. Must be positive.	required

Returns:

Type	Description
`scipy.stats frozen distribution`	A frozen `scipy.stats.lognorm` for the terminal price :math:`S_T`, with `.cdf`, `.ppf`, `.rvs`, `.mean` etc.

Raises:

Type	Description
`ValueError`	If `vol` or `maturity` is not positive.

Notes

:math:S_T = F \exp(\sigma\sqrt{T}\,Z - \sigma^2 T/2), so :math:\mathbb{E}[S_T] = F -- the martingale property that makes the resulting prices arbitrage-free.

Examples:

>>> import numpy as np
>>> from rcopula.derivatives import lognormal_terminal
>>> m = lognormal_terminal(100.0, 0.25, 1.0)
>>> bool(abs(m.mean() - 100.0) < 1e-9)
True

black76 ¶

black76(forward: float, strike: float, vol: float, maturity: float, rate: float = 0.0, kind: str = 'call') -> float

Price of a plain European call or put on one asset, using the Black-76 formula.

The single-asset building block: quote a forward, a strike and a vol, get a price. Used here to check the multi-asset pricers and to convert prices to implied vols.

Parameters:

Name	Type	Description	Default
`forward`	`float`	Forward price of the underlying for delivery at `maturity`.	required
`strike`	`float`	Option strike, same units as `forward`.	required
`vol`	`float`	Annualised Black-76 (lognormal) volatility as a decimal. If zero or negative, the discounted intrinsic value is returned.	required
`maturity`	`float`	Time to expiry in years. If zero or negative, the discounted intrinsic value is returned.	required
`rate`	`float`	Continuously compounded discount rate, as a decimal.	`0.0`
`kind`	`('call', 'put')`	Option type.	`"call"`

Returns:

Type	Description
`float`	Discounted option price, in the units of `forward`.

Raises:

Type	Description
`ValueError`	If `kind` is not `"call"` or `"put"`.

Notes

:math:C = e^{-rT}[F N(d_1) - K N(d_2)] with :math:d_{1,2} = [\ln(F/K) \pm \sigma^2 T/2] / (\sigma\sqrt{T}) (Black, 1976).

Examples:

>>> from rcopula.derivatives import black76
>>> float(round(black76(100.0, 100.0, 0.2, 1.0), 6))
7.965567

Put-call parity holds:

>>> c = black76(100.0, 90.0, 0.2, 1.0)
>>> p = black76(100.0, 90.0, 0.2, 1.0, kind="put")
>>> bool(abs((c - p) - (100.0 - 90.0)) < 1e-10)
True

margrabe ¶

margrabe(forward1: float, forward2: float, vol1: float, vol2: float, correlation: float, maturity: float, rate: float = 0.0) -> float

Exact price of the right to swap asset 2 for asset 1 at expiry (Margrabe formula).

The payoff is whatever asset 1 is worth above asset 2, or nothing: :math:\max(S_1 - S_2, 0), a spread option with zero strike. It is exact when both prices are lognormal with a constant correlation, which makes it the benchmark for checking :func:spread_option.

Parameters:

Name	Type	Description	Default
`forward1`	`float`	Forward price of asset 1 (the one you receive).	required
`forward2`	`float`	Forward price of asset 2 (the one you give up), same units.	required
`vol1`	`float`	Annualised lognormal volatility of asset 1, as a decimal.	required
`vol2`	`float`	Annualised lognormal volatility of asset 2, as a decimal.	required
`correlation`	`float`	Correlation of the two assets' log-returns, in [-1, 1].	required
`maturity`	`float`	Time to expiry in years.	required
`rate`	`float`	Continuously compounded discount rate, as a decimal.	`0.0`

Returns:

Type	Description
`float`	Discounted option price, in the units of the forwards. If the spread volatility is zero, the discounted intrinsic value.

Notes

Margrabe (1978) showed this is Black-Scholes with the spread volatility

.. math:: \sigma^2 = \sigma_1^2 + \sigma_2^2 - 2\rho\sigma_1\sigma_2,

and no strike. Because it is exact under jointly lognormal dynamics -- which is precisely a Gaussian copula with lognormal margins -- it is the validation anchor for :func:spread_option.

Examples:

>>> from rcopula.derivatives import margrabe
>>> float(round(margrabe(100.0, 95.0, 0.2, 0.3, 0.5, 1.0), 6))
12.952273

Perfect correlation with equal vols leaves only the intrinsic difference:

>>> float(round(margrabe(100.0, 95.0, 0.25, 0.25, 1.0, 1.0), 10))
5.0

kirk_spread ¶

kirk_spread(forward1: float, forward2: float, strike: float, vol1: float, vol2: float, correlation: float, maturity: float, rate: float = 0.0) -> float

Fast approximate price of a spread option with a strike (Kirk's formula).

Payoff :math:\max(S_1 - S_2 - K, 0): pays when asset 1 beats asset 2 by more than strike. The standard closed-form quick price for crack, spark and calendar spreads; compare it with :func:spread_option to see what a non-Gaussian dependence changes.

Parameters:

Name	Type	Description	Default
`forward1`	`float`	Forward price of the long asset.	required
`forward2`	`float`	Forward price of the short asset, same units.	required
`strike`	`float`	Spread strike, same units. `0` reproduces :func:`margrabe`.	required
`vol1`	`float`	Annualised lognormal volatility of asset 1, as a decimal.	required
`vol2`	`float`	Annualised lognormal volatility of asset 2, as a decimal.	required
`correlation`	`float`	Correlation of the two assets' log-returns, in [-1, 1].	required
`maturity`	`float`	Time to expiry in years.	required
`rate`	`float`	Continuously compounded discount rate, as a decimal.	`0.0`

Returns:

Type	Description
`float`	Approximate discounted option price, in the units of the forwards.

Notes

Kirk (1995) treats :math:S_2 + K as a single lognormal asset, which is exact at :math:K = 0 (where it reduces to :func:margrabe) and stays accurate for moderate strikes. Widely used in energy markets, where spread options are the standard product.

Examples:

At zero strike it coincides with Margrabe:

>>> from rcopula.derivatives import kirk_spread, margrabe
>>> a = kirk_spread(100.0, 95.0, 0.0, 0.2, 0.3, 0.5, 1.0)
>>> b = margrabe(100.0, 95.0, 0.2, 0.3, 0.5, 1.0)
>>> bool(abs(a - b) < 1e-10)
True

implied_volatility ¶

implied_volatility(price: float, forward: float, strike: float, maturity: float, rate: float = 0.0, kind: str = 'call') -> float

The volatility that makes Black-76 reproduce a given option price.

The usual way to quote an option price as a vol. Here it is also how :func:basket_implied_vol turns simulated basket prices into a smile.

Parameters:

Name	Type	Description	Default
`price`	`float`	Discounted option price to match, in the units of `forward`.	required
`forward`	`float`	Forward price of the underlying.	required
`strike`	`float`	Option strike, same units.	required
`maturity`	`float`	Time to expiry in years.	required
`rate`	`float`	Continuously compounded discount rate, as a decimal.	`0.0`
`kind`	`('call', 'put')`	Option type of the quoted price.	`"call"`

Returns:

Type	Description
`float`	Annualised implied volatility as a decimal. Returns `0.0` when the price is at or below discounted intrinsic value.

Raises:

Type	Description
`ValueError`	If `kind` is not `"call"` or `"put"`, or no volatility in (1e-8, 10] reproduces `price` -- typically because the price exceeds the no-arbitrage upper bound.

Notes

Solved with Brent's method on black76(vol) - price.

Examples:

>>> from rcopula.derivatives import black76, implied_volatility
>>> p = black76(100.0, 110.0, 0.27, 1.5)
>>> float(round(implied_volatility(p, 100.0, 110.0, 1.5), 10))
0.27

basket_option ¶

basket_option(copula: Copula, margins: Margin | list[Margin], strike: float, maturity: float, weights: ArrayLike | None = None, rate: float = 0.0, kind: str = 'call', n: int = 200000, random_state: Generator | int | None = None) -> MonteCarloPrice

Monte-Carlo price of a call or put on a weighted basket of assets.

Payoff :math:\max(\sum_i w_i S_i - K, 0) for a call. Each asset's price at expiry comes from its own margin; the copula decides how they move together. Swap the copula (e.g. Gaussian for Student-t or Gumbel) to see how much of the price is the dependence assumption.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Dependence between the `d` underlyings, e.g. `GaussianCopula(0.4, dim=3)`.	required
`margins`	`scipy.stats frozen distribution or list of them, length d`	Terminal price distribution of each underlying, e.g. from :func:`lognormal_terminal` or :class:`SmileMargin`. A single margin is reused for every asset.	required
`strike`	`float`	Basket strike, in the units of the weighted basket value.	required
`maturity`	`float`	Time to expiry in years (used only for discounting).	required
`weights`	`array_like of float, shape (d,)`	Basket weights (number of units of each asset). Default is equal weights `1 / d`.	`None`
`rate`	`float`	Continuously compounded discount rate, as a decimal.	`0.0`
`kind`	`('call', 'put')`	Option type.	`"call"`
`n`	`int`	Number of simulated scenarios.	`200_000`
`random_state`	`(int, Generator or None)`	Seed or generator, for reproducible prices.	`None`

Returns:

Type	Description
`MonteCarloPrice`	`(price, standard_error, n)`: discounted price, its simulation standard error, and the scenario count.

Raises:

Type	Description
`ValueError`	If `kind` is not `"call"` or `"put"`, or `weights` does not have one entry per copula dimension.

Notes

The basket is where the dependence model earns its keep. A basket is less volatile than its components, by an amount that depends entirely on how they co-move -- and on how they co-move in the tail, which correlation alone does not capture.

Examples:

>>> from rcopula import GaussianCopula
>>> from rcopula.derivatives import basket_option, lognormal_terminal
>>> margins = [lognormal_terminal(100.0, 0.25, 1.0)] * 3
>>> price = basket_option(GaussianCopula(0.4, dim=3), margins, 100.0, 1.0,
...                       n=80_000, random_state=0)
>>> bool(price.price > 0 and price.standard_error < 0.2)
True

Higher dependence means less diversification and a more valuable option:

>>> low = basket_option(GaussianCopula(0.0, dim=3), margins, 100.0, 1.0,
...                     n=80_000, random_state=0).price
>>> high = basket_option(GaussianCopula(0.9, dim=3), margins, 100.0, 1.0,
...                      n=80_000, random_state=0).price
>>> bool(high > low)
True

rainbow_option ¶

rainbow_option(copula: Copula, margins: Margin | list[Margin], strike: float, maturity: float, rate: float = 0.0, kind: str = 'call', on: str = 'best', n: int = 200000, random_state: Generator | int | None = None) -> MonteCarloPrice

Monte-Carlo price of an option on the best (or worst) performer of several assets.

Payoff :math:\max(\max_i S_i - K, 0) for on="best", or with :math:\min_i for on="worst" (calls; puts flip the sign). Worst-of structures are common in structured notes, so this is where an over-optimistic dependence assumption is most expensive.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Dependence between the `d` underlyings.	required
`margins`	`scipy.stats frozen distribution or list of them, length d`	Terminal price distribution of each underlying. A single margin is reused for every asset. For best/worst comparisons to make sense the assets should be on a common scale (e.g. all with the same forward, or normalised to performance).	required
`strike`	`float`	Strike applied to the best or worst terminal price.	required
`maturity`	`float`	Time to expiry in years (used only for discounting).	required
`rate`	`float`	Continuously compounded discount rate, as a decimal.	`0.0`
`kind`	`('call', 'put')`	Option type.	`"call"`
`on`	`('best', 'worst')`	Whether the payoff references the highest or lowest terminal price.	`"best"`
`n`	`int`	Number of simulated scenarios.	`200_000`
`random_state`	`(int, Generator or None)`	Seed or generator, for reproducible prices.	`None`

Returns:

Type	Description
`MonteCarloPrice`	`(price, standard_error, n)`.

Raises:

Type	Description
`ValueError`	If `on` is not `"best"`/`"worst"` or `kind` is not `"call"`/`"put"`.

Notes

These are the payoffs most sensitive to dependence, and in opposite directions: a best-of is worth most when the assets are independent (many chances for one to finish high), a worst-of when they move together (less chance that any one drags the minimum down).

Examples:

>>> from rcopula import GaussianCopula
>>> from rcopula.derivatives import lognormal_terminal, rainbow_option
>>> margins = [lognormal_terminal(100.0, 0.3, 1.0)] * 3
>>> free, tied = GaussianCopula(0.0, dim=3), GaussianCopula(0.9, dim=3)
>>> best_free = rainbow_option(free, margins, 100.0, 1.0, n=80_000,
...                            random_state=0, on="best").price
>>> best_tied = rainbow_option(tied, margins, 100.0, 1.0, n=80_000,
...                            random_state=0, on="best").price
>>> bool(best_free > best_tied)
True

and the worst-of ordering is reversed:

>>> worst_free = rainbow_option(free, margins, 100.0, 1.0, n=80_000,
...                             random_state=0, on="worst").price
>>> worst_tied = rainbow_option(tied, margins, 100.0, 1.0, n=80_000,
...                             random_state=0, on="worst").price
>>> bool(worst_tied > worst_free)
True

spread_option ¶

spread_option(copula: Copula, margins: list[Margin], strike: float, maturity: float, rate: float = 0.0, n: int = 200000, random_state: Generator | int | None = None, kind: str = 'call') -> MonteCarloPrice

Monte-Carlo price of a call or put on the spread between two assets.

Call payoff :math:\max(S_1 - S_2 - K, 0): pays when asset 1 beats asset 2 by more than the strike (e.g. a crack or spark spread). Put payoff :math:\max(K - (S_1 - S_2), 0). Unlike :func:kirk_spread, any margins and any copula can be used.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Bivariate (`dim=2`) dependence between the two assets.	required
`margins`	`list of scipy.stats frozen distributions, length 2`	Terminal price distributions of asset 1 (long) and asset 2 (short).	required
`strike`	`float`	Spread strike, in price units.	required
`maturity`	`float`	Time to expiry in years (used only for discounting).	required
`rate`	`float`	Continuously compounded discount rate, as a decimal.	`0.0`
`n`	`int`	Number of simulated scenarios.	`200_000`
`random_state`	`(int, Generator or None)`	Seed or generator, for reproducible prices.	`None`
`kind`	`('call', 'put')`	Option type.	`"call"`

Returns:

Type	Description
`MonteCarloPrice`	`(price, standard_error, n)`.

Raises:

Type	Description
`ValueError`	If the copula is not bivariate, or `kind` is not `"call"` or `"put"`.

Notes

Put-call parity holds scenario by scenario: :math:C - P = e^{-rT}(\mathbb{E}[S_1 - S_2] - K).

At :math:K = 0 with a Gaussian copula and lognormal margins this is the Margrabe exchange option, which has an exact price -- so :func:margrabe is the check that this simulation is right.

Examples:

>>> import numpy as np
>>> from rcopula import GaussianCopula
>>> from rcopula.derivatives import lognormal_terminal, margrabe, spread_option
>>> margins = [lognormal_terminal(100.0, 0.2, 1.0), lognormal_terminal(95.0, 0.3, 1.0)]
>>> mc = spread_option(GaussianCopula(0.5), margins, 0.0, 1.0,
...                    n=400_000, random_state=0)
>>> exact = margrabe(100.0, 95.0, 0.2, 0.3, 0.5, 1.0)
>>> bool(abs(mc.price - exact) < 4 * mc.standard_error)
True

cms_convexity_adjustment ¶

cms_convexity_adjustment(forward: float, vol: float, maturity: float, tenor: float, frequency: int = 2, model: str = 'lognormal') -> float

How much to add to a forward swap rate to get the expected CMS rate actually paid.

A constant-maturity swap (CMS) pays a swap rate on a single date, not over the swap's own schedule, and that timing mismatch makes the expected paid rate higher than the forward swap rate. This returns that gap (the convexity adjustment), which grows with vol, time to fixing and swap tenor.

Parameters:

Name	Type	Description	Default
`forward`	`float`	Forward swap rate, as a decimal (`0.05` for 5%).	required
`vol`	`float`	Volatility of that rate: lognormal (a decimal, e.g. `0.20`) by default, or normal/absolute (in rate units, e.g. `0.01` for 100 bp) if `model="normal"`. Must be non-negative.	required
`maturity`	`float`	Time to the fixing, in years. Must be non-negative.	required
`tenor`	`float`	Tenor of the underlying swap, in years.	required
`frequency`	`int`	Fixed-leg payments per year.	`2`
`model`	`('lognormal', 'normal')`	Dynamics assumed for the swap rate. Normal (Bachelier) is the market standard for rates near or below zero, where a lognormal vol is meaningless.	`"lognormal"`

Returns:

Type	Description
`float`	The additive adjustment, in the same units as `forward` (multiply by 1e4 for basis points). `0.0` if `vol` or `maturity` is zero.

Raises:

Type	Description
`ValueError`	If `model` is unknown, `vol` or `maturity` is negative, `tenor` is shorter than one coupon period, or the rate is so negative that the bond price diverges.

Notes

The forward swap rate is a martingale under the annuity measure, not the payment measure, so the expected rate that actually gets paid is higher than the forward. The correction follows from the fact that the bond's forward price is the martingale: expanding :math:G(y_T) to second order in :math:\mathbb{E}[G(y_T)] = G(y_0) gives

.. math:: \mathbb{E}^T[y_T] - y_0 \;\approx\; -\tfrac{1}{2}\,\mathrm{Var}[y_T]\,\frac{G''(y_0)}{G'(y_0)},

with :math:\mathrm{Var}[y_T] = y_0^2\sigma^2 T for a lognormal rate and :math:\sigma^2 T for a normal one. Since :math:G' < 0 and :math:G'' > 0, the adjustment is positive.

Ignoring it is not a rounding error: at a 5% rate, 20% vol and 5 years to expiry on a 10-year swap it is worth several basis points, and a CMS spread option is a difference of two such rates, so the two adjustments do not cancel -- the longer tenor carries the larger one, which is the whole directional content of the trade.

This is the second-order approximation. Hagan's replication approach prices the same quantity as a strip of swaptions and so captures the smile; that is more accurate but needs a full swaption surface. What this function gives is the standard first-cut number, and :func:cms_margin builds on it.

Examples:

Around 24 basis points on a 10-year rate fixing in 5 years -- far too large to leave out of a spread that is itself only ~150 bp wide:

>>> from rcopula.derivatives import cms_convexity_adjustment
>>> ca = cms_convexity_adjustment(0.05, 0.20, 5.0, 10.0)
>>> float(round(ca * 1e4, 2))
23.62

It is positive, and grows with volatility, expiry and tenor:

>>> longer = cms_convexity_adjustment(0.05, 0.20, 5.0, 30.0)
>>> bool(longer > ca > 0)
True

With no volatility there is nothing to adjust:

>>> cms_convexity_adjustment(0.05, 0.0, 5.0, 10.0)
0.0

cms_margin ¶

cms_margin(forward: float, vol: float, maturity: float, tenor: float, frequency: int = 2, model: str = 'lognormal') -> Margin

The distribution of a CMS rate at its fixing, centred on the convexity-adjusted rate.

Use it as one margin of a copula model of two or more swap rates; this is what :func:cms_spread_option does for each leg.

Parameters:

Name	Type	Description	Default
`forward`	`float`	Forward swap rate, as a decimal.	required
`vol`	`float`	Volatility of the rate: lognormal (decimal) for `model="lognormal"`, absolute (rate units) for `model="normal"`. Must be positive for either model.	required
`maturity`	`float`	Time to the fixing, in years. Must be positive.	required
`tenor`	`float`	Tenor of the underlying swap, in years.	required
`frequency`	`int`	Fixed-leg payments per year.	`2`
`model`	`('lognormal', 'normal')`	Shape of the distribution.	`"lognormal"`

Returns:

Type	Description
`scipy.stats frozen distribution`	`scipy.stats.lognorm` (via :func:`lognormal_terminal`) or `scipy.stats.norm`, with mean `forward + cms_convexity_adjustment(...)`.

Raises:

Type	Description
`ValueError`	As :func:`cms_convexity_adjustment`; and, for either model, if `vol` or `maturity` is not positive. A zero vol or zero time to fixing would make the rate a known constant (a point mass), which is not a continuous margin a copula can use.

Notes

The mean is the forward swap rate plus :func:cms_convexity_adjustment; the shape is lognormal or normal according to model. Plugging these into a copula is how a CMS spread option gets priced with each leg keeping its own dynamics.

Examples:

The mean is the adjusted rate, not the forward:

>>> from rcopula.derivatives import cms_convexity_adjustment, cms_margin
>>> m = cms_margin(0.05, 0.20, 5.0, 10.0)
>>> ca = cms_convexity_adjustment(0.05, 0.20, 5.0, 10.0)
>>> bool(abs(m.mean() - (0.05 + ca)) < 1e-12)
True

cms_spread_option ¶

cms_spread_option(copula: Copula, legs: list[CmsLeg], strike: float, maturity: float, rate: float = 0.0, notional: float = 1.0, kind: str = 'call', n: int = 200000, random_state: Generator | int | None = None) -> MonteCarloPrice

Monte-Carlo price of an option on the gap between two swap rates (a curve trade).

Payoff :math:N \max(y_1 - y_2 - K, 0) on the two convexity-adjusted swap rates. The classic trade is 10-year minus 2-year: a bet on the slope of the curve rather than on its level.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Bivariate (`dim=2`) dependence between the two rates.	required
`legs`	`list of CmsLeg, length 2`	The two legs, in the order they appear in the spread (`legs[0] - legs[1]`).	required
`strike`	`float`	Spread strike, as a decimal (`0.015` for 150 bp).	required
`maturity`	`float`	Fixing date, in years; also used for discounting.	required
`rate`	`float`	Continuously compounded discount rate, as a decimal.	`0.0`
`notional`	`float`	Notional the rate spread is paid on.	`1.0`
`kind`	`('call', 'put')`	A call pays when the curve steepens beyond the strike; a put when it flattens.	`"call"`
`n`	`int`	Number of simulated scenarios.	`200_000`
`random_state`	`(int, Generator or None)`	Seed or generator, for reproducible prices.	`None`

Returns:

Type	Description
`MonteCarloPrice`	`(price, standard_error, n)`, in units of `notional`.

Raises:

Type	Description
`ValueError`	If the copula is not bivariate, `legs` does not have exactly two entries, or `kind` is not `"call"`/`"put"`.

Notes

This is a case where the copula is doing real work. The two rates are strongly dependent but not lognormally so; the spread is a small difference of two large numbers, so its distribution is extremely sensitive to how the joint tails behave; and the market quotes each leg's smile separately. A single correlation applied to two lognormals gets all three wrong at once. Choosing the dependence structure explicitly at least makes the assumption visible -- and, in the Gaussian-copula/lognormal case, reduces to the Margrabe/Kirk answer that desks already use.

Examples:

A 10y-2y steepener struck at the current spread:

>>> import rcopula as rc
>>> from rcopula.derivatives import CmsLeg, cms_spread_option
>>> legs = [CmsLeg(0.045, 0.22, 10.0), CmsLeg(0.030, 0.28, 2.0)]
>>> mc = cms_spread_option(rc.GaussianCopula(0.85), legs, 0.015, 5.0,
...                        n=200_000, random_state=0)
>>> bool(mc.price > 0.0)
True

Higher correlation squeezes the spread's distribution and so cheapens the option -- the correlation exposure that makes these trades interesting:

>>> tight = cms_spread_option(rc.GaussianCopula(0.97), legs, 0.015, 5.0,
...                           n=200_000, random_state=0)
>>> bool(tight.price < mc.price)
True

Lower-tail dependence, which pins the two rates together when both fall, prices differently again even at the same rank correlation:

>>> clayton = cms_spread_option(rc.ClaytonCopula.from_tau(rc.GaussianCopula(0.85).tau()),
...                             legs, 0.015, 5.0, n=200_000, random_state=0)
>>> bool(abs(clayton.price - mc.price) > clayton.standard_error)
True

basket_implied_vol ¶

basket_implied_vol(copula: Copula, margins: Margin | list[Margin], strikes: ArrayLike, maturity: float, weights: ArrayLike | None = None, rate: float = 0.0, n: int = 200000, random_state: Generator | int | None = None) -> tuple[NDArray[np.float64], NDArray[np.float64]]

The implied-volatility smile of a basket, given its components and a copula.

Prices a basket call at each strike by simulation, then converts each price to a Black-76 implied vol. Use it to see what smile your component smiles plus a dependence assumption imply for the basket, e.g. to compare with quoted index options.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Dependence between the `d` components.	required
`margins`	`scipy.stats frozen distribution or list of them, length d`	Terminal price distribution of each component (e.g. :class:`SmileMargin`). A single margin is reused for every asset.	required
`strikes`	`array_like of float, shape (k,)`	Basket strikes at which to compute the implied vol.	required
`maturity`	`float`	Time to expiry in years.	required
`weights`	`array_like of float, shape (d,)`	Basket weights. Default is equal weights `1 / d`.	`None`
`rate`	`float`	Continuously compounded discount rate, as a decimal.	`0.0`
`n`	`int`	Number of simulated scenarios (shared across all strikes).	`200_000`
`random_state`	`(int, Generator or None)`	Seed or generator, for reproducible results.	`None`

Returns:

Name	Type	Description
`strikes`	`numpy.ndarray of float, shape (k,)`	The input strikes.
`vols`	`numpy.ndarray of float, shape (k,)`	Annualised implied volatility at each strike, as decimals.

Raises:

Type	Description
`ValueError`	If `weights` does not have length `d`, or a simulated price cannot be inverted to a volatility (see :func:`implied_volatility`).

Notes

This is the calculation a single correlation number cannot do: give each component the marginal its own smile implies (see :class:SmileMargin), choose a dependence structure, and the basket smile falls out.

The forward used in the inversion is the simulated basket mean, not the analytic one, so Monte-Carlo noise in the mean does not show up as a spurious skew. All strikes share the same simulated scenarios.

Examples:

Lognormal components under a Gaussian copula give a nearly flat basket smile -- the basket is then close to lognormal itself:

>>> import numpy as np
>>> from rcopula import GaussianCopula
>>> from rcopula.derivatives import basket_implied_vol, lognormal_terminal
>>> margins = [lognormal_terminal(100.0, 0.25, 1.0)] * 3
>>> k, v = basket_implied_vol(GaussianCopula(0.5, dim=3), margins,
...                           [90, 100, 110], 1.0, n=200_000, random_state=0)
>>> bool(v.std() < 0.02)
True

A tail-dependent copula bends it, because the basket is no longer lognormal:

>>> from rcopula import GumbelCopula
>>> _, vg = basket_implied_vol(GumbelCopula.from_tau(1 / 3, dim=3), margins,
...                            [90, 100, 110], 1.0, n=200_000, random_state=0)
>>> bool(vg.std() > v.std())
True

rcopula.portfolio ¶

Portfolio construction and copula pairs trading.

Two applications that both turn on the same idea: the conditional distribution is the signal.

Pairs trading. The classical approach standardises a price spread and trades its deviations, which implicitly assumes the two assets are jointly normal with a stable linear relationship. A copula replaces that with :math:h(u_1 \mid u_2) = P(U_1 \le u_1 \mid U_2 = u_2): the probability of seeing asset 1 this low given where asset 2 actually is. When that probability is 2%, asset 1 is cheap relative to its partner in a sense that survives non-normal margins, asymmetric dependence and tail clustering -- none of which a z-scored spread can express. The quantity is uniform under the fitted model, so the same thresholds mean the same thing across every pair.

Mean-CVaR optimisation. Mean-variance treats upside and downside alike and is blind to tail dependence. Minimising conditional value at risk instead targets the loss that actually matters, and Rockafellar & Uryasev showed the problem is a linear program once returns are represented by scenarios -- which is exactly what a copula model produces.

============================== ============================================== :func:mispricing_index Conditional probabilities for a pair. :func:pairs_signal Entry and exit signals from those. :func:backtest_pairs Walk-forward backtest of the strategy. :func:simulate_returns Scenario returns from a copula model. :func:mean_cvar_weights CVaR-optimal weights, by linear programming. :func:efficient_frontier The mean-CVaR frontier. :func:min_variance_weights Markowitz benchmark, for comparison. ============================== ==============================================

References

Rockafellar, R. T. and Uryasev, S. (2000). Optimization of conditional value-at-risk. Journal of Risk 2(3), 21-42. The linear-programming formulation used by :func:mean_cvar_weights. Liew, R. Q. and Wu, Y. (2013). Pairs trading: a copula approach. Journal of Derivatives & Hedge Funds 19(1), 12-30. Stander, Y., Marais, D. and Botha, I. (2013). Trading strategies with copulas. Journal of Economic and Financial Sciences 6(1), 83-107. Xie, W., Liew, R. Q., Wu, Y. and Zou, X. (2016). Pairs trading with copulas. Journal of Trading 11(3), 41-52.

BacktestResult ¶

Bases: NamedTuple

Headline numbers and daily detail from :func:backtest_pairs.

A named tuple, so fields can be read by name (result.hit_rate) or unpacked in order. All figures are frictionless (no costs or slippage).

Attributes:

Name	Type	Description
`total_return`	`float`	Compounded return of the strategy over the whole sample, as a fraction (0.05 means +5%).
`annualised_sharpe`	`float`	Mean over standard deviation (`ddof=1`) of the per-period strategy return over every live period -- all periods after the first `train + 1` warm-up rows, flat periods counting as a zero return -- scaled by `sqrt(periods_per_year)`. `0.0` if the returns have zero variance.
`n_trades`	`int`	Number of trades opened: each move from flat into a position, and each direct flip from long to short or back, counts once. Closing a position does not count, so this is the number of round trips.
`hit_rate`	`float`	Fraction of in-position periods with a positive strategy return, between 0 and 1. `0.0` if no position was ever held.
`returns`	`numpy.ndarray of float, shape (n,)`	Per-period strategy return, `position * (r1 - r2)`; zero when flat.
`positions`	`numpy.ndarray of int, shape (n,)`	Position held in each period: `+1`, `-1` or `0`. The first `train + 1` entries are always `0`.

EfficientFrontier ¶

Bases: NamedTuple

Points on the mean-CVaR frontier, as returned by :func:efficient_frontier.

A named tuple, so it unpacks like the plain 3-tuple returned by earlier versions: mu, cvar, w = efficient_frontier(r).

Attributes:

Name	Type	Description
`achieved_returns`	`numpy.ndarray of float, shape (m,)`	Mean portfolio return achieved by the optimal weights at each point, `weights @ scenarios.mean(axis=0)`. It is at least the target return that was asked for, and can exceed it where the return constraint does not bind.
`cvars`	`numpy.ndarray of float, shape (m,)`	CVaR (expected shortfall of the portfolio loss) at each point, as a positive loss.
`weights`	`numpy.ndarray of float, shape (m, d)`	Portfolio weights at each point; each row sums to one.

mispricing_index ¶

mispricing_index(copula: Copula, u: ArrayLike) -> tuple[NDArray, NDArray]

How cheap or rich each asset of a pair looks, given where its partner is.

For each day, returns two probabilities between 0 and 1. h1 is the chance that asset 1 would be this low given today's level of asset 2; h2 is the same for asset 2 given asset 1. A small h1 (say 0.02) means asset 1 is unusually cheap relative to asset 2; a value near 0.5 means nothing unusual. These are the inputs to :func:pairs_signal.

Technically these are the copula conditional distribution functions (h-functions) :math:h_1 = P(U_1 \le u_1 \mid U_2 = u_2) and :math:h_2 = P(U_2 \le u_2 \mid U_1 = u_1).

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	A bivariate copula with its parameters already set, usually fitted to a training window of the pair (see :func:`rcopula.fit`).	required
`u`	`array_like of float, shape (n, 2) or (2,)`	Pseudo-observations: each asset's return converted to a rank in (0, 1), for example with :func:`rcopula.pseudo_obs`. Column 0 is asset 1, column 1 is asset 2. A single row of length 2 is accepted.	required

Returns:

Name	Type	Description
`h1`	`numpy.ndarray of float, shape (n,)`	:math:`P(U_1 \le u_1 \mid U_2 = u_2)` for each row: small means asset 1 is low relative to asset 2.
`h2`	`numpy.ndarray of float, shape (n,)`	:math:`P(U_2 \le u_2 \mid U_1 = u_1)` for each row: small means asset 2 is low relative to asset 1.

Raises:

Type	Description
`ValueError`	If `u` does not have exactly two columns.

Notes

Both are uniform under the fitted copula, so a value of 0.02 always means "only a 2% chance of being this low" whatever the pair, the margins or the dependence shape. That comparability is what a z-scored price spread cannot give you.

Examples:

>>> import numpy as np
>>> from rcopula import ClaytonCopula
>>> from rcopula.portfolio import mispricing_index
>>> cop = ClaytonCopula(3.0)
>>> h1, h2 = mispricing_index(cop, cop.rvs(5000, random_state=0))
>>> bool(abs(h1.mean() - 0.5) < 0.02 and abs(h2.mean() - 0.5) < 0.02)
True

Asset 1 unusually low given a high asset 2 shows up as a small h1:

>>> h1, h2 = mispricing_index(cop, [[0.05, 0.90]])
>>> bool(h1[0] < 0.05 and h2[0] > 0.95)
True

pairs_signal ¶

pairs_signal(copula: Copula, u: ArrayLike, entry: float = 0.05, exit_band: float = 0.5) -> NDArray[np.int_]

Long, short or flat for a pair, based on how mispriced one leg looks.

Returns +1 to go long the spread (long asset 1, short asset 2), -1 for the reverse, and 0 to stand aside, one value per row of u.

The rows of u are read in time order and the signal has memory:

Entry. From flat, a long position (+1) opens when both legs agree that asset 1 is cheap relative to asset 2 -- h1 at or below entry and h2 at or above 1 - entry (see :func:mispricing_index); a short position (-1) opens on the mirror condition. Requiring both is what distinguishes relative mispricing from a common move: if the whole market falls, both conditionals stay near 0.5 and no signal fires.
Hold. An open position is kept on later rows until the mispricing has reverted: a long closes once h1 >= 0.5 - exit_band or h2 <= 0.5 + exit_band (either leg is back inside the band [0.5 - exit_band, 0.5 + exit_band] or beyond it); a short closes on the mirror condition.
Flip. An opposite entry signal while a position is open reverses it directly.

With the default exit_band=0.5 the band is all of [0, 1], so every position closes on the next row unless the entry condition still holds: a memoryless, one-row signal.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	A bivariate copula with its parameters already set (usually fitted to a training window of the pair).	required
`u`	`array_like of float, shape (n, 2) or (2,)`	Pseudo-observations (ranks in (0, 1)) of the two assets; column 0 is asset 1, column 1 is asset 2.	required
`entry`	`float`	Probability threshold for opening, strictly between 0 and 0.5. Smaller means rarer, higher-conviction trades.	`0.05`
`exit_band`	`float`	Half-width of the exit band around 0.5, in `[0, 0.5]`. A position is held until the conditional probability that triggered it crosses back inside `[0.5 - exit_band, 0.5 + exit_band]`. `0.0` holds until the conditional crosses 0.5 itself (full reversion); `0.5` (the default) closes on the next row unless the entry rule fires again. Values between `0.5 - entry` and `0.5` behave like `0.5`.	`0.5`

Returns:

Type	Description
`numpy.ndarray of int, shape (n,)`	`+1` (long asset 1 / short asset 2), `-1` (short asset 1 / long asset 2) or `0` (no position) for each row.

Raises:

Type	Description
`ValueError`	If `entry` is not strictly between 0 and 0.5, `exit_band` is not in `[0, 0.5]`, or `u` does not have two columns.

Examples:

>>> import numpy as np
>>> from rcopula import ClaytonCopula
>>> from rcopula.portfolio import pairs_signal
>>> cop = ClaytonCopula(3.0)
>>> u = np.array([[0.02, 0.95], [0.95, 0.02], [0.50, 0.50]])
>>> pairs_signal(cop, u)
array([ 1, -1,  0])

With a narrower exit band the long is held while asset 1 is still cheap (row 2) and closed once it has reverted (row 3):

>>> u = np.array([[0.02, 0.95], [0.20, 0.80], [0.60, 0.50]])
>>> pairs_signal(cop, u, exit_band=0.1)
array([1, 1, 0])

Signals are rare by construction -- a few percent of observations:

>>> signals = pairs_signal(cop, cop.rvs(5000, random_state=0))
>>> bool(0.0 < np.mean(signals != 0) < 0.15)
True

backtest_pairs ¶

backtest_pairs(returns: ArrayLike, copula: Copula, train: int = 250, entry: float = 0.05, periods_per_year: int = 252, refit_every: int = 0, exit_band: float = 0.5) -> BacktestResult

Replay the copula pairs strategy through history, without look-ahead.

Use it to see how often :func:pairs_signal would have fired on a real pair and whether the trades made money before costs.

At each period t the copula is fitted on the trailing train returns, today's return is ranked within that window, and the :func:pairs_signal entry/exit rule (with entry and exit_band) updates the position, which is earned in the next period (t + 1). The position earned in period t + 1 therefore never depends on the return of period t + 1, which is the part most easily got wrong.

Parameters:

Name	Type	Description	Default
`returns`	`array_like or pandas.DataFrame of float, shape (n, 2)`	Period (e.g. daily) simple returns of the two assets, as fractions (0.01 = 1%). Column 0 is asset 1, column 1 is asset 2.	required
`copula`	`Copula`	Copula family to use, e.g. `ClaytonCopula()`. Its starting parameter is ignored; parameters are re-estimated from each training window.	required
`train`	`int`	Length of the trailing fitting window, in periods (250 is roughly one trading year of daily data).	`250`
`entry`	`float`	Entry threshold passed to :func:`pairs_signal`; must lie in (0, 0.5).	`0.05`
`periods_per_year`	`int`	Used only to annualise the Sharpe ratio (252 for daily, 52 weekly, 12 monthly).	`252`
`refit_every`	`int`	Refit cadence in periods. `0` refits every period; larger values are faster and make the position depend on a slightly staler model.	`0`
`exit_band`	`float`	Exit band passed to the :func:`pairs_signal` rule; in `[0, 0.5]`. The default closes every position after one period unless the entry rule fires again; smaller values hold positions until the pair has reverted.	`0.5`

Returns:

Type	Description
`BacktestResult`	Named tuple with `total_return`, `annualised_sharpe`, `n_trades`, `hit_rate`, `returns` (shape (n,)) and `positions` (shape (n,)); see :class:`BacktestResult`.

Raises:

Type	Description
`ValueError`	If `returns` is not two-column, has `train + 1` rows or fewer, or `entry` / `exit_band` are out of range.

Notes

Deliberately frictionless: no transaction costs, no borrow cost, no slippage, no capacity limit, and a spread traded at equal notional rather than a hedge ratio. Pairs strategies live or die on costs, so treat the Sharpe here as an upper bound on what the signal could deliver, not as a strategy result.

Examples:

>>> import numpy as np
>>> from rcopula import ClaytonCopula, CopulaDistribution
>>> from rcopula.portfolio import backtest_pairs
>>> from scipy import stats
>>> mv = CopulaDistribution(ClaytonCopula(4.0), [stats.norm(0, 0.01)] * 2)
>>> r = mv.rvs(600, random_state=0)
>>> result = backtest_pairs(r, ClaytonCopula(), train=250, refit_every=25)
>>> bool(result.n_trades >= 0 and np.isfinite(result.annualised_sharpe))
True

simulate_returns ¶

simulate_returns(copula: Copula, margins: Margin | list[Margin], n: int = 20000, random_state: Generator | int | None = None) -> NDArray[np.float64]

Simulate many joint return scenarios from a copula plus per-asset margins.

The copula sets how the assets move together (including in the tails); the margins set each asset's own return distribution. The output is a table of scenarios to feed :func:mean_cvar_weights or :func:efficient_frontier.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	A copula of dimension `d` with its parameters set.	required
`margins`	`scipy.stats frozen distribution or list of them, length d`	Return distribution of each asset, e.g. `stats.norm(0.0005, 0.012)` for a daily mean of 0.05% and volatility of 1.2%. A single distribution is used for every asset.	required
`n`	`int`	Number of scenarios (rows) to draw.	`20000`
`random_state`	`(int, Generator or None)`	Seed or generator for reproducible draws.	`None`

Returns:

Type	Description
`numpy.ndarray of float, shape (n, d)`	Simulated returns, one scenario per row and one asset per column, in the units of the margins.

Examples:

>>> from scipy import stats
>>> from rcopula import StudentCopula
>>> from rcopula.portfolio import simulate_returns
>>> r = simulate_returns(StudentCopula(0.4, dim=4, df=5), stats.norm(0.0005, 0.012),
...                      n=5000, random_state=0)
>>> r.shape
(5000, 4)

mean_cvar_weights ¶

mean_cvar_weights(scenarios: ArrayLike, alpha: float = 0.95, target_return: float | None = None, bounds: tuple[float, float] = (0.0, 1.0)) -> NDArray[np.float64]

Portfolio weights that minimise the average loss on the worst days.

Finds fully invested weights that minimise conditional value at risk (CVaR, also called expected shortfall): the average loss in the worst 1 - alpha share of scenarios. Unlike minimum variance, it looks only at the downside and sees tail dependence. Optionally require a minimum mean return to trace out a risk/return trade-off.

Parameters:

Name	Type	Description	Default
`scenarios`	`array_like of float, shape (n, d)`	Scenario returns (not losses), one scenario per row and one asset per column, e.g. from :func:`simulate_returns` or historical data.	required
`alpha`	`float`	CVaR confidence level in (0, 1). 0.95 averages the worst 5% of scenarios.	`0.95`
`target_return`	`float`	Minimum required mean portfolio return per scenario, in the same units as `scenarios`. Omit (`None`) for the global CVaR minimum.	`None`
`bounds`	`tuple of (float, float)`	`(lower, upper)` bound applied to every weight. The default forbids shorting; e.g. `(-1.0, 1.0)` allows it.	`(0.0, 1.0)`

Returns:

Type	Description
`numpy.ndarray of float, shape (d,)`	Weights, summing to one.

Raises:

Type	Description
`ValueError`	If `alpha` is not in (0, 1), or the problem is infeasible (most often a `target_return` above what any allowed portfolio can earn).

Notes

Rockafellar & Uryasev (2000) showed that minimising

.. math:: \mathrm{CVaR}\alpha(w) = \min\ z + \frac{1}{(1-\alpha)n}\sum_i \bigl(-r_i^{\top}w - z\bigr)^{+}

is a linear program in :math:(w, z, s) once returns are given as scenarios. That makes it exact and fast -- no gradient descent, no local optima -- and it is why scenario-based CVaR optimisation is practical at all.

The program has one auxiliary variable per scenario. Its constraint matrix is built as a scipy.sparse matrix (n * (d + 2) non-zeros) and solved with HiGHS, so memory grows linearly in n: tens of thousands of scenarios are routine.

Examples:

The optimiser avoids the asset with the fat left tail, even when its mean is identical:

>>> import numpy as np
>>> from rcopula.portfolio import mean_cvar_weights
>>> rng = np.random.default_rng(0)
>>> safe = rng.normal(0.001, 0.01, 4000)
>>> risky = rng.standard_t(2.5, 4000) * 0.006 + 0.001
>>> w = mean_cvar_weights(np.column_stack([safe, risky]))
>>> bool(w[0] > w[1])
True
>>> bool(abs(w.sum() - 1.0) < 1e-8)
True

min_variance_weights ¶

min_variance_weights(scenarios: ArrayLike, bounds: tuple[float, float] = (0.0, 1.0)) -> NDArray[np.float64]

Portfolio weights with the lowest volatility, as a Markowitz benchmark.

Provided for comparison with :func:mean_cvar_weights: mean-variance treats gains and losses symmetrically and sees only the covariance matrix, so it cannot distinguish two portfolios with the same covariance but different tail dependence.

Parameters:

Name	Type	Description	Default
`scenarios`	`array_like of float, shape (n, d)`	Scenario or historical returns, one row per scenario, one column per asset. Only their sample covariance is used.	required
`bounds`	`tuple of (float, float)`	`(lower, upper)` bound applied to every weight. The default forbids shorting.	`(0.0, 1.0)`

Returns:

Type	Description
`numpy.ndarray of float, shape (d,)`	Weights, summing to one, that minimise portfolio variance.

Raises:

Type	Description
`ValueError`	If the optimiser (SLSQP) fails to converge.

Examples:

>>> import numpy as np
>>> from rcopula.portfolio import min_variance_weights
>>> rng = np.random.default_rng(0)
>>> r = np.column_stack([rng.normal(0, 0.01, 2000), rng.normal(0, 0.03, 2000)])
>>> w = min_variance_weights(r)
>>> bool(w[0] > w[1])          # the quieter asset gets the weight
True

efficient_frontier ¶

efficient_frontier(scenarios: ArrayLike, alpha: float = 0.95, n_points: int = 20, bounds: tuple[float, float] = (0.0, 1.0)) -> EfficientFrontier

The best tail risk achievable at each level of expected return.

Solves :func:mean_cvar_weights for a ladder of target returns, from the minimum-CVaR portfolio's mean up to the highest single-asset mean, and reports the achieved mean return, CVaR and weights at each point -- the mean-CVaR efficient frontier. Plot cvars against achieved_returns to see how much extra tail loss each step of extra return costs.

Parameters:

Name	Type	Description	Default
`scenarios`	`array_like of float, shape (n, d)`	Scenario returns (not losses), one row per scenario, one column per asset.	required
`alpha`	`float`	CVaR confidence level in (0, 1).	`0.95`
`n_points`	`int`	Number of target returns to try along the frontier.	`20`
`bounds`	`tuple of (float, float)`	`(lower, upper)` bound applied to every weight. The default forbids shorting.	`(0.0, 1.0)`

Returns:

Type Description

EfficientFrontier

Named tuple (achieved_returns, cvars, weights):

achieved_returns : numpy.ndarray of float, shape (m,) Mean portfolio return achieved at each point (not the target that was requested; it is at least the target), in the units of scenarios. cvars : numpy.ndarray of float, shape (m,) CVaR at each point, as a positive loss. weights : numpy.ndarray of float, shape (m, d) Portfolio weights at each point; each row sums to one.

m <= n_points: infeasible targets are dropped rather than raising.

Raises:

Type	Description
`ValueError`	If `alpha` is not in (0, 1).

Examples:

>>> import numpy as np
>>> from rcopula.portfolio import efficient_frontier
>>> rng = np.random.default_rng(0)
>>> r = np.column_stack([rng.normal(0.0004, 0.01, 3000),
...                      rng.normal(0.0010, 0.02, 3000)])
>>> mu, cvar, w = efficient_frontier(r, n_points=8)
>>> bool(np.all(np.diff(cvar) >= -1e-9))     # more return costs more risk
True
>>> front = efficient_frontier(r, n_points=8)
>>> bool(np.allclose(front.achieved_returns, front.weights @ r.mean(axis=0)))
True

rcopula.statarb ¶

Choosing what to trade: pair and partner selection for statistical arbitrage.

:mod:rcopula.portfolio trades a pair once you have one. This picks it. Out of 500 names there are 124,750 pairs and 2.6 billion quadruples, so the selection rule is not a preliminary to the strategy -- for a copula strategy it largely is the strategy, and different rules surface genuinely different candidates.

Pairs. Six rules, of which the last two are the copula-specific ones:

=================== ========================================================= distance Sum of squared deviations between normalised cumulative returns (raw values, not ranks). The classic (Gatev et al.), and it selects for level agreement, which is not what a copula uses. pearson Linear correlation of the raw values. Included to be argued with: it is not invariant to the marginal transforms a copula discards. spearman Rank correlation. Cheap, and close to kendall. kendall Rank correlation, concordance-based. tail Lower-tail dependence, estimated nonparametrically. Picks pairs that crash together, which is where a copula strategy makes and loses its money. qq Mean distance from the diagonal of the QQ plot. Sensitive to disagreement anywhere in the distribution, not just the middle. =================== =========================================================

Partners. For a vine or a 4-dimensional copula you need a quadruple: a target and three partners. Exhaustive search over all triples of the remaining names is :math:O(n^3), so the four approaches below all reduce the candidate set first, exactly as Stübinger, Mangold and Krauss do.

=================== ========================================================= traditional Largest sum of pairwise Spearman correlations. extended Largest multivariate Spearman's rho -- three generalisations of the bivariate coefficient, averaged (Schmid and Schmidt 2007). Unlike the traditional sum, this is a genuine :math:d-dimensional measure rather than a pile of bivariate ones. geometric Smallest total distance from the hyper-diagonal of the unit hypercube. Purely geometric, no distributional assumption at all. extremal Strongest joint-tail concentration. Targets what Mangold's chi-squared test targets -- see the note under :func:select_partners about what this is and is not. =================== =========================================================

Almost everything here works on ranks, which is the point: a copula does not see the margins, so neither should the rule that chooses its inputs. The two exceptions are deliberate and are the conventional definitions of those rules: distance works on the raw values (it treats each column as a return series, cumulates it into a price path and normalises that path, as Gatev et al. do with normalised prices), and pearson is the linear correlation of the raw values. Ranking first would turn them into different rules -- pearson on ranks is spearman -- so they are kept as the literature defines them, as baselines to compare the rank-based rules with. Every partner-selection rule uses ranks.

============================ ================================================ :func:select_pairs Rank every pair by one of six criteria. :func:select_partners Find the best quadruple for a target. :func:multivariate_spearman Schmid and Schmidt's d-dimensional rho. :func:diagonal_distance Mean distance from the hyper-diagonal. :func:tail_concentration Joint-tail co-occurrence, relative to chance. ============================ ================================================

Examples:

>>> import numpy as np, pandas as pd, rcopula as rc
>>> from rcopula.statarb import select_pairs
>>> rng = np.random.default_rng(0)
>>> factor = rng.standard_normal(500)
>>> returns = pd.DataFrame({
...     "AAA": factor + 0.3 * rng.standard_normal(500),
...     "BBB": factor + 0.3 * rng.standard_normal(500),
...     "CCC": rng.standard_normal(500),
... })
>>> select_pairs(returns, method="kendall").iloc[0][["first", "second"]].tolist()
['AAA', 'BBB']

References

Gatev, E., Goetzmann, W. N. and Rouwenhorst, K. G. (2006). Pairs trading: performance of a relative-value arbitrage rule. Review of Financial Studies 19(3), 797-827. The distance rule. Stübinger, J., Mangold, B. and Krauss, C. (2018). Statistical arbitrage with vine copulas. Quantitative Finance 18(11), 1831-1849. The four partner-selection approaches. Schmid, F. and Schmidt, R. (2007). Multivariate extensions of Spearman's rho and related statistics. Statistics and Probability Letters 77(4), 407-416. The three multivariate rho formulas. Liew, R. Q. and Wu, Y. (2013). Pairs trading: a copula approach. J. Derivatives and Hedge Funds 19(1), 12-30. Krauss, C. (2017). Statistical arbitrage pairs trading strategies: review and outlook. J. Economic Surveys 31(2), 513-545.

multivariate_spearman ¶

multivariate_spearman(u: ArrayLike) -> float

Measure how strongly a whole group of columns moves together, as one number.

This is Schmid and Schmidt's multivariate Spearman's rho: a rank correlation for three or more variables at once. It is 0 when the columns are independent and 1 when they move in perfect lockstep (in rank terms), and equals ordinary Spearman's rho for two columns. :func:select_partners uses it for method="extended".

Parameters:

Name	Type	Description	Default
`u`	`array_like or pandas.DataFrame of float, shape (n, d)`	Observations, one column per variable, `d >= 2`. Raw data or pseudo-observations both work: ranks are always taken first.	required

Returns:

Type	Description
`float`	The average of the three coefficients below. 1 means perfect co-movement, 0 independence; negative values mean the columns tend to move against each other (the lower bound depends on `d`).

Raises:

Type	Description
`ValueError`	If `u` has fewer than 2 columns.

Notes

Three separate generalisations of the bivariate coefficient, averaged. Each reduces to ordinary Spearman's rho when :math:d = 2, and each measures something slightly different in higher dimensions -- the first two the concentration towards a corner, the third the average pairwise concordance -- so the mean is more robust than any one of them.

.. math::

\rho_1 &= h(d)\Bigl\{-1 + \tfrac{2^d}{n}
          \textstyle\sum_j \prod_i (1 - U_{ij})\Bigr\} \\
\rho_2 &= h(d)\Bigl\{-1 + \tfrac{2^d}{n}
          \textstyle\sum_j \prod_i U_{ij}\Bigr\} \\
\rho_3 &= -3 + \tfrac{12}{n\binom{d}{2}}
          \textstyle\sum_{k<l}\sum_j (1-U_{kj})(1-U_{lj})

with :math:h(d) = (d+1)/(2^d - d - 1).

Examples:

It agrees with the bivariate coefficient in two dimensions:

>>> import numpy as np, rcopula as rc
>>> from rcopula.statarb import multivariate_spearman
>>> u = rc.GaussianCopula(0.7).rvs(20000, random_state=0)
>>> bool(abs(multivariate_spearman(u) - rc.cor_spearman(u)[0, 1]) < 0.02)
True

It is zero under independence and one under comonotonicity, in any dimension:

>>> independent = rc.IndependenceCopula(4).rvs(20000, random_state=0)
>>> bool(abs(multivariate_spearman(independent)) < 0.02)
True
>>> comonotone = np.tile(np.random.default_rng(0).uniform(size=(20000, 1)), (1, 4))
>>> bool(multivariate_spearman(comonotone) > 0.99)
True

diagonal_distance ¶

diagonal_distance(u: ArrayLike) -> float

Measure how far ranked data sit from the "all move together" line (small = related).

Technically, the mean Euclidean distance of the pseudo-observations from the hyper-diagonal of the unit hypercube. :func:select_partners uses it (negated) for method="geometric".

Parameters:

Name	Type	Description	Default
`u`	`array_like or pandas.DataFrame of float, shape (n, d)`	Observations, one column per variable. Raw data or pseudo-observations both work: ranks are always taken first.	required

Returns:

Type	Description
`float`	Mean distance, `>= 0`. Zero means every row lies on the diagonal (perfect co-movement); independent columns give a clearly positive value (about 0.47 for `d = 4`).

Notes

Perfectly concordant observations lie exactly on the line from :math:(0,\dots,0) to :math:(1,\dots,1); independent ones scatter away from it. Small means strongly related, which is the opposite convention to every other measure here -- :func:select_partners handles the sign.

The appeal is that it assumes nothing: no distribution, no copula family, not even that the dependence is monotone in the usual sense.

Examples:

>>> import numpy as np, rcopula as rc
>>> from rcopula.statarb import diagonal_distance
>>> comonotone = np.tile(np.random.default_rng(0).uniform(size=(5000, 1)), (1, 4))
>>> float(diagonal_distance(comonotone)) < 1e-12
True
>>> independent = rc.IndependenceCopula(4).rvs(5000, random_state=0)
>>> bool(diagonal_distance(independent) > 0.3)
True

tail_concentration ¶

tail_concentration(u: ArrayLike, quantile: float | None = None) -> float

Measure how much more often all columns are extreme at the same time than chance would give.

Counts observations in the lower and upper corners of the unit hypercube and divides by what independence would put there. One means no more joint extremes than chance; larger means the names crash and spike together. This is the quantity a copula strategy lives or dies on, and no correlation measures it.

Parameters:

Name	Type	Description	Default
`u`	`array_like or pandas.DataFrame of float, shape (n, d)`	Observations, one column per variable. Raw data or pseudo-observations both work: ranks are always taken first.	required
`quantile`	`float or None`	Corner size `q`, strictly between 0 and 0.5: a row counts as jointly low if every rank is `<= q` and jointly high if every rank is `>= 1 - q`. `None` chooses it from the data as `min(0.25, (40 / n) ** (1 / d))`, so the corner is expected to hold about 40 rows under independence -- which matters more than it sounds; see the note below.	`None`

Returns:

Type	Description
`float`	Both corners combined, relative to independence: `(P(all low) + P(all high)) / (2 q**d)`. 1 means no more joint extremes than chance, larger means more, 0 means none observed.

Raises:

Type	Description
`ValueError`	If `quantile` is not strictly between 0 and 0.5.

Notes

A fixed corner cannot work beyond two dimensions. The corner holds :math:q^d of the mass, so at the conventional :math:q = 0.05 a 200,000-row sample expects 500 observations in it at :math:d = 2, 25 at :math:d = 3, and 1.3 at :math:d = 4 -- at which point the statistic is counting noise and scatters between 0.4 and 1.2 under independence.

The default therefore solves :math:q = (c/n)^{1/d} for a target count :math:c, so the corner always holds enough observations to mean something.

What that costs, stated plainly. A wide corner stops measuring the tail. At :math:d = 2 and 200,000 observations, Clayton and Gaussian at the same Kendall's tau give 2.15 against 2.13 at :math:q = 0.3 -- indistinguishable -- and only separate at :math:q \le 0.02, where the ratio reaches 1.2 and then 1.4. But :math:q = 0.02 at :math:d = 4 needs 250 million observations to keep 40 in the corner. Even a million rows only supports :math:q = 0.08.

So in four dimensions at any realistic sample size this ranks by joint co-movement in general, not by tail dependence specifically, and it cannot tell a heavy-tailed copula from a Gaussian one. That is a property of the data, not of the implementation: the observations required to resolve a four-dimensional tail simply are not there. Use it in two dimensions with an explicit small quantile when the tail is the question.

Examples:

>>> import rcopula as rc
>>> from rcopula.statarb import tail_concentration
>>> independent = rc.IndependenceCopula(4).rvs(40000, random_state=0)
>>> bool(abs(tail_concentration(independent) - 1.0) < 0.25)
True
>>> dependent = rc.ClaytonCopula(4.0, dim=4).rvs(40000, random_state=0)
>>> bool(tail_concentration(dependent) > 5)
True

select_pairs ¶

select_pairs(data: ArrayLike, method: PairMethod = 'kendall', *, top: int | None = None, names: list[str] | None = None) -> pd.DataFrame

Score every possible pair of columns and list them from most to least related.

The first step of a pairs-trading strategy: given returns for many instruments, find the pairs worth trading. Pick the criterion with method; the copula-specific ones are "tail" (pairs that crash together) and "qq".

Parameters:

Name	Type	Description	Default
`data`	`pandas.DataFrame or array_like of float, shape (n, k)`	Returns, one row per period and one column per instrument, `k >= 2`. DataFrame column names are carried through to the output.	required
`method`	`('distance', 'pearson', 'spearman', 'kendall', 'tail', 'qq')`	The selection criterion; see the module docstring for what each measures. `"distance"` and `"pearson"` use the raw values; the others use ranks.	`"distance"`
`top`	`int or None`	Keep only the best `top` pairs (a non-negative integer; `0` returns an empty frame with the usual columns). `None` keeps all `k(k-1)/2`.	`None, keyword-only`
`names`	`list of str or None`	Column names to use when `data` is a plain array (length `k`). Ignored for a DataFrame. `None` names the columns `x0, x1, ...`.	`None, keyword-only`

Returns:

Type Description

DataFrame

One row per pair, sorted best first, with columns:

first (column label): the first instrument of the pair.
second (column label): the second instrument.
score (float): the criterion value. Oriented so that larger is always better, whichever criterion was used -- the distance and QQ rules are negated for this reason, so their scores are non-positive.
rank (int): 1 for the best pair, 2 for the next, and so on.

out.attrs["method"] records the criterion used.

Raises:

Type	Description
`ValueError`	If `method` is not one of the six names, `data` has fewer than two columns, or `top` is negative or not an integer.

Notes

Cost is :math:O(k^2) pairs, each :math:O(n) for most rules. For kendall the whole matrix is computed at once by :func:~rcopula.cor_kendall, which for hundreds of names is seconds rather than the minutes that separate pairwise calls would take.

Examples:

>>> import numpy as np, pandas as pd
>>> from rcopula.statarb import select_pairs
>>> rng = np.random.default_rng(0)
>>> shared = rng.standard_normal(400)
>>> data = pd.DataFrame({
...     "A": shared + 0.2 * rng.standard_normal(400),
...     "B": shared + 0.2 * rng.standard_normal(400),
...     "C": rng.standard_normal(400),
...     "D": rng.standard_normal(400),
... })
>>> ranked = select_pairs(data, method="spearman")
>>> len(ranked)
6
>>> ranked.iloc[0]["first"], ranked.iloc[0]["second"]
('A', 'B')

The criteria disagree, which is the reason to have more than one:

>>> by_tail = select_pairs(data, method="tail").iloc[0]
>>> bool(by_tail["score"] > 0)
True

select_partners ¶

select_partners(data: ArrayLike, target: str | int, *, method: PartnerMethod = 'extended', n_partners: int = 3, n_candidates: int = 50, names: list[str] | None = None) -> dict[str, Any]

Find the few other columns that move most closely with a chosen target column.

Used to build a group (typically a target plus three partners) for a vine or four-dimensional copula trading strategy.

A vine or a :math:d-dimensional copula strategy needs a group, not a pair. Searching every triple of 499 remaining names is 20 million combinations per target, so this pre-screens to the n_candidates most correlated names and searches exhaustively within those -- which is what Stübinger, Mangold and Krauss do, and what makes the problem finite.

Parameters:

Name	Type	Description	Default
`data`	`pandas.DataFrame or array_like of float, shape (n, k)`	Returns, one row per period and one column per instrument. Needs at least `n_partners + 1` columns.	required
`target`	`str, int or other column label`	The column to find partners for. It is first looked up as a column label -- so with integer column labels an `int` names the label -- and only if it is not a label is an integer (Python or NumPy) read as a position, `0 <= target < k`.	required
`method`	`('traditional', 'extended', 'geometric', 'extremal')`	How a candidate group is scored: `"traditional"` sums the pairwise Spearman correlations, `"extended"` uses :func:`multivariate_spearman`, `"geometric"` uses :func:`diagonal_distance` (smaller is better), `"extremal"` uses :func:`tail_concentration`. See the module docstring.	`"traditional"`
`n_partners`	`int`	How many partners to choose, at least 1. Three gives the quadruple the literature uses.	`3, keyword-only`
`n_candidates`	`int`	Pre-screening width: only the `n_candidates` columns with the largest absolute rank correlation to the target are searched (never fewer than `n_partners`). Larger is more thorough and the work grows as :math:`\binom{n_{\text{candidates}}}{n_{\text{partners}}}`.	`50, keyword-only`
`names`	`list of str or None`	Column names to use when `data` is a plain array. Ignored for a DataFrame. `None` names the columns `x0, x1, ...`.	`None, keyword-only`

Returns:

Type Description

dict

With keys:

"target" (column label): the target, as a label even if given by position.
"partners" (list of column labels): the best partners, length n_partners.
"score" (float): the winning group's score, oriented so larger is better ("geometric" is negated, so it is non-positive).
"method" (str): the method used.
"considered" (int): how many combinations were actually evaluated, so the pre-screening is visible rather than implicit.

Raises:

Type	Description
`ValueError`	If `method` is not one of the four names, `target` is neither a column label nor a valid position, `n_partners < 1`, or there are fewer than `n_partners` columns besides the target.

Notes

On extremal, two caveats, both worth reading before relying on it.

Stübinger et al. use Mangold's (2015) chi-squared test against independence within a specific parametric copula family. This uses joint-tail concentration instead -- the same property, a different statistic, arrived at without that paper. It is not a reimplementation of Mangold's test and produces no p-value.

More importantly, in four dimensions it is not really an extremal rule at all. A four-dimensional corner deep enough to isolate the tail is empty at any realistic sample size -- see :func:tail_concentration -- so the corner widens until it holds data, and a wide corner measures joint co-movement in general. Expect it to behave like a more robust traditional, and do not expect it to find tail dependence that the other rules miss.

Examples:

>>> import numpy as np, pandas as pd
>>> from rcopula.statarb import select_partners
>>> rng = np.random.default_rng(0)
>>> factor = rng.standard_normal(600)
>>> frame = pd.DataFrame({
...     "TGT": factor + 0.3 * rng.standard_normal(600),
...     "P1": factor + 0.3 * rng.standard_normal(600),
...     "P2": factor + 0.3 * rng.standard_normal(600),
...     "P3": factor + 0.3 * rng.standard_normal(600),
...     "N1": rng.standard_normal(600),
...     "N2": rng.standard_normal(600),
... })
>>> found = select_partners(frame, "TGT", method="extended")
>>> sorted(found["partners"])
['P1', 'P2', 'P3']

Copula-GARCH margins¶

fit_garch and CopulaGarch.fit take two model switches. vol="gjr" adds the GJR leverage term (volatility rises more after a fall than after a rally of the same size); mean="zero", "ar1" or "arma11" replaces the constant mean. The defaults (vol="garch", mean="constant") are the GARCH(1,1) of earlier releases, with identical results. Fits are checked against R's rugarch (tools/rgolden/10_garch.R).

from rcopula.garch import fit_garch

res = fit_garch(returns, dist="t", vol="gjr", mean="ar1")
res.gamma, res.phi  # leverage and AR(1) coefficients
res.forecast_vol(10)  # accounts for gamma: persistence = alpha + gamma/2 + beta
res.forecast_mean(10)  # AR(1) mean decaying to mu / (1 - phi)

rcopula.garch ¶

Copula-GARCH: dependence in the cross-section, volatility in time.

Fitting a copula directly to asset returns is almost always wrong. Returns are not identically distributed -- they arrive in volatile and quiet regimes -- and a copula fitted to the raw series confuses volatility clustering with dependence. Two assets that are independent but share a calm month and a wild month will look strongly tail-dependent, because both were large at the same time for a reason that has nothing to do with either.

The copula-GARCH model separates the two:

.. math::

r_{j,t} = \mu_j + \sigma_{j,t} z_{j,t}, \qquad
\sigma_{j,t}^2 = \omega_j + \alpha_j \varepsilon_{j,t-1}^2
                 + \beta_j \sigma_{j,t-1}^2,

with the innovation vector :math:(z_{1,t}, \dots, z_{d,t}) i.i.d. across time and coupled by a copula. Each margin gets its own volatility dynamics; the copula then describes what is left, which is dependence proper.

This is a two-step estimator (Patton 2006): fit each GARCH by quasi-maximum-likelihood, take the standardised residuals, and fit the copula to their pseudo-observations. The second step inherits the first step's estimation error, which is why the copula standard errors reported here are conditional on the fitted margins.

Forecasting is where the model earns its place. Simulating forward gives a joint predictive distribution of returns over any horizon, from which portfolio VaR, expected shortfall or an option payoff follows directly -- with tail dependence and volatility persistence both present, which no single correlation number can deliver.

============================ ================================================ :func:fit_garch GARCH(1,1) or GJR-GARCH(1,1) by QMLE, with a constant, zero, AR(1) or ARMA(1,1) mean, in pure NumPy. :class:GarchResult A fitted margin, with forecasting. :class:CopulaGarch The joint model: GARCH margins plus a copula. ============================ ================================================

The GARCH implementation here is deliberately small -- GARCH(1,1) or its leverage-effect variant GJR-GARCH(1,1) (Glosten, Jagannathan & Runkle 1993), a constant, zero, AR(1) or ARMA(1,1) mean, and normal or Student-t innovations -- because that is what the copula literature uses and it keeps rcopula dependency-free. With vol="gjr" the variance reacts more to bad news than to good,

.. math::

\sigma_t^2 = \omega + (\alpha + \gamma\,\mathbf{1}[\varepsilon_{t-1} < 0])
             \,\varepsilon_{t-1}^2 + \beta\sigma_{t-1}^2,

which is the asymmetry almost every equity index shows. For EGARCH, long-memory or regime-switching margins, fit them with the arch package and pass the standardised residuals to :func:~rcopula.fit.fit yourself; the second step is unchanged.

References

Patton, A. J. (2006). Modelling asymmetric exchange rate dependence. International Economic Review 47(2), 527-556. The two-step copula-GARCH estimator. Jondeau, E. and Rockinger, M. (2006). The copula-GARCH model of conditional dependencies: an international stock market application. Journal of International Money and Finance 25(5), 827-853. Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics 31(3), 307-327. Glosten, L. R., Jagannathan, R. and Runkle, D. E. (1993). On the relation between the expected value and the volatility of the nominal excess return on stocks. Journal of Finance 48(5), 1779-1801. The GJR (threshold) GARCH model, vol="gjr". Bollerslev, T. and Wooldridge, J. M. (1992). Quasi-maximum likelihood estimation and inference in dynamic models with time-varying covariances. Econometric Reviews 11(2), 143-172. Why the normal-innovation GARCH estimate stays consistent under misspecification -- the "Q" in QMLE. Barone-Adesi, G., Giannopoulos, K. and Vosper, L. (1999). VaR without correlations for portfolios of derivative securities. Journal of Futures Markets 19(5), 583-602. Filtered historical simulation, the innovations="empirical" option.

GarchResult `dataclass` ¶

GarchResult(mu: float, omega: float, alpha: float, beta: float, df: float | None, sigma: NDArray[float64], resid: NDArray[float64], loglik: float, dist: str, name: str = '', gamma: float = 0.0, phi: float = 0.0, theta: float = 0.0, vol: str = 'garch', mean: str = 'constant', cond_mean: NDArray[float64] | None = None)

One asset's fitted volatility model: parameters, per-period volatility, forecasts.

This is what :func:fit_garch returns. It tells you how volatile the asset is today, how volatile it is on average, how quickly a volatility spike fades, and it holds the "de-volatilised" returns (the standardised residuals) that the copula step is fitted to. You normally get one from :func:fit_garch rather than building it yourself.

The model is :math:r_t = m_t + \varepsilon_t, :math:\varepsilon_t = \sigma_t z_t, with conditional mean :math:m_t = \mu + \phi r_{t-1} + \theta\varepsilon_{t-1} and variance :math:\sigma_t^2 = \omega + (\alpha + \gamma\mathbf{1}[\varepsilon_{t-1}<0]) \varepsilon_{t-1}^2 + \beta\sigma_{t-1}^2. The default model has :math:\phi = \theta = \gamma = 0: a constant-mean GARCH(1,1).

Parameters:

Name	Type	Description	Default
`mu`	`float`	Intercept of the mean equation per period, in the units of the input returns. For the constant-mean model it is the mean return; with an AR term the long-run mean is `mu / (1 - phi)` (see :attr:`unconditional_mean`). `0.0` for `mean="zero"`.	required
`omega`	`float`	Variance intercept :math:`\omega`, in squared return units.	required
`alpha`	`float`	Reaction to the previous period's shock (ARCH coefficient), between 0 and 1.	required
`beta`	`float`	Weight on the previous period's variance (GARCH coefficient), between 0 and 1.	required
`df`	`float or None`	Student-t degrees of freedom of the innovations; `None` for normal innovations.	required
`sigma`	`numpy.ndarray of float, shape (n,)`	Fitted conditional standard deviation for each observation.	required
`resid`	`numpy.ndarray of float, shape (n,)`	Standardised residuals, unitless.	required
`loglik`	`float`	Maximised log-likelihood on the original data scale.	required
`dist`	`str`	`"normal"` or `"t"`.	required
`name`	`str`	Series label.	`""`
`gamma`	`float`	Extra reaction to a negative previous shock (the GJR leverage coefficient). `0.0` for plain GARCH.	`0.0`
`phi`	`float`	AR(1) coefficient on the previous return, between -1 and 1.	`0.0`
`theta`	`float`	MA(1) coefficient on the previous shock, between -1 and 1.	`0.0`
`vol`	`str`	Variance model: `"garch"` or `"gjr"`.	`"garch"`
`mean`	`str`	Mean model: `"constant"`, `"zero"`, `"ar1"` or `"arma11"`.	`"constant"`
`cond_mean`	`numpy.ndarray of float, shape (n,)`	Fitted conditional mean :math:`m_t` for each observation. Needed to forecast an AR/ARMA mean; when omitted it is taken to be `mu` throughout.	`None`

Attributes:

Name	Type	Description
`mu, omega, alpha, beta, gamma`	`float`	Mean intercept and variance-equation parameters, on the original scale of the data (`mu` in return units, `omega` in squared return units; `alpha`, `beta` and `gamma` are unitless).
`phi, theta`	`float`	AR(1) and MA(1) coefficients of the mean, unitless; `0.0` when the mean model does not have them.
`df`	`float or None`	Innovation degrees of freedom; `None` for normal innovations.
`sigma`	`numpy.ndarray of float, shape (n,)`	Fitted conditional standard deviations, one per observation, in the units of the input returns (e.g. daily volatility for daily returns).
`resid`	`numpy.ndarray of float, shape (n,)`	Standardised residuals :math:`z_t = (r_t - m_t)/\sigma_t`: each period's surprise divided by that period's volatility. These are the input to the copula step.
`loglik`	`float`	Maximised log-likelihood.
`dist`	`str`	`"normal"` or `"t"`.
`name`	`str`	Series label, carried through from a `pandas` column name; `""` if none was given.
`vol`	`str`	`"garch"` or `"gjr"`.
`mean`	`str`	`"constant"`, `"zero"`, `"ar1"` or `"arma11"`.
`cond_mean`	`numpy.ndarray of float, shape (n,), or None`	Fitted conditional means, in the units of the input returns.

persistence `property` ¶

persistence: float

How long volatility shocks linger: :math:\alpha + \gamma/2 + \beta.

Values near 1 (typical for daily equity returns: 0.97-0.99) mean a volatility spike fades slowly; values well below 1 mean it dies out within a few periods. For plain GARCH (gamma = 0) this is :math:\alpha + \beta. The GJR term counts half because a symmetric innovation is negative half the time.

Returns:

Type	Description
`float`	:math:`\alpha + \gamma/2 + \beta`, between 0 and 1 (capped at 0.9999 by the fit).

unconditional_mean `property` ¶

unconditional_mean: float

The long-run average return the mean reverts to: :math:\mu/(1-\phi).

Equal to mu unless the mean has an AR term.

Returns:

Type	Description
`float`	Per-period mean return, in the units of the input returns.

unconditional_vol `property` ¶

unconditional_vol: float

The long-run average volatility the model reverts to.

Technically the unconditional standard deviation of the shocks, :math:\sqrt{\omega/(1-\alpha-\gamma/2-\beta)}. Forecasts drift towards it as the horizon grows.

Returns:

Type	Description
`float`	Per-period standard deviation, in the units of the input returns (e.g. daily volatility for daily returns; multiply by :math:`\sqrt{252}` for an annualised figure).

half_life `property` ¶

half_life: float

Number of periods for a volatility shock to fade by half.

Computed as log(0.5) / log(persistence). The unit is one period of the input data: days for daily returns, weeks for weekly returns, and so on.

Returns:

Type	Description
`float`	Half-life in periods of the input data. `0.0` when `persistence == 0` (no persistence: a shock is gone by the next period), and `inf` when `persistence >= 1` (shocks never fade).

n_params `property` ¶

n_params: int

Number of estimated parameters, as used by :attr:aic and :attr:bic.

The default model has 4 (mu, omega, alpha, beta), plus 1 for Student-t df. The GJR gamma adds 1; the mean adds 0 ("zero"), 1 ("constant"), 2 ("ar1": mu, phi) or 3 ("arma11": mu, phi, theta) in place of the constant's 1.

Returns:

Type	Description
`int`	The parameter count, between 3 and 8.

aic `property` ¶

aic: float

Akaike information criterion; lower is better when comparing fits.

Use it to compare, say, a normal and a Student-t fit of the same series. Computed as 2 * n_params - 2 * loglik.

Returns:

Type	Description
`float`	The AIC.

bic `property` ¶

bic: float

Bayesian information criterion; lower is better, and stricter than AIC.

Computed as n_params * log(n) - 2 * loglik, with n the number of observations.

Returns:

Type	Description
`float`	The BIC.

innovation ¶

innovation() -> Any

The distribution of each period's volatility-adjusted "surprise" (variance 1).

Standard normal for dist="normal"; a Student-t rescaled to unit variance for dist="t". It is a frozen scipy.stats distribution, so it plugs straight into :class:~rcopula.distribution.CopulaDistribution.

Returns:

Type	Description
`scipy.stats frozen distribution`	Mean 0, variance 1, with `.ppf`, `.cdf`, `.rvs` and so on.

forecast_variance ¶

forecast_variance(horizon: int = 1) -> NDArray[np.float64]

Forecast the variance (volatility squared) for each of the next periods.

Use it to see how today's elevated (or depressed) volatility is expected to drift back to its long-run level over the coming periods.

Parameters:

Name	Type	Description	Default
`horizon`	`int`	Number of periods ahead to forecast; must be at least 1.	`1`

Returns:

Type	Description
`numpy.ndarray of float, shape (horizon,)`	Element `h - 1` is the expected variance `h` periods after the last observation, in squared return units. These are per-period variances, not cumulative ones.

Raises:

Type	Description
`ValueError`	If `horizon < 1`.

Notes

One step ahead is exact (for GJR it uses the sign of the last shock); beyond that the forecast decays geometrically towards the unconditional variance,

.. math:: \mathbb{E}[\sigma_{n+h}^2] = \bar\sigma^2 + (\alpha+\gamma/2+\beta)^{h-1}\bigl(\sigma_{n+1}^2 - \bar\sigma^2\bigr),

because a future shock is negative with probability 1/2 whatever its size, so :math:\mathbb{E}[\mathbf{1}[\varepsilon<0]\varepsilon^2] = \sigma^2/2 for normal and Student-t innovations alike.

Examples:

>>> import numpy as np
>>> from rcopula.garch import fit_garch
>>> rng = np.random.default_rng(0)
>>> x = rng.standard_normal(1500) * 0.01
>>> res = fit_garch(x)
>>> v = res.forecast_variance(250)
>>> bool(abs(np.sqrt(v[-1]) / res.unconditional_vol - 1) < 0.05)
True

forecast_mean ¶

forecast_mean(horizon: int = 1) -> NDArray[np.float64]

Forecast the expected return for each of the next periods.

Constant (or zero) for the constant- and zero-mean models. With an AR(1) or ARMA(1,1) mean the first forecast uses the last observed return and shock, and later ones decay geometrically (at rate phi) towards :attr:unconditional_mean.

Parameters:

Name	Type	Description	Default
`horizon`	`int`	Number of periods ahead to forecast; must be at least 1.	`1`

Returns:

Type	Description
`numpy.ndarray of float, shape (horizon,)`	Element `h - 1` is the expected return `h` periods after the last observation, in the units of the input returns.

Raises:

Type	Description
`ValueError`	If `horizon < 1`.

Notes

.. math:: \mathbb{E}[r_{n+1}] = \mu + \phi r_n + \theta\varepsilon_n, \qquad \mathbb{E}[r_{n+h}] = \mu + \phi\,\mathbb{E}[r_{n+h-1}] \quad (h \ge 2).

Examples:

>>> import numpy as np
>>> from rcopula.garch import fit_garch
>>> rng = np.random.default_rng(0)
>>> x = np.zeros(3000)
>>> for t in range(1, 3000):
...     x[t] = 0.1 + 0.5 * x[t - 1] + rng.standard_normal()
>>> res = fit_garch(x, mean="ar1")
>>> m = res.forecast_mean(200)
>>> bool(abs(m[-1] - res.unconditional_mean) < 1e-9)
True
>>> bool(abs(res.unconditional_mean - 0.2) < 0.1)
True

forecast_vol ¶

forecast_vol(horizon: int = 1) -> NDArray[np.float64]

Forecast the volatility (standard deviation) for each of the next periods.

The square root of :meth:forecast_variance.

Parameters:

Name	Type	Description	Default
`horizon`	`int`	Number of periods ahead to forecast; must be at least 1.	`1`

Returns:

Type	Description
`numpy.ndarray of float, shape (horizon,)`	Per-period volatility forecasts, in the units of the input returns.

Raises:

Type	Description
`ValueError`	If `horizon < 1`.

Examples:

>>> import numpy as np
>>> from rcopula.garch import fit_garch
>>> rng = np.random.default_rng(0)
>>> res = fit_garch(rng.standard_normal(1000) * 0.01)
>>> res.forecast_vol(5).shape
(5,)

CopulaGarch ¶

CopulaGarch(margins: list[GarchResult], copula: Copula, innovations: Literal['empirical', 'parametric'] = 'empirical')

Multi-asset returns: each asset's own volatility, plus a copula linking the shocks.

Each asset gets a GARCH(1,1) model -- optionally GJR, with an AR or ARMA mean (see :func:fit_garch) -- so that its volatility can rise and fall over time; the copula then describes how the assets' de-volatilised shocks move together, including in the tails. Use it to simulate future joint returns and to produce forward-looking portfolio VaR and expected shortfall that reflect both today's volatility and crash co-movement.

Most users build one with :meth:fit from a returns table. The constructor is for assembling a model from margins and a copula you have already fitted (or chosen) yourself.

Parameters:

Name	Type	Description	Default
`margins`	`list of GarchResult, length d`	One fitted GARCH per series, e.g. from :func:`fit_garch`.	required
`copula`	`Copula`	Copula for the standardised innovations, with `copula.dim == d`. It is used as given (not refitted).	required
`innovations`	`('empirical', 'parametric')`	How to invert the copula's uniforms when simulating. `"empirical"` draws from the observed standardised residuals (filtered historical simulation) and so inherits their skew and kurtosis without assuming a shape; `"parametric"` uses the fitted normal or Student-t. Empirical cannot produce an innovation larger than the largest one seen, so use parametric for long-horizon or deep-tail work.	`"empirical"`

Attributes:

Name	Type	Description
`margins`	`list of GarchResult, length d`	The per-asset volatility models.
`copula`	`Copula`	The copula coupling the standardised innovations.
`innovations`	`str`	`"empirical"` or `"parametric"`.

Raises:

Type	Description
`ValueError`	If `len(margins) != copula.dim`, or `innovations` is not one of the two allowed strings.

Examples:

See :meth:fit.

dim `property` ¶

dim: int

Number of assets in the model.

Returns:

Type	Description
`int`	`d`, the number of margins.

names `property` ¶

names: list[str]

Asset labels, falling back to "x0", "x1", ... for unnamed margins.

Returns:

Type	Description
`list of str, length d`	One label per asset, in column order.

fit `classmethod` ¶

fit(returns: ArrayLike, copula: Copula, dist: Literal['normal', 't'] = 'normal', innovations: Literal['empirical', 'parametric'] = 'empirical', method: str = 'mpl', vol: Literal['garch', 'gjr'] = 'garch', mean: Literal['constant', 'zero', 'ar1', 'arma11'] = 'constant') -> CopulaGarch

Fit the whole model to a table of asset returns.

Two-step estimation (Patton 2006): first a GARCH(1,1) is fitted to each column to remove volatility clustering, then the copula is fitted to the resulting standardised residuals. This is the usual entry point.

Parameters:

Name	Type	Description	Default
`returns`	`array_like of float or pandas.DataFrame, shape (n, d)`	Returns, one row per period (oldest first) and one column per asset, with at least 50 rows. Column names are kept if a frame is passed.	required
`copula`	`Copula`	Family to fit, with `dim` matching the number of columns. Any starting parameters are ignored -- it is refitted.	required
`dist`	`('normal', 't')`	Innovation distribution for the marginal GARCH models; see :func:`fit_garch`.	`"normal"`
`innovations`	`('empirical', 'parametric')`	Simulation margin; see the class docstring.	`"empirical"`
`method`	`str`	Copula estimation method, passed to :func:`~rcopula.fit.fit`. The default `"mpl"` is the standard choice here, since the residuals' distribution is not being claimed to be exactly the fitted one.	`"mpl"`
`vol`	`('garch', 'gjr')`	Variance model for every margin; `"gjr"` adds the leverage effect. See :func:`fit_garch`.	`"garch"`
`mean`	`('constant', 'zero', 'ar1', 'arma11')`	Mean model for every margin. See :func:`fit_garch`.	`"constant"`

Returns:

Type	Description
`CopulaGarch`	The fitted joint model.

Raises:

Type	Description
`ValueError`	If `returns` is not 2-d, its column count differs from `copula.dim`, or any column fails :func:`fit_garch` (too short, non-finite or constant).

Examples:

Two independent series driven by a common volatility process. Filtering first is what stops the shared volatility being read as tail dependence:

>>> import numpy as np
>>> import rcopula as rc
>>> from rcopula.garch import CopulaGarch
>>> rng = np.random.default_rng(1)
>>> h, e = np.zeros(3000), rng.standard_normal(3000)
>>> for t in range(1, 3000):
...     h[t] = 0.98 * h[t - 1] + 0.25 * e[t]
>>> r = rng.standard_normal((3000, 2)) * (0.01 * np.exp(h))[:, None]
>>> model = CopulaGarch.fit(r, rc.StudentCopula(0.0, df=8.0, dim=2))
>>> raw = rc.fit(rc.StudentCopula(0.0, df=8.0), rc.pseudo_obs(r), method="mpl")
>>> bool(raw.copula.lambda_().upper > 3 * model.copula.lambda_().upper)
True

Note that the artifact is invisible to rank correlation -- both are near zero -- which is precisely why it goes unnoticed:

>>> bool(abs(raw.copula.params[0]) < 0.06 and abs(model.copula.params[0]) < 0.06)
True

The fitted model then forecasts a joint distribution of returns:

>>> paths = model.simulate(horizon=10, n=2000, random_state=0)
>>> paths.shape
(2000, 10, 2)

simulate ¶

simulate(horizon: int = 1, n: int = 10000, random_state: Generator | int | None = None) -> NDArray[np.float64]

Generate n possible future return paths for all assets, starting from today.

Use it for scenario analysis or any risk measure the built-in :meth:forecast_risk does not cover (e.g. a path-dependent payoff).

Each step draws an innovation vector from the copula -- so the cross-sectional dependence is the fitted one -- and pushes it through each margin's own GARCH recursion, so volatility keeps clustering along the path. Both effects are present simultaneously, which is the reason to build the model at all. Paths start from the last observed volatility, shock and (for an AR/ARMA mean) return.

Each margin is simulated with the model it was fitted with: the GJR recursion when vol="gjr" (negative shocks feed more into the next period's variance) and the AR(1)/ARMA(1,1) mean recursion :math:r_t = \mu + \phi r_{t-1} + \theta\varepsilon_{t-1} + \varepsilon_t when there is one. With innovations="empirical" a GJR margin's variance follows the observed residuals' own share of negative squared shocks rather than exactly one half, which is the point of filtered historical simulation.

Parameters:

Name	Type	Description	Default
`horizon`	`int`	Number of future periods per path; must be at least 1.	`1`
`n`	`int`	Number of simulated paths.	`10_000`
`random_state`	`(int, Generator or None)`	Seed or generator, for reproducible draws.	`None`

Returns:

Type	Description
`numpy.ndarray of float, shape (n, horizon, d)`	`out[i, t, j]` is the return of asset `j` in period `t + 1` of path `i`, in the units of the input returns.

Raises:

Type	Description
`ValueError`	If `horizon < 1` or `n < 1`.

Examples:

Simulated innovations carry the copula's dependence:

>>> import numpy as np
>>> from scipy import stats
>>> import rcopula as rc
>>> from rcopula.garch import CopulaGarch, fit_garch
>>> rng = np.random.default_rng(0)
>>> r = rng.standard_normal((1200, 2)) * 0.01
>>> margins = [fit_garch(r[:, j]) for j in range(2)]
>>> model = CopulaGarch(margins, rc.ClaytonCopula.from_tau(0.5))
>>> paths = model.simulate(horizon=1, n=8000, random_state=0)
>>> tau = stats.kendalltau(paths[:, 0, 0], paths[:, 0, 1]).statistic
>>> bool(abs(tau - 0.5) < 0.03)
True

forecast ¶

forecast(horizon: int = 1, n: int = 10000, random_state: Generator | int | None = None) -> NDArray[np.float64]

Simulate each asset's total return over the next horizon periods.

Runs :meth:simulate and adds up each path's per-period returns. Sums the simulated log-returns, which is the usual convention. For simple returns compound them instead (simulate and use np.prod(1 + paths, axis=1) - 1).

Parameters:

Name	Type	Description	Default
`horizon`	`int`	Number of future periods to accumulate over; must be at least 1.	`1`
`n`	`int`	Number of simulated scenarios.	`10_000`
`random_state`	`(int, Generator or None)`	Seed or generator, for reproducible draws.	`None`

Returns:

Type	Description
`numpy.ndarray of float, shape (n, d)`	Cumulative return of each asset in each scenario.

Raises:

Type	Description
`ValueError`	If `horizon < 1` or `n < 1`.

Examples:

>>> import numpy as np
>>> import rcopula as rc
>>> from rcopula.garch import CopulaGarch
>>> rng = np.random.default_rng(0)
>>> r = rng.standard_normal((1000, 2)) * 0.01
>>> model = CopulaGarch.fit(r, rc.GaussianCopula(0.0, dim=2))
>>> model.forecast(horizon=5, n=100, random_state=0).shape
(100, 2)

forecast_risk ¶

forecast_risk(weights: ArrayLike | None = None, alpha: float = 0.99, horizon: int = 1, n: int = 50000, random_state: Generator | int | None = None) -> dict[str, float]

Forecast a portfolio's Value-at-Risk and expected shortfall over the horizon.

This is the payoff of the whole construction: a forward-looking risk number that respects both current volatility -- which a static copula ignores -- and tail dependence, which a GARCH-only model ignores. Estimated by Monte Carlo from :meth:forecast.

Parameters:

Name	Type	Description	Default
`weights`	`array_like of float, shape (d,)`	Portfolio exposures, one per asset in column order (e.g. `[0.6, 0.4]`). They are used as given and not rescaled to sum to 1, so they can describe a leveraged, long-short or partly invested book: the portfolio return is `asset_returns @ weights`, and doubling the weights doubles every statistic. Equal-weighted (`1 / d` each) if omitted.	`None`
`alpha`	`float`	Confidence level, e.g. 0.99 for a 99% VaR.	`0.99`
`horizon`	`int`	Forecast horizon in periods; must be at least 1.	`1`
`n`	`int`	Number of simulated scenarios. More scenarios give a less noisy estimate, especially at high `alpha`.	`50_000`
`random_state`	`(int, Generator or None)`	Seed or generator, for reproducible results.	`None`

Returns:

Type Description

dict of str to float

Statistics of the portfolio return over the horizon, in the units of the input returns:

"var" -- Value-at-Risk, as a loss (positive means losing money): exceeded with probability 1 - alpha.
"expected_shortfall" -- average loss given the VaR is exceeded, as a loss.
"mean" -- mean portfolio return (a return, not a loss).
"volatility" -- sample standard deviation (ddof=1) of the simulated portfolio return.

Raises:

Type	Description
`ValueError`	If `weights` does not have length `d`, `horizon < 1` or `n < 2`.

Examples:

Tail dependence raises the risk number at identical Kendall tau, which is exactly the comparison a correlation-based model cannot make:

>>> import numpy as np
>>> import rcopula as rc
>>> from rcopula.garch import CopulaGarch, fit_garch
>>> rng = np.random.default_rng(0)
>>> r = rng.standard_normal((1500, 2)) * 0.01
>>> margins = [fit_garch(r[:, j]) for j in range(2)]
>>> gauss = CopulaGarch(margins, rc.GaussianCopula.from_tau(0.5))
>>> clayton = CopulaGarch(margins, rc.ClaytonCopula.from_tau(0.5))
>>> a = gauss.forecast_risk(alpha=0.99, n=40_000, random_state=0)
>>> b = clayton.forecast_risk(alpha=0.99, n=40_000, random_state=0)
>>> bool(b["var"] > a["var"])
True

summary ¶

summary() -> pd.DataFrame

A table of each asset's fitted volatility parameters, one row per asset.

Returns:

Type	Description
`pandas.DataFrame, shape (d, 8) to (d, 11)`	Indexed by :attr:`names`, with float columns `mu`, `omega`, `alpha`, `beta`, `df` (NaN for normal innovations), `persistence`, `half_life` (periods) and `loglik`. Columns `phi` and `theta` (after `mu`) appear when any margin has an AR/ARMA mean, and `gamma` (after `beta`) when any margin is GJR.

Examples:

>>> import numpy as np
>>> import rcopula as rc
>>> from rcopula.garch import CopulaGarch
>>> rng = np.random.default_rng(0)
>>> r = rng.standard_normal((1000, 2)) * 0.01
>>> model = CopulaGarch.fit(r, rc.GaussianCopula(0.0, dim=2))
>>> list(model.summary().columns)
['mu', 'omega', 'alpha', 'beta', 'df', 'persistence', 'half_life', 'loglik']

fit_garch ¶

fit_garch(x: ArrayLike, dist: Literal['normal', 't'] = 'normal', name: str = '', vol: Literal['garch', 'gjr'] = 'garch', mean: Literal['constant', 'zero', 'ar1', 'arma11'] = 'constant') -> GarchResult

Estimate one asset's time-varying volatility from its return history.

Fits a GARCH(1,1) model -- by default with a constant mean -- by (quasi-)maximum likelihood. Use it to measure how volatile a series is today versus on average, to forecast volatility, or -- the main use in this package -- to strip volatility clustering out of returns before fitting a copula (see :class:CopulaGarch). Options add a leverage effect (GJR-GARCH) and an AR(1) or ARMA(1,1) mean; the mean and variance are estimated jointly.

Parameters:

Name	Type	Description	Default
`x`	`array_like of float, shape (n,)`	A single return series (e.g. daily log-returns), oldest first, with at least 50 finite observations. Any units work (0.01 or 1.0 for 1%); the fitted `mu`, `omega` and `sigma` come back in the same units. Multi-dimensional input is flattened.	required
`dist`	`('normal', 't')`	Innovation distribution. `"normal"` is quasi-MLE -- consistent for the variance parameters even when returns are fat-tailed (Bollerslev & Wooldridge 1992), which is why it remains the default. `"t"` estimates the degrees of freedom as well and gives a better fit when you intend to simulate from the margin rather than only filter with it.	`"normal"`
`name`	`str`	Label stored on the result (used as a row name by :meth:`CopulaGarch.summary`).	`""`
`vol`	`('garch', 'gjr')`	Variance model. `"garch"` is the symmetric GARCH(1,1). `"gjr"` (Glosten-Jagannathan-Runkle) adds a coefficient `gamma` on squared negative shocks, so bad news can raise volatility more than good news of the same size -- the leverage effect typical of equities. It is constrained to `alpha >= 0`, `alpha + gamma >= 0` (`gamma` may be negative, for series where good news is the more volatile) and `alpha + gamma/2 + beta < 1` (stationarity).	`"garch"`
`mean`	`('constant', 'zero', 'ar1', 'arma11')`	Mean model for :math:`m_t` in :math:`r_t = m_t + \varepsilon_t`. `"constant"`: :math:`m_t = \mu`. `"zero"`: :math:`m_t = 0` (no mean is estimated; common for daily returns whose mean is indistinguishable from zero). `"ar1"`: :math:`m_t = \mu + \phi r_{t-1}`. `"arma11"`: :math:`m_t = \mu + \phi r_{t-1} + \theta\varepsilon_{t-1}`. The AR and MA coefficients are kept inside `(-1, 1)` so the mean is stationary and invertible.	`"constant"`

Returns:

Type	Description
`GarchResult`	Fitted parameters, the conditional volatility path `sigma` (shape `(n,)`), the standardised residuals `resid` (shape `(n,)`) and forecasting methods.

Raises:

Type	Description
`ValueError`	If `x` has fewer than 50 observations, contains NaN or infinite values, is constant, or if `dist`, `vol` or `mean` is not one of the allowed strings.

Notes

The series is internally rescaled to unit variance before optimising and the parameters are mapped back afterwards. GARCH is exactly scale-equivariant -- :math:x \mapsto cx sends :math:(\mu,\omega) \mapsto (c\mu, c^2\omega) and leaves :math:\alpha,\beta alone -- so this changes nothing statistically, but it keeps every quantity at order 1. Without it, daily returns give :math:\omega \approx 10^{-6}, which sits below the optimiser's convergence tolerance and produces silently unconverged fits.

The variance recursion starts from the sample variance of x and an AR or ARMA mean from its long-run level (the first shock is :math:r_0 - \mu/(1-\phi)), so every model's likelihood sums over all n observations and nested models -- "constant" inside "ar1", "garch" inside "gjr" -- can be compared by AIC or a likelihood-ratio test.

Examples:

Parameters are recovered from a simulated series:

>>> import numpy as np
>>> from rcopula.garch import fit_garch
>>> rng = np.random.default_rng(0)
>>> n, omega, alpha, beta = 8000, 0.05, 0.10, 0.85
>>> s2, e = 1.0, 0.0
>>> x = np.empty(n)
>>> z = rng.standard_normal(n)
>>> for i in range(n):
...     s2 = omega + alpha * e**2 + beta * s2
...     e = np.sqrt(s2) * z[i]
...     x[i] = e
>>> res = fit_garch(x)
>>> bool(abs(res.alpha - 0.10) < 0.03 and abs(res.beta - 0.85) < 0.05)
True

Volatility clustering is removed by the filter -- the whole point:

>>> raw = np.corrcoef(x[1:] ** 2, x[:-1] ** 2)[0, 1]
>>> filtered = np.corrcoef(res.resid[1:] ** 2, res.resid[:-1] ** 2)[0, 1]
>>> bool(filtered < 0.25 * raw)
True

A leverage effect -- negative shocks raising volatility more -- is picked up by vol="gjr", which then beats plain GARCH on AIC:

>>> s2, e = 1.0, 0.0
>>> for i in range(n):
...     s2 = 0.05 + (0.03 + 0.12 * (e < 0)) * e**2 + 0.86 * s2
...     e = np.sqrt(s2) * z[i]
...     x[i] = e
>>> gjr = fit_garch(x, vol="gjr")
>>> bool(abs(gjr.gamma - 0.12) < 0.05 and abs(gjr.beta - 0.86) < 0.05)
True
>>> bool(gjr.aic < fit_garch(x).aic)
True

rcopula.factor ¶

Factor copulas: a few common factors drive hundreds of variables.

An unstructured Gaussian or Student-t copula on d variables has d(d-1)/2 correlations -- 319,600 for 800 stocks, far more than a few years of data can pin down. A factor copula assumes instead that a handful of hidden factors drive all the co-movement. Each variable's latent score is

.. math::

Z_i = a_i M + b_i G_{g(i)} + \sqrt{1 - a_i^2 - b_i^2}\,\varepsilon_i ,

where :math:M is one common ("market") factor, :math:G_g is the factor of variable :math:i's group :math:g(i) (a sector, a region, a rating bucket), and the :math:\varepsilon_i are independent noise. All of them are standard normal, so :math:Z is multivariate normal with correlation matrix

.. math::

R = A A^\top + \operatorname{diag}(1 - a_i^2 - b_i^2),

where row :math:i of :math:A is :math:a_i in the market column and :math:b_i in its group's column. That is 2d parameters instead of d(d-1)/2. The Gaussian factor copula is the copula of :math:Z. The Student-t factor copula divides every :math:Z_i by one shared :math:\sqrt{W/\nu} with :math:W \sim \chi^2_\nu -- a common "panic" variable that makes every pair crash (and rally) together -- which adds one degrees-of-freedom parameter and gives every pair tail dependence.

Both are exactly a Gaussian or Student-t copula with the structured matrix :math:R, so their density is known in closed form. The point of the structure is speed: :math:R^{-1} and :math:\det R follow from the Woodbury identity and the matrix-determinant lemma using only a :math:k \times k matrix (:math:k = number of factors), so the density costs :math:O(dk) per observation and no :math:d \times d matrix is ever inverted. 800 variables are as easy as 8.

============================ ================================================ :class:FactorCopula The model: density, sampler, distribution function, pairwise dependence measures. :func:fit_factor Estimate the loadings (and df) from data. ============================ ================================================

References

Krupskii, P. and Joe, H. (2013). Factor copula models for multivariate data. Journal of Multivariate Analysis 120, 85-101. Hull, J. and White, A. (2004). Valuation of a CDO and an n-th to default CDS without Monte Carlo simulation. Journal of Derivatives 12(2), 8-23. The conditional-independence quadrature used by :meth:FactorCopula.cdf. Lindskog, F., McNeil, A. and Schmock, U. (2003). Kendall's tau for elliptical distributions. In Credit Risk, Physica-Verlag, 149-156. The identity :math:\tau = (2/\pi)\arcsin r, used both ways here. Demarta, S. and McNeil, A. J. (2005). The t copula and related copulas. International Statistical Review 73(1), 111-129. Tail dependence of the t copula and the profile likelihood for df.

FactorCopula ¶

FactorCopula(market: ArrayLike, group_loadings: ArrayLike | None = None, groups: ArrayLike | None = None, *, family: str | None = None, df: float | None = None, free: ArrayLike | None = None)

Bases: Copula

A copula in which one market factor (plus optional group factors) drives every variable.

Each variable i has a hidden score Z_i = a_i * M + b_i * G[g(i)] + sqrt(1 - a_i**2 - b_i**2) * e_i: a loading a_i on a market factor M shared by all variables, a loading b_i on the factor of its group g(i) (for example its sector), and noise of its own. Given the factors, the variables are independent. That needs 2d numbers instead of the d(d-1)/2 correlations of a full matrix, which is what makes 800 stocks tractable. Use :func:fit_factor to estimate one from data.

With family="student" all scores are divided by one shared random scale (a chi-square "panic" variable with df degrees of freedom), so every pair of variables tends to crash -- and rally -- together.

Parameters:

Name	Type	Description	Default
`market`	`array_like of float, shape (d,)`	Market loadings `a_i`, each in `[-1, 1]`. Must have `d >= 2` entries. `nan` means "to be estimated" (see :meth:`unfitted`).	required
`group_loadings`	`array_like of float, shape (d,), or None`	Group loadings `b_i`, each in `[-1, 1]`. Give it together with `groups`, or leave both as `None` for a one-factor model.	`None`
`groups`	`array_like, shape (d,), or None`	Group label of each variable (integers or strings). Labels are re-coded to `0, 1, ..., n_groups - 1` in sorted order; read the codes back from :attr:`groups`.	`None`
`family`	`('gaussian', 'student')`	Keyword-only. `None` means `"student"` if `df` is given and `"gaussian"` otherwise.	`"gaussian"`
`df`	`float or None`	Keyword-only. Degrees of freedom of the Student-t family, at least `0.01`; `None` gives `4.0` when `family="student"`. Not allowed for the Gaussian family.	`None`
`free`	`array_like of bool, shape (n_params,), or None`	Keyword-only. Free/fixed mask, as for every :class:`~rcopula.core.base.Copula`. :func:`rcopula.fit` honours a pinned `df`; loadings cannot be pinned.	`None`

Attributes:

Name	Type	Description
`family`	`str`	`"gaussian"` or `"student"`.
`market`	`numpy.ndarray of float, shape (d,)`	Market loadings (read-only view of :attr:`params`).
`group_loadings`	`numpy.ndarray of float, shape (d,), or None`	Group loadings, or `None` for a one-factor model.
`groups`	`numpy.ndarray of int, shape (d,), or None`	Group code of each variable, in `0..n_groups - 1`.
`group_labels`	`numpy.ndarray, shape (n_groups,), or None`	The original label behind each code.
`n_groups`	`int`	Number of group factors (`0` for a one-factor model).
`n_factors`	`int`	`1 + n_groups`.
`df`	`float`	Degrees of freedom; `inf` for the Gaussian family.
`idiosyncratic`	`numpy.ndarray of float, shape (d,)`	Each variable's own-noise variance, `1 - a_i2 - b_i2`.
`params`	`numpy.ndarray of float, shape (n_params,)`	`market`, then `group_loadings` (if any), then `df` (Student-t only). :attr:`n_params` therefore counts every loading plus `df`.
`dim`	`int`	Number of variables `d`.

Raises:

Type	Description
`ValueError`	If `d < 2`; if exactly one of `group_loadings` and `groups` is given or either has the wrong length; if `family` is unknown or `df` is given with `family="gaussian"`; if a loading is outside `[-1, 1]`; or if some `a_i2 + b_i2 >= 1` (no noise left).

Notes

The copula is exactly a :class:~rcopula.GaussianCopula or :class:~rcopula.StudentCopula whose correlation matrix is R = A A' + diag(1 - a**2 - b**2) (:meth:sigma); :meth:to_elliptical builds that dense equivalent. The density uses the Woodbury identity and the matrix-determinant lemma, so it costs O(d k) per observation for k = 1 + n_groups factors and never forms a d x d inverse: 1,000 days of 800 stocks evaluate in a fraction of a second. The sampler is O(d) per draw.

Pairwise measures follow from R: Kendall's tau is (2/pi) arcsin(r_ij) for both families, Spearman's rho is (6/pi) arcsin(r_ij/2) for the Gaussian (and computed numerically for the Student-t), and the Student-t tail dependence is 2 t_{df+1}(-sqrt((df+1)(1-r_ij)/(1+r_ij))). :meth:tau, :meth:rho and :meth:lambda_ follow the elliptical classes' conventions; the *_matrix methods return full d x d matrices.

Examples:

Six stocks in two sectors, with a Student-t panic factor:

>>> import numpy as np
>>> import rcopula as rc
>>> cop = rc.FactorCopula(
...     [0.6, 0.5, 0.55, 0.6, 0.5, 0.45],
...     group_loadings=[0.4, 0.3, 0.35, 0.3, 0.4, 0.45],
...     groups=[0, 0, 0, 1, 1, 1],
...     df=4.0,
... )
>>> cop.dim, cop.n_factors, cop.n_params
(6, 3, 13)
>>> R = cop.sigma()
>>> round(float(R[0, 1]), 4), round(float(R[0, 3]), 4)   # same vs different sector
(0.42, 0.36)
>>> u = cop.rvs(1000, random_state=0)
>>> u.shape
(1000, 6)

It is exactly the Student-t copula with that matrix:

>>> dense = cop.to_elliptical()
>>> bool(np.allclose(cop.logpdf(u[:5]), dense.logpdf(u[:5]), atol=1e-9))
True
>>> bool(cop.lambda_matrix()[0, 1] > 0)               # crashes are shared
True

family `property` ¶

family: str

The copula family: "gaussian" or "student".

Returns:

Type	Description
`str`	Family name.

market `property` ¶

market: NDArray[float64]

Market loadings a_i, one per variable.

Returns:

Type	Description
`numpy.ndarray of float, shape (d,)`	Read-only view of the first `d` parameters.

group_loadings `property` ¶

group_loadings: NDArray[float64] | None

Group loadings b_i, one per variable, or None without groups.

Returns:

Type	Description
`numpy.ndarray of float, shape (d,), or None`	Read-only view of parameters `d .. 2d - 1`.

groups `property` ¶

groups: NDArray[intp] | None

Group code of each variable, in 0 .. n_groups - 1, or None.

Returns:

Type	Description
`numpy.ndarray of int, shape (d,), or None`	A copy of the codes; :attr:`group_labels` maps codes to labels.

group_labels `property` ¶

group_labels: NDArray[Any] | None

The original group label behind each code, or None.

Returns:

Type	Description
`numpy.ndarray, shape (n_groups,), or None`	Sorted distinct labels; label `group_labels[c]` has code `c`.

n_groups `property` ¶

n_groups: int

Number of group factors; 0 for a one-factor model.

Returns:

Type	Description
`int`	Count of distinct group labels.

n_factors `property` ¶

n_factors: int

Total number of common factors: the market plus one per group.

Returns:

Type	Description
`int`	`1 + n_groups`.

df `property` ¶

df: float

Degrees of freedom of the Student-t family; inf for the Gaussian.

Returns:

Type	Description
`float`	`df` (`nan` while unfitted), or `inf`.

idiosyncratic `property` ¶

idiosyncratic: NDArray[float64]

Each variable's own-noise variance, 1 - a_i**2 - b_i**2.

Returns:

Type	Description
`numpy.ndarray of float, shape (d,)`	Values in `(0, 1]`; small means the variable is almost fully explained by the factors.

param_bounds `property` ¶

param_bounds: list[tuple[float, float]]

Allowed range of each parameter, as (lower, upper) pairs.

Returns:

Type	Description
`list of tuple of (float, float), length n_params`	`(-1, 1)` for every loading and `(0.01, inf)` for `df`. The joint condition `a_i2 + b_i2 < 1` is checked separately.

unfitted `classmethod` ¶

unfitted(dim: int, groups: ArrayLike | None = None, family: str = 'student') -> FactorCopula

Create a factor copula whose loadings are still to be estimated.

Every loading (and df) is nan. Pass it to :func:rcopula.fit or :meth:rcopula.garch.CopulaGarch.fit, which estimate it with :func:fit_factor; it cannot be evaluated before.

Parameters:

Name	Type	Description	Default
`dim`	`int`	Number of variables `d`, at least 2.	required
`groups`	`array_like, shape (d,), or None`	Group label of each variable; `None` for one market factor only.	`None`
`family`	`('gaussian', 'student')`	Copula family.	`"gaussian"`

Returns:

Type	Description
`FactorCopula`	A template with `nan` parameters.

Raises:

Type	Description
`ValueError`	If `dim < 2`, `groups` has the wrong length or `family` is unknown.

Examples:

>>> import rcopula as rc
>>> template = rc.FactorCopula.unfitted(4, groups=[0, 0, 1, 1])
>>> template.n_params
9

loglik ¶

loglik(data: ArrayLike) -> float

Sum the log-density over all rows: how well the model explains the data.

Higher is better. Useful for comparing fitted models on the same data (for example a Gaussian against a Student-t factor copula on held-out days), or for building an AIC/BIC with :attr:n_params.

Parameters:

Name	Type	Description	Default
`data`	`array_like of float, shape (n, d)`	Observations. If every value already lies strictly inside `(0, 1)` they are used as they are; otherwise each column is converted to pseudo-observations (ranks scaled into `(0, 1)`).	required

Returns:

Type	Description
`float`	The total log-likelihood, `sum(log c(u_t))`.

Raises:

Type	Description
`ValueError`	If the copula is unfitted or `data` does not have `d` columns.

Examples:

>>> import rcopula as rc
>>> cop = rc.FactorCopula([0.7, 0.6, 0.5])
>>> u = cop.rvs(500, random_state=1)
>>> bool(cop.loglik(u) > 0)
True

cdf ¶

cdf(u: ArrayLike) -> NDArray[np.float64]

Probability that every variable is at or below the given point, C(u).

Computed by conditioning on the factors: given the market factor and the group factors the variables are independent, so C(u) is a product of normal probabilities averaged over the factors.

Parameters:

Name	Type	Description	Default
`u`	`array_like of float, shape (n, d) or (d,)`	Points, one row per point; values outside `[0, 1]` are clipped.	required

Returns:

Type	Description
`numpy.ndarray of float, shape (n,)`	Probabilities in `[0, 1]`.

Raises:

Type	Description
`ValueError`	If the copula is unfitted or `u` does not have `d` columns.

Notes

The integral over the market factor and over each group factor uses 48-node Gauss-Hermite quadrature, nested (group factors are independent given the market, so the cost is O(48^2 d) per point, not exponential in the number of groups). The Student-t family adds a 40-node generalised Gauss-Laguerre rule over the chi-square scale, multiplying the cost by 40 (below 8 degrees of freedom, a 256-node Gauss-Legendre rule on the chi variable's probability scale, which handles the heavy tails). Everything is summed in logs, so tiny joint probabilities do not underflow. Agreement with the dense :class:~rcopula.GaussianCopula / :class:~rcopula.StudentCopula distribution functions is typically better than 1e-4.

Examples:

>>> import rcopula as rc
>>> cop = rc.FactorCopula([0.6, 0.6])
>>> dense = cop.to_elliptical()          # Gaussian copula, r = 0.36
>>> bool(abs(cop.cdf([0.3, 0.4])[0] - dense.cdf([0.3, 0.4])[0]) < 1e-4)
True

sigma ¶

sigma() -> NDArray[np.float64]

Return the implied d x d correlation matrix R = A A' + diag(noise).

Entry [i, j] is a_i a_j plus b_i b_j when i and j share a group; the diagonal is 1. This is the correlation matrix of the latent scores (before the Student-t scaling), matching :meth:rcopula.GaussianCopula.sigma.

Returns:

Type	Description
`numpy.ndarray of float, shape (d, d)`	Symmetric positive-definite correlation matrix.

Raises:

Type	Description
`ValueError`	If the copula is unfitted.

Examples:

>>> import rcopula as rc
>>> rc.FactorCopula([0.8, 0.5]).sigma()
array([[1. , 0.4],
       [0.4, 1. ]])

to_elliptical ¶

to_elliptical() -> GaussianCopula | StudentCopula

Build the equivalent dense Gaussian or Student-t copula.

The factor copula is an elliptical copula with correlation matrix :meth:sigma; this returns that copula with an unstructured (dispstr="un") matrix. Handy for checks and for functions that need an elliptical copula, but it has d(d-1)/2 parameters and a O(d^3) density, so prefer the factor form for large d.

Returns:

Type	Description
`GaussianCopula or StudentCopula`	`StudentCopula` (with the same `df`) for the Student-t family.

Raises:

Type	Description
`ValueError`	If the copula is unfitted.

Examples:

>>> import rcopula as rc
>>> rc.FactorCopula([0.8, 0.5], df=5.0).to_elliptical().describe()
'Student copula, dim 2, rho.12=0.4, df=5'

tau_matrix ¶

tau_matrix() -> NDArray[np.float64]

Kendall's tau of every pair, as a d x d matrix.

Uses tau_ij = (2/pi) arcsin(r_ij), exact for both families.

Returns:

Type	Description
`numpy.ndarray of float, shape (d, d)`	Symmetric, ones on the diagonal.

Raises:

Type	Description
`ValueError`	If the copula is unfitted.

Examples:

>>> import rcopula as rc
>>> float(round(rc.FactorCopula([0.8, 0.5]).tau_matrix()[0, 1], 6))
0.26198

rho_matrix ¶

rho_matrix() -> NDArray[np.float64]

Spearman's rho of every pair, as a d x d matrix.

Returns:

Type	Description
`numpy.ndarray of float, shape (d, d)`	Symmetric, ones on the diagonal.

Raises:

Type	Description
`ValueError`	If the copula is unfitted.

Notes

Gaussian: (6/pi) arcsin(r_ij / 2), exact. Student-t: there is no closed form, so each distinct correlation is integrated numerically as in :meth:rcopula.StudentCopula.rho; with more than 15 distinct values, a polynomial through 15 integrated Chebyshev points spanning them is used instead (error about 1e-8), which keeps 800 variables to several seconds rather than hours.

Examples:

>>> import rcopula as rc
>>> float(round(rc.FactorCopula([0.8, 0.5]).rho_matrix()[0, 1], 6))
0.384565

lambda_matrix ¶

lambda_matrix() -> NDArray[np.float64]

Tail dependence of every pair, as a d x d matrix (same in both tails).

Returns:

Type	Description
`numpy.ndarray of float, shape (d, d)`	Zero off the diagonal for the Gaussian family; `2 t_{df+1}(-sqrt((df+1)(1-r)/(1+r)))` for the Student-t. Ones on the diagonal.

Raises:

Type	Description
`ValueError`	If the copula is unfitted.

Examples:

>>> import rcopula as rc
>>> cop = rc.FactorCopula([0.8, 0.5], df=4.0)
>>> lam = cop.lambda_matrix()[0, 1]
>>> bool(abs(lam - cop.to_elliptical().lambda_().upper) < 1e-12)
True

tau ¶

tau() -> Any

Kendall's tau implied by the copula: one value per pair.

Returns:

Type	Description
`float or numpy.ndarray of float, shape (d * (d - 1) / 2,)`	A float when `d == 2`; otherwise every pair in :func:`~rcopula.P2p` order (as :meth:`rcopula.GaussianCopula.tau` does). :meth:`tau_matrix` gives the same values as a matrix.

Raises:

Type	Description
`ValueError`	If the copula is unfitted.

Examples:

>>> import rcopula as rc
>>> round(rc.FactorCopula([0.8, 0.5]).tau(), 6)
0.26198

rho ¶

rho() -> Any

Spearman's rho implied by the copula: one value per pair.

Returns:

Type	Description
`float or numpy.ndarray of float, shape (d * (d - 1) / 2,)`	A float when `d == 2`; otherwise every pair in :func:`~rcopula.P2p` order. See :meth:`rho_matrix`.

Raises:

Type	Description
`ValueError`	If the copula is unfitted.

lambda_ ¶

lambda_() -> TailDependence

Tail dependence, when every pair shares one value.

Returns:

Type	Description
`TailDependence`	`(lower, upper)`, equal to each other; zero for the Gaussian.

Raises:

Type	Description
`ValueError`	If the copula is unfitted, or if it is a Student-t factor copula whose pairs have different tail dependence (the usual case for `d > 2`); use :meth:`lambda_matrix` then.

Examples:

>>> import rcopula as rc
>>> rc.FactorCopula([0.8, 0.5]).lambda_()
TailDependence(lower=0.0, upper=0.0)

describe ¶

describe() -> str

Short human-readable summary: family, size, factors and loadings.

Returns:

Type	Description
`str`	One line, e.g. `"Factor copula (Student-t, df=4), dim 800, 1 market + 10 group factors, market loadings 0.45..0.65 (mean 0.55), group loadings 0.30..0.50 (mean 0.40), 1601 parameters"`.

Examples:

>>> import rcopula as rc
>>> rc.FactorCopula([0.6, 0.5, 0.4]).describe()
'Factor copula (Gaussian), dim 3, 1 market factor,
 market loadings 0.4..0.6 (mean 0.5), 3 parameters'

fit_factor ¶

fit_factor(data: ArrayLike, groups: ArrayLike | None = None, family: str = 'student', *, method: str = 'itau', df: float | None = None, df_grid: Sequence[float] = DEFAULT_DF_GRID, max_iter: int = 2000, tol: float = 1e-07, max_communality: float = 0.98) -> FactorCopula

Fit a factor copula to data: loadings by correlation matching, df by likelihood.

Give it a table of observations (one column per variable) and, optionally, each variable's group (for example its sector). It returns a :class:FactorCopula with one market factor plus one factor per group. Built for many variables: 800 columns of 1,000 rows take a few seconds.

Parameters:

Name	Type	Description	Default
`data`	`array_like or pandas.DataFrame of float, shape (n, d)`	Observations, one row per observation and one column per variable, with `d >= 2`. If every value lies strictly inside `(0, 1)` the data are taken as pseudo-observations; otherwise each column is rank-transformed first.	required
`groups`	`array_like, shape (d,), or None`	Group label of each column. `None` fits one market factor only.	`None`
`family`	`('student', 'gaussian')`	Copula family. The Student-t adds a degrees-of-freedom parameter and tail dependence.	`"student"`
`method`	`('itau', 'normal')`	Which correlation matrix the loadings are matched to. `"itau"` inverts Kendall's tau, `r = sin(pi * tau / 2)`, which is consistent for both families (heavy joint tails do not bias it). `"normal"` uses the Pearson correlation of normal scores, which is faster but consistent only for the Gaussian family.	`"itau"`
`df`	`float or None`	Keyword-only. Hold the Student-t degrees of freedom at this value instead of estimating them. Not allowed with `family="gaussian"`.	`None`
`df_grid`	`sequence of float`	Keyword-only. Degrees of freedom at which the profile log-likelihood is evaluated before the best one is refined. The estimate is confined to `[min(df_grid), max(df_grid)]`.	`DEFAULT_DF_GRID`
`max_iter`	`int`	Keyword-only. Maximum number of correlation-matching sweeps.	`2000`
`tol`	`float`	Keyword-only. The sweeps stop when no loading moves by more than this.	`1e-7`
`max_communality`	`float`	Keyword-only. Upper limit on `a_i2 + b_i2`, so every variable keeps at least `1 - max_communality` of noise of its own (the density needs it to be positive). Must lie in `(0, 1)`.	`0.98`

Returns:

Type	Description
`FactorCopula`	The fitted copula. Its :meth:`~FactorCopula.loglik` on the fitting data is the maximised (profile) log-likelihood.

Raises:

Type	Description
`ValueError`	If `data` is not 2-D with at least two columns and two rows, if `groups` has the wrong length, if `family` or `method` is unknown, if `df` is given for the Gaussian family, if `max_communality` is not in `(0, 1)`, or if the correlation matrix has missing values (a constant column, for instance).

Notes

Loadings. The model implies r_ij = a_i a_j + b_i b_j [same group] off the diagonal. The loadings minimise the squared distance between that and the target correlation matrix over all off-diagonal pairs, by alternating least squares: each sweep updates every market loading given the rest, then every group loading given the new market loadings (the iterated principal-factor idea, restricted to the factor pattern). Because each loading is fitted to d - 1 correlations, the noise in the individual correlations largely averages out. Group loadings are kept non-negative (a group factor adds correlation within its group); a group with a single member gets loading 0, since it cannot be identified.

Degrees of freedom. With the loadings fixed, the Student-t log-likelihood is evaluated at every df_grid value and the best one is refined by a bounded one-dimensional search on log(df) between its grid neighbours. Each evaluation is O(n d k).

Accuracy. Simulated from a Student-t factor copula with df = 4, 50 variables in 5 groups, market loadings in [0.4, 0.7] and group loadings in [0.3, 0.5], 1,000 rows: mean absolute loading error about 0.02 (market) and 0.03 (group), and df within about 0.5 of the truth. At 800 variables in 10 groups the loading errors shrink to about 0.015 and the whole fit takes 5-10 seconds, most of it the 319,600 Kendall's taus (:func:~rcopula.cor_kendall).

This is a two-stage, moment-based estimator, not full maximum likelihood, so there are no standard errors; bootstrap the rows if you need them.

Examples:

>>> import numpy as np
>>> import rcopula as rc
>>> rng = np.random.default_rng(0)
>>> groups = np.repeat([0, 1, 2], 6)
>>> truth = rc.FactorCopula(
...     rng.uniform(0.4, 0.7, 18), rng.uniform(0.3, 0.5, 18), groups, df=4.0
... )
>>> u = truth.rvs(2000, random_state=1)
>>> fitted = rc.fit_factor(u, groups=groups)
>>> bool(np.mean(np.abs(fitted.market - truth.market)) < 0.05)
True
>>> bool(2.5 < fitted.df < 6.5)
True
>>> fitted.n_params
37

rcopula.dynamic ¶

Time-varying copulas: dependence that moves.

Every copula elsewhere in this package is constant -- one parameter for the whole sample. That is often indefensible. Dependence between assets rises in a crisis and falls afterwards; dependence between a river's peak and its volume changes with land use; dependence between two sensors changes with the weather. Fitting one number to all of it produces an average that describes no particular period.

The fix is to let the copula parameter follow an observation-driven recursion, so :math:\theta_t is a deterministic function of the past. The likelihood stays exact -- there is no filtering approximation and no latent state to integrate out -- which is the whole appeal of the observation-driven class.

Two recursions are provided.

Patton (2006). The parameter responds to a moving average of a forcing term computed from the last :math:q observations:

.. math::

\theta_t = \Lambda\!\left(\omega + \beta\,\theta_{t-1}
           + \alpha\, \frac{1}{q}\sum_{j=1}^{q} m(u_{t-j}, v_{t-j})\right),

with :math:\Lambda a link keeping :math:\theta_t inside its domain. Patton used :math:m = \Phi^{-1}(u)\Phi^{-1}(v) for correlation parameters and :math:m = |u - v| otherwise. The second is decreasing in dependence, so a fitted :math:\alpha is normally negative for it -- a sign flip that is easy to misread as a bug.

GAS (Creal, Koopman and Lucas, 2013). The parameter responds to the score of its own likelihood, which is the direction that would most improve the fit at the last observation:

.. math::

f_{t+1} = \omega + \beta f_t + \alpha s_t, \qquad
s_t = \frac{\partial \log c(u_t, v_t; \Lambda(f_t))}{\partial f_t},
\qquad \theta_t = \Lambda(f_t).

The score is family-specific and automatic: no forcing term has to be chosen, and a tail observation moves a Clayton copula differently from a Gumbel one. This usually fits better than Patton's recursion and is harder to explain, which is a fair summary of the trade-off between them.

For :math:d > 2 elliptical copulas the natural object is a whole correlation matrix, and the standard recursion for that is Engle's DCC -- see :func:fit_dcc.

============================ ================================================ :class:DynamicCopula A family plus a parameter recursion. :func:fit_dynamic Estimate the recursion by maximum likelihood. :class:DynamicFitResult What came back, with a test against constancy. :func:fit_dcc Engle's DCC for a time-varying correlation matrix. :class:DccResult The filtered correlation path. ============================ ================================================

.. warning::

A dynamic copula is fitted to pseudo-observations of i.i.d. innovations, not to raw data. Volatility clustering left in the margins will be picked up by this recursion and reported as time-varying dependence, which is precisely the confusion :mod:rcopula.garch exists to prevent. Filter the margins first.

Examples:

>>> import numpy as np, rcopula as rc
>>> from rcopula.dynamic import DynamicCopula
>>> model = DynamicCopula(rc.GaussianCopula(0.0), coefficients=(0.02, 0.1, 0.9))
>>> u = model.simulate(500, random_state=0)
>>> path = model.filter(u).path
>>> float(path.min()) > -1.0 and float(path.max()) < 1.0
True

References

Patton, A. J. (2006). Modelling asymmetric exchange rate dependence. International Economic Review 47(2), 527-556. Creal, D., Koopman, S. J. and Lucas, A. (2013). Generalized autoregressive score models with applications. Journal of Applied Econometrics 28(5), 777-795. Engle, R. F. (2002). Dynamic conditional correlation: a simple class of multivariate generalized autoregressive conditional heteroskedasticity models. Journal of Business and Economic Statistics 20(3), 339-350. Aielli, G. P. (2013). Dynamic conditional correlation: on properties and estimation. Journal of Business and Economic Statistics 31(3), 282-299. Why the two-step DCC estimator is inconsistent, and the cDCC correction. Manner, H. and Reznikova, O. (2012). A survey on time-varying copulas: specification, simulations, and application. Econometric Reviews 31(6), 654-687.

DynamicPath `dataclass` ¶

DynamicPath(path: NDArray[float64], linked: NDArray[float64], loglik_contributions: NDArray[float64])

The copula parameter at every time step, as worked out by a dynamic copula.

This is what :meth:DynamicCopula.filter returns: one parameter value per observation, plus how well each observation fits the copula in force at that moment. You do not normally build one yourself.

Parameters:

Name	Type	Description	Default
`path`	`numpy.ndarray of float, shape (n,)`	See `Attributes`.	required
`linked`	`numpy.ndarray of float, shape (n,)`	See `Attributes`.	required
`loglik_contributions`	`numpy.ndarray of float, shape (n,)`	See `Attributes`.	required

Attributes:

Name	Type	Description
`path`	`numpy.ndarray of float, shape (n,)`	The copula parameter at each observation, on its natural scale (for a Gaussian copula, the correlation at each time step).
`linked`	`numpy.ndarray of float, shape (n,)`	The same path before the link -- the unconstrained scale the recursion runs on.
`loglik_contributions`	`numpy.ndarray of float, shape (n,)`	`log c(u_t; theta_t)` term by term. Summing gives the log-likelihood; the series itself shows which observations the model finds surprising. A non-finite term is replaced by `-1e10`.
`loglik`	`float`	Total log-likelihood, the sum of `loglik_contributions` (a read-only property).

loglik `property` ¶

loglik: float

Total log-likelihood: how well the whole path fits the data (higher is better).

Returns:

Type	Description
`float`	The sum of `loglik_contributions`.

DynamicCopula ¶

DynamicCopula(family: Copula, *, coefficients: ArrayLike, driver: Driver = 'patton', forcing: Forcing = 'auto', lags: int = 10, index: int = 0, bounds: tuple[float, float] | None = None)

A two-variable copula whose strength of dependence changes over time.

An ordinary copula has one fixed parameter for the whole sample. Here the parameter is updated after every observation by a simple rule (an observation-driven recursion), so it can rise in turbulent periods and fall in calm ones. Use it to evaluate, simulate or forecast a model whose coefficients you already know; to estimate the coefficients from data, use :func:fit_dynamic.

Parameters:

Name	Type	Description	Default
`family`	`Copula`	A bivariate copula (`family.dim == 2`), for example `rcopula.GaussianCopula(0.3)`. Its current parameter values set the starting point of the recursion and, for a multi-parameter family such as :class:`~rcopula.StudentCopula`, fix everything that does not vary (the degrees of freedom stay constant).	required
`coefficients`	`array_like of float, shape (3,)`	`(omega, alpha, beta)`, required and keyword-only. `omega` is the intercept, `alpha` how strongly the parameter reacts to new data, and `beta` how persistent it is (values near 1 change slowly). They act on the linked (unconstrained) scale, so any real numbers are valid.	required
`driver`	`('patton', 'gas')`	Which recursion. `"patton"` uses a moving average of a forcing term; `"gas"` uses the likelihood score. See the module docstring.	`"patton"`
`forcing`	`('auto', 'normal-product', 'abs-difference')`	Patton's forcing term. `"auto"` picks the normal product for elliptical families and the absolute difference otherwise, which is what Patton did. The GAS driver needs no such choice: passing anything other than `"auto"` with `driver="gas"` has no effect and issues a `UserWarning`.	`"auto"`
`lags`	`int`	Length of the moving-average window, in observations; must be at least 1. Patton used 10. Only used by the `"patton"` driver.	`10`
`index`	`int`	Which entry of `family.params` varies, when the family has more than one. The default is the first, which is the correlation for a t copula.	`0`
`bounds`	`tuple of (float, float) or None`	`(lower, upper)` working range for the varying parameter, with `lower < upper`. `None` uses the family's own domain, with an infinite side replaced by a finite span (25 for a one-sided parameter such as Clayton's, `(-40, 40)` for Frank's).	`None`

Attributes:

Name	Type	Description
`family`	`Copula`	The copula passed in.
`coefficients`	`numpy.ndarray of float, shape (3,)`	`(omega, alpha, beta)`.
`driver`	`str`	`"patton"` or `"gas"`.
`forcing`	`str`	The resolved forcing term, `"normal-product"` or `"abs-difference"` (never `"auto"`). Unused by the GAS driver.
`lags`	`int`	Moving-average window length.
`index`	`int`	Position of the varying parameter in `family.params`.
`link`	`object`	The map from the unconstrained scale onto `(lower, upper)`; its `lower` and `upper` attributes give the working range.

Raises:

Type	Description
`ValueError`	If `family` is not bivariate, `coefficients` does not have exactly three entries, `driver` is not `"patton"` or `"gas"`, `forcing` is not one of the listed names, `lags` is not a whole number of at least 1, or `bounds` is not increasing.

Warns:

Type	Description
`UserWarning`	If `forcing` is set explicitly with `driver="gas"`, where it is ignored.

Notes

The likelihood is evaluated observation by observation, because the parameter differs at every one. Expect roughly 50 microseconds per observation per likelihood evaluation -- a few seconds to fit a decade of daily data, not milliseconds.

Examples:

>>> import rcopula as rc
>>> from rcopula.dynamic import DynamicCopula
>>> model = DynamicCopula(rc.ClaytonCopula(1.0), coefficients=(0.1, -0.5, 0.9))
>>> model.driver
'patton'
>>> model.forcing
'abs-difference'

filter ¶

filter(u: ArrayLike) -> DynamicPath

Work out the copula parameter at every time step of u.

Runs the recursion forwards through the data with the model's fixed coefficients, and reports the parameter path and how well each observation fits.

Parameters:

Name	Type	Description	Default
`u`	`array_like of float, shape (n, 2)`	Pseudo-observations, each value strictly between 0 and 1, in time order (oldest first). Unlike :func:`fit_dynamic`, raw data are not converted here.	required

Returns:

Type	Description
`DynamicPath`	`path` (parameter per observation, shape `(n,)`), `linked` (the same on the unconstrained scale) and `loglik_contributions` (per-observation log-density); `.loglik` gives the total.

Raises:

Type	Description
`ValueError`	If `u` is not two-dimensional with exactly two columns.

Examples:

>>> import numpy as np, rcopula as rc
>>> from rcopula.dynamic import DynamicCopula
>>> model = DynamicCopula(rc.GaussianCopula(0.5), coefficients=(0.0, 0.0, 0.0))
>>> u = rc.GaussianCopula(0.5).rvs(50, random_state=0)
>>> path = model.filter(u).path
>>> bool(np.allclose(path, 0.0))     # omega = alpha = beta = 0 pins it at Lambda(0)
True

loglik ¶

loglik(u: ArrayLike) -> float

Score how well the model fits u: the log-likelihood (higher is better).

Shorthand for self.filter(u).loglik.

Parameters:

Name	Type	Description	Default
`u`	`array_like of float, shape (n, 2)`	Pseudo-observations strictly between 0 and 1, in time order.	required

Returns:

Type	Description
`float`	Sum over observations of `log c(u_t; theta_t)`.

Raises:

Type	Description
`ValueError`	If `u` does not have shape `(n, 2)`.

with_coefficients ¶

with_coefficients(coefficients: ArrayLike) -> DynamicCopula

Make a copy of this model with different (omega, alpha, beta).

Everything else -- family, driver, forcing, lags, varying parameter and working range -- is carried over unchanged. The original is not modified.

Parameters:

Name	Type	Description	Default
`coefficients`	`array_like of float, shape (3,)`	The new `(omega, alpha, beta)`.	required

Returns:

Type	Description
`DynamicCopula`	A new model.

Raises:

Type	Description
`ValueError`	If `coefficients` does not have exactly three entries.

Examples:

>>> import rcopula as rc
>>> from rcopula.dynamic import DynamicCopula
>>> model = DynamicCopula(rc.GaussianCopula(0.0), coefficients=(0.0, 0.1, 0.9))
>>> model.with_coefficients((0.0, 0.2, 0.5)).coefficients.tolist()
[0.0, 0.2, 0.5]

simulate ¶

simulate(size: int, *, random_state: Any = None, burn_in: int = 200) -> NDArray[np.float64]

Generate a random time series of pseudo-observations from the model.

Useful for checking a fit (simulate from known coefficients, refit, compare) or for scenario analysis. The recursion feeds on its own output, so this is genuinely sequential: each observation is drawn from the copula the previous ones implied.

Parameters:

Name	Type	Description	Default
`size`	`int`	Number of observations to return; a non-negative whole number.	required
`random_state`	`(int, Generator or None)`	Seed or generator for reproducible draws. `None` uses fresh randomness.	`None`
`burn_in`	`int`	Draws generated and thrown away first, so the returned path does not depend on where `family.params` happened to start. Non-negative.	`200`

Returns:

Type	Description
`numpy.ndarray of float, shape (size, 2)`	Pseudo-observations in time order, each value between 0 and 1.

Examples:

>>> import numpy as np, rcopula as rc
>>> from rcopula.dynamic import DynamicCopula
>>> model = DynamicCopula(rc.GaussianCopula(0.0), coefficients=(0.05, 0.2, 0.9))
>>> u = model.simulate(200, random_state=1)
>>> u.shape
(200, 2)

forecast ¶

forecast(u: ArrayLike, horizon: int, *, draws: int = 2000, random_state: Any = None) -> dict[str, NDArray[np.float64]]

Predict where the copula parameter is heading over the next few steps.

The model is first run over u, then many possible futures are simulated from its last state. The recursion is driven by future observations, so beyond one step ahead the parameter is a random variable rather than a number; this returns a summary of its distribution at each step.

Step 1 is not random: it is the filter's own one-step-ahead parameter, determined by the observations in u (for the GAS driver this includes the score of the last row), so all four summaries coincide there and equal what :meth:filter would assign to the next observation. Each later step uses the observation simulated at the step before, with its score taken at the state it was drawn from.

Parameters:

Name	Type	Description	Default
`u`	`array_like of float, shape (n, 2)`	Pseudo-observations strictly between 0 and 1, in time order; the forecast starts after the last row.	required
`horizon`	`int`	Number of steps ahead to forecast; at least 1.	required
`draws`	`int`	Number of simulated futures. More gives smoother percentiles at proportionally higher cost.	`2000`
`random_state`	`(int, Generator or None)`	Seed or generator for reproducible simulation.	`None`

Returns:

Type	Description
`dict of str to numpy.ndarray of float`	Four arrays, each of shape `(horizon,)`, with entry `h` for the step `h + 1` ahead: `"mean"`: average simulated parameter. `"median"`: median simulated parameter. `"lower"`: 5th percentile. `"upper"`: 95th percentile.

Raises:

Type	Description
`ValueError`	If `horizon` is less than 1, or `u` does not have shape `(n, 2)`.

Examples:

>>> import rcopula as rc
>>> from rcopula.dynamic import DynamicCopula
>>> model = DynamicCopula(rc.GaussianCopula(0.4), coefficients=(0.05, 0.1, 0.85))
>>> u = model.simulate(300, random_state=0)
>>> ahead = model.forecast(u, horizon=5, draws=200, random_state=1)
>>> sorted(ahead)
['lower', 'mean', 'median', 'upper']
>>> ahead["mean"].shape
(5,)

DynamicFitResult `dataclass` ¶

DynamicFitResult(model: DynamicCopula, coefficients: NDArray[float64], loglik: float, path: NDArray[float64], constant_loglik: float, constant_param: float, n_obs: int, converged: bool, message: str = '', _u: NDArray[float64] = (lambda: np.empty((0, 2)))())

The result of fitting a time-varying copula with :func:fit_dynamic.

Holds the estimated coefficients, the fitted parameter path, and the fit of an ordinary constant-parameter copula for comparison, so you can judge whether letting dependence move was worth it. You do not normally build one yourself.

Every constructor argument is the field of the same name listed under Attributes; message defaults to "".

Attributes:

Name	Type	Description
`model`	`DynamicCopula`	The model at the estimated coefficients; use it to filter new data, simulate or forecast.
`coefficients`	`numpy.ndarray of float, shape (3,)`	Estimated `(omega, alpha, beta)`.
`loglik`	`float`	Maximised log-likelihood of the time-varying model.
`path`	`numpy.ndarray of float, shape (n_obs,)`	The filtered parameter path: the copula parameter at each observation.
`constant_loglik`	`float`	Log-likelihood of the same family with a constant parameter, fitted by maximum likelihood. The comparison is the point of the exercise.
`constant_param`	`float`	The fitted constant value of the varying parameter.
`n_obs`	`int`	Number of observations used.
`converged`	`bool`	Whether the optimiser reported success.
`message`	`str`	The optimiser's status message.
`n_params, aic, bic, persistence`		Read-only properties; see each one below.

n_params `property` ¶

n_params: int

Number of estimated coefficients: three, always (omega, alpha and beta).

Returns:

Type	Description
`int`	`3`.

aic `property` ¶

aic: float

Akaike information criterion: fit penalised for complexity (lower is better).

Returns:

Type	Description
`float`	`2 * n_params - 2 * loglik`.

bic `property` ¶

bic: float

Bayesian information criterion: like AIC with a stiffer penalty (lower is better).

Returns:

Type	Description
`float`	`n_params * log(n_obs) - 2 * loglik`.

persistence `property` ¶

persistence: float

How slowly the dependence changes: the autoregressive coefficient beta.

Near 1 the parameter drifts slowly and the recursion is close to a unit root; near 0 it is essentially the constant model with noise.

Returns:

Type	Description
`float`	`coefficients[2]`.

constancy_test ¶

constancy_test() -> tuple[float, float]

Check whether a time-varying copula fits clearly better than a constant one.

Computes the likelihood-ratio statistic comparing this fit with the constant-parameter fit of the same family. A large statistic (small p-value) suggests the dependence really does move -- but read the caveat in the Notes before relying on the p-value.

Returns:

Name	Type	Description
`statistic`	`float`	`2 * (loglik - constant_loglik)`, floored at 0.
`pvalue`	`float`	Nominal p-value from a chi-squared distribution with 2 degrees of freedom. Approximate only; see Notes.

Notes

The null sits on the boundary of the parameter space in beta (a constant copula is alpha = 0 with beta unidentified), so the chi-squared reference distribution is not exact -- this is the Davies problem, and the reported p-value is conservative-to-wrong in the usual way for a nuisance parameter identified only under the alternative. Treat it as a diagnostic, not a test. A bootstrap under the constant null is the honest version and is left to the caller.

Examples:

>>> import rcopula as rc
>>> from rcopula.dynamic import DynamicCopula, fit_dynamic
>>> model = DynamicCopula(rc.GaussianCopula(0.0), coefficients=(0.1, 0.3, 0.8))
>>> u = model.simulate(600, random_state=0)
>>> res = fit_dynamic(u, rc.GaussianCopula(0.0))
>>> statistic, pvalue = res.constancy_test()
>>> bool(statistic > 0)
True

summary ¶

summary() -> str

Produce a plain-text report of the fit, ready to print.

Lists the coefficients, log-likelihood, AIC/BIC, the constant-copula comparison and the range of the fitted parameter path, plus a warning line if the optimiser did not converge.

Returns:

Type	Description
`str`	A multi-line report.

Examples:

>>> import rcopula as rc
>>> from rcopula.dynamic import DynamicCopula, fit_dynamic
>>> model = DynamicCopula(rc.GaussianCopula(0.0), coefficients=(0.1, 0.3, 0.8))
>>> u = model.simulate(400, random_state=2)
>>> text = fit_dynamic(u, rc.GaussianCopula(0.0)).summary()
>>> print(text.splitlines()[0])
Time-varying Gaussian copula, dim 2, rho=0.151034 (patton recursion)
>>> "LR vs constant" in text
True

DccResult `dataclass` ¶

DccResult(correlations: NDArray[float64], a: float, b: float, loglik: float, unconditional: NDArray[float64], df: float | None, converged: bool, n_obs: int)

The result of :func:fit_dcc: a correlation matrix for every time step.

Holds the estimated DCC coefficients and the whole path of correlation matrices they imply, so you can see how the dependence between every pair of variables moved through the sample. You do not normally build one yourself.

Every constructor argument is the field of the same name listed under Attributes.

Attributes:

Name	Type	Description
`correlations`	`numpy.ndarray of float, shape (n, d, d)`	The correlation matrix at each observation; `correlations[t]` is symmetric with a unit diagonal.
`a`	`float`	The DCC news-impact coefficient: how strongly the latest observation moves the correlations. Non-negative.
`b`	`float`	The DCC persistence coefficient: how much of the previous correlation carries over. Non-negative, with `a + b < 1`.
`loglik`	`float`	Copula log-likelihood, i.e. of the Gaussian (or t) copula density at the filtered correlations -- not the joint likelihood of the data.
`unconditional`	`numpy.ndarray of float, shape (d, d)`	The target correlation matrix the recursion reverts to.
`df`	`float or None`	Degrees of freedom, if a t copula was used; `None` for Gaussian.
`converged`	`bool`	Whether the optimiser reported success.
`n_obs`	`int`	Number of observations, `n`.
`persistence`	`float`	`a + b` (a read-only property; see below).

persistence `property` ¶

persistence: float

How long shocks to correlation last: a + b.

At 1 the recursion has a unit root and never reverts; values close to 1 (0.95-0.99 is common for daily data) mean slow mean reversion.

Returns:

Type	Description
`float`	`a + b`, between 0 and 1.

pair ¶

pair(i: int, j: int) -> NDArray[np.float64]

Get the correlation between variables i and j at every time step.

Parameters:

Name	Type	Description	Default
`i`	`int`	Column index of the first variable, `0 <= i < d`.	required
`j`	`int`	Column index of the second variable, `0 <= j < d`.	required

Returns:

Type	Description
`numpy.ndarray of float, shape (n,)`	`correlations[:, i, j]`.

Raises:

Type	Description
`IndexError`	If `i` or `j` is out of range.

Examples:

>>> import rcopula as rc
>>> from rcopula.dynamic import fit_dcc
>>> u = rc.GaussianCopula(0.5, dim=3, dispstr="ex").rvs(200, random_state=0)
>>> fit_dcc(u).pair(0, 2).shape
(200,)

copulas ¶

copulas() -> list[Copula]

Turn the correlation path into one ordinary copula object per time step.

Each entry is a Gaussian copula (or a t copula with df degrees of freedom) with unstructured correlation equal to correlations[t], so you can sample from or evaluate the model in force at any date. Building n copula objects is not free; do it when the models are wanted individually, not to compute a likelihood.

Returns:

Type	Description
`list of Copula, length n`	:class:`~rcopula.GaussianCopula` objects if `df` is `None`, otherwise :class:`~rcopula.StudentCopula`, each of dimension `d`.

summary ¶

summary() -> str

Produce a plain-text report of the DCC fit, ready to print.

Lists a, b, their sum, the copula log-likelihood, and the mean and range of the correlation between the first two variables.

Returns:

Type	Description
`str`	A multi-line report.

fit_dynamic ¶

fit_dynamic(u: ArrayLike, family: Copula, *, driver: Driver = 'patton', forcing: Forcing = 'auto', lags: int = 10, index: int = 0, bounds: tuple[float, float] | None = None, start: ArrayLike | None = None) -> DynamicFitResult

Estimate, from data, a two-variable copula whose dependence changes over time.

Finds the (omega, alpha, beta) coefficients of a :class:DynamicCopula that make the observed series most likely (maximum likelihood), and also fits an ordinary constant copula of the same family so the two can be compared. Use it when you suspect the link between two series is stronger in some periods than others.

Parameters:

Name	Type	Description	Default
`u`	`array_like of float, shape (n, 2)`	Two series in time order (oldest first). Ideally pseudo-observations strictly between 0 and 1. If any value is 0, 1 or outside that range, the data are converted with :func:`~rcopula.pseudo_obs` first, so raw innovations are accepted -- but see the warning in the module docstring about which series to pass.	required
`family`	`Copula`	The bivariate family whose parameter varies, for example `rcopula.GaussianCopula(0.0)`. Its starting parameter values do not matter; they are replaced by the constant fit.	required
`driver`	`('patton', 'gas')`	Which recursion; as for :class:`DynamicCopula`.	`"patton"`
`forcing`	`('auto', 'normal-product', 'abs-difference')`	Patton's forcing term; as for :class:`DynamicCopula`.	`"auto"`
`lags`	`int`	Moving-average window for the Patton driver; at least 1.	`10`
`index`	`int`	Which entry of `family.params` varies.	`0`
`bounds`	`tuple of (float, float) or None`	Working range of the varying parameter; `None` uses the family's domain, as for :class:`DynamicCopula`.	`None`
`start`	`array_like of float, shape (3,), or None`	Starting `(omega, alpha, beta)` for the optimiser. The default tries three starts, led by a persistent recursion centred on the constant maximum-likelihood estimate, which is the standard choice and matters: this likelihood is not concave.	`None`

Returns:

Type	Description
`DynamicFitResult`	Estimated `coefficients`, fitted `model` and parameter `path`, `loglik`, the constant fit's `constant_loglik` and `constant_param`, `n_obs`, `converged` and `message`. Call `.summary()` for a report or `.constancy_test()` for the comparison.

Raises:

Type	Description
`ValueError`	If `u` does not have shape `(n, 2)`. A `family` that is not bivariate is also rejected (by the constant fit or by :class:`DynamicCopula`).

Notes

Optimisation is Nelder-Mead on the negative log-likelihood. Coefficients with |beta| >= 0.9999 are rejected as non-stationary. Fitting is slow for families without a vectorised density (see :class:DynamicCopula).

Examples:

>>> import rcopula as rc
>>> from rcopula.dynamic import DynamicCopula, fit_dynamic
>>> truth = DynamicCopula(rc.GaussianCopula(0.0), coefficients=(0.1, 0.4, 0.8))
>>> u = truth.simulate(800, random_state=0)
>>> res = fit_dynamic(u, rc.GaussianCopula(0.0))
>>> bool(res.loglik > res.constant_loglik)
True

fit_dcc ¶

fit_dcc(u: ArrayLike, *, df: float | None = None, start: tuple[float, float] = (0.03, 0.95)) -> DccResult

Estimate a correlation matrix that changes over time, for two or more series.

This is Engle's dynamic conditional correlation (DCC) model, applied to the copula of u. It is the multivariate answer to :class:DynamicCopula: instead of one parameter moving, the whole correlation matrix does, driven by the outer product of the last observation. Only two numbers, a and b, are estimated, however many series there are.

Parameters:

Name	Type	Description	Default
`u`	`array_like of float, shape (n, d)`	`d >= 2` series in time order (oldest first), ideally pseudo-observations strictly between 0 and 1. If any value is 0, 1 or outside that range, the data are converted with :func:`~rcopula.pseudo_obs` first.	required
`df`	`float or None`	Fit a t copula with this many degrees of freedom (a positive number) instead of a Gaussian one. Not estimated -- profile it by calling this over a grid, which is what practitioners do anyway because the df likelihood is flat.	`None`
`start`	`tuple of (float, float)`	Starting `(a, b)` for the optimiser; both non-negative with `a + b < 1`.	`(0.03, 0.95)`

Returns:

Type	Description
`DccResult`	`correlations` (shape `(n, d, d)`), coefficients `a` and `b`, `loglik`, the `unconditional` target matrix, `df`, `converged` and `n_obs`.

Raises:

Type	Description
`ValueError`	If `u` is not two-dimensional with at least two columns.

Notes

a and b are constrained to be non-negative with a + b < 1, which is sufficient (not necessary) for :math:Q_t to stay positive definite.

Examples:

>>> import numpy as np, rcopula as rc
>>> from rcopula.dynamic import fit_dcc
>>> rng = np.random.default_rng(0)
>>> u = rc.GaussianCopula(0.5, dim=3, dispstr="ex").rvs(400, random_state=0)
>>> res = fit_dcc(u)
>>> res.correlations.shape
(400, 3, 3)
>>> bool(0.0 <= res.persistence < 1.0)
True

rcopula.discrete ¶

Copulas for discrete and mixed margins.

Everything else in this package assumes continuous margins, where Sklar's theorem gives a unique copula and rank-based inference is exact. Neither holds for counts, ordinal scales, or a continuous variable paired with a binary one -- and the literature that quietly applies continuous machinery to them is large.

What actually breaks. Sklar's theorem still says a copula exists, but it is unique only on :math:\mathrm{Ran}\,F_1 \times \cdots \times \mathrm{Ran}\,F_d. For a Bernoulli margin that range is three points, so all a copula can be identified from is what happens at three points: infinitely many copulas give exactly the same joint distribution. This is not a small-sample problem and more data does not fix it. It means a fitted parameter is interpretable within a chosen family and comparisons across families are on much weaker ground than they look (Genest and Nešlehová 2007, which is worth reading before using any of this).

What still works. The joint distribution is perfectly well defined, and so is its likelihood -- as a finite difference of the copula rather than a derivative of it:

.. math::

P(X = x) = \sum_{j \in \{0,1\}^d} (-1)^{|j|}
           C\bigl(u_1^{(j_1)}, \dots, u_d^{(j_d)}\bigr),
\qquad u_k^{(0)} = F_k(x_k), \quad u_k^{(1)} = F_k(x_k^-).

That is exact, it is what :func:discrete_pmf computes, and maximising it is what :func:fit_discrete does. For mixed margins the two operations combine: differentiate along the continuous coordinates, difference along the discrete ones -- :func:mixed_pdf.

Ranks, and why they mislead. Ties break the correspondence between a sample rank correlation and the copula's own. Kendall's tau-b divides out the ties within each margin, so it still reaches 1 for comonotone identical margins -- but when the margins differ, no coupling can align their atoms and the ceiling drops: two Bernoullis at 0.1 and 0.9 cannot exceed 0.111 however strongly they are coupled. :func:tau_upper_bound computes it. The practical consequence is that inverting a sample tau to get a copula parameter, which is exact for continuous margins, is simply wrong here -- which is why :func:fit_discrete offers likelihood only.

The distributional transform (Ferguson 1967; Rüschendorf 2009) turns a discrete variable into an exactly uniform one by randomising within each atom. It is the honest version of jittering: :func:distributional_transform gives pseudo-observations that any continuous-margin method can consume, at the cost of the randomisation being part of the answer. Average over several draws.

============================ ================================================ :func:discrete_pmf Exact probability mass, by inclusion-exclusion. :func:mixed_pdf Density for any mix of discrete and continuous. :func:fit_discrete Maximum likelihood on the exact mass function. :func:distributional_transform Randomised pseudo-observations. :func:tau_upper_bound The largest Kendall tau these margins allow. :func:checkerboard The canonical member of the identified class. ============================ ================================================

Examples:

>>> import numpy as np, rcopula as rc
>>> from scipy import stats
>>> from rcopula.discrete import discrete_pmf
>>> margins = [stats.poisson(3.0), stats.poisson(2.0)]
>>> x = np.array([[3, 2], [0, 0], [5, 4]])
>>> mass = discrete_pmf(rc.GaussianCopula(0.6), x, margins)
>>> bool(np.all(mass > 0))
True

References

Genest, C. and Nešlehová, J. (2007). A primer on copulas for count data. ASTIN Bulletin 37(2), 475-515. The paper to read first; the source of the identifiability caveat. Rüschendorf, L. (2009). On the distributional transform, Sklar's theorem, and the empirical copula process. J. Statistical Planning and Inference 139(11), 3921-3927. Ferguson, T. S. (1967). Mathematical Statistics: A Decision Theoretic Approach. Academic Press. The distributional transform. Nikoloulopoulos, A. K. (2013). Copula-based models for multivariate discrete response data. In Copulae in Mathematical and Quantitative Finance, 231-249. The inclusion-exclusion likelihood. Song, P. X.-K., Li, M. and Yuan, Y. (2009). Joint regression analysis of correlated data using Gaussian copulas. Biometrics 65(1), 60-68. Denuit, M. and Lambert, P. (2005). Constraints on concordance measures in bivariate discrete data. J. Multivariate Analysis 93(1), 40-57. Where the attainable range of Kendall's tau comes from. Sun, T., Song, X. and Zhang, X. (2021). scDesign2: a transparent simulator for single-cell RNA sequencing data. Genome Biology 22, 163. A Gaussian copula with negative-binomial margins, at scale.

DiscreteMargin ¶

Bases: Protocol

The methods a discrete marginal distribution needs for the functions in this module.

This is a typing protocol, not a class to instantiate: any object with cdf, pmf and ppf methods qualifies. It is satisfied by every discrete scipy.stats frozen distribution (poisson, nbinom, binom, geom, randint, ...), e.g. scipy.stats.poisson(3.0).

The methods are, for an array of values x or probabilities q:

cdf(x) -- P(X <= x), array of float in [0, 1];
pmf(x) -- P(X = x), array of float in [0, 1];
ppf(q) -- the quantile function, the smallest x with cdf(x) >= q.

Because this is a runtime_checkable protocol, isinstance(margin, DiscreteMargin) checks only that the three methods exist.

DiscreteFitResult `dataclass` ¶

DiscreteFitResult(copula: Copula, params: NDArray[float64], loglik: float, n_obs: int, converged: bool, independent_loglik: float, message: str = '')

The result of fitting a copula to discrete data: the fitted copula plus fit statistics.

Returned by :func:fit_discrete; you normally do not construct it yourself. Call :meth:summary for a readable report.

Attributes:

Name	Type	Description
`copula`	`Copula`	The copula at the estimated parameters.
`params`	`numpy.ndarray of float, shape (p,)`	All of the copula's parameters (free and fixed), in the order of `copula.param_names`.
`loglik`	`float`	Log-likelihood at the estimate.
`n_obs`	`int`	Number of observations (rows) used.
`converged`	`bool`	Whether the optimiser reported success.
`independent_loglik`	`float`	The same margins under independence, for a likelihood ratio.
`message`	`str, default ""`	The optimiser's status message; useful when `converged` is False.

n_params `property` ¶

n_params: int

Number of copula parameters that were estimated (not held fixed).

Returns:

Type	Description
`int`

aic `property` ¶

aic: float

Akaike information criterion, 2 * n_params - 2 * loglik; lower is better.

Returns:

Type	Description
`float`

bic `property` ¶

bic: float

Bayesian information criterion, n_params * log(n_obs) - 2 * loglik; lower is better.

Returns:

Type	Description
`float`

independence_test ¶

independence_test() -> tuple[float, float]

Test whether the variables are dependent at all, against the independence model.

This is a likelihood-ratio test. A small p-value is evidence of dependence.

Unlike the constancy test in :mod:rcopula.dynamic, this null is interior for every family here, so the chi-squared reference is the usual asymptotic one.

Returns:

Name	Type	Description
`statistic`	`float`	`2 * (loglik - independent_loglik)`, floored at zero.
`pvalue`	`float`	Upper-tail probability of a chi-squared distribution with `n_params` degrees of freedom.

summary ¶

summary() -> str

Return a printable text report of the fit.

Includes the parameter estimates, log-likelihood, AIC/BIC, the likelihood-ratio test against independence, and a warning if the optimiser did not converge.

Returns:

Type	Description
`str`

discrete_pmf ¶

discrete_pmf(copula: Copula, x: ArrayLike, margins: list[Any]) -> NDArray[np.float64]

Compute the exact probability of each observed row of counts under a copula model.

Use this when every variable is discrete (counts, categories coded as integers, ...). It returns :math:P(X_1 = x_1, \dots, X_d = x_d) for each row, combining the copula (the dependence) with the given marginal distributions.

Computes the :math:2^d-term inclusion-exclusion sum in the module docstring. This is the C-volume of the rectangle :math:\prod_k (F_k(x_k^-), F_k(x_k)], which is what a copula's :math:d-increasing property guarantees is non-negative -- so a negative result here means the copula is not one.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Any `d`-dimensional copula from this package, e.g. `GaussianCopula(0.6)`.	required
`x`	`array_like of float or int, shape (n, d) or (d,)`	Observed values, on the margins' own scale (e.g. the counts themselves, not ranks). A single row may be given as shape `(d,)`.	required
`margins`	`list of DiscreteMargin, length d`	One frozen discrete distribution per column, e.g. `[scipy.stats.poisson(3.0), scipy.stats.poisson(2.0)]`.	required

Returns:

Type	Description
`numpy.ndarray of float, shape (n,)`	Probability of each row, in `[0, 1]`. Tiny negative values from floating-point cancellation are set to zero.

Raises:

Type	Description
`ValueError`	If `x` does not have `copula.dim` columns, or `margins` does not have `copula.dim` entries.

Notes

Cost is :math:2^d copula CDF evaluations, each vectorised over observations. Past about :math:d = 15 that is the binding constraint and a composite-likelihood approach (pairs, or a vine) is the usual answer.

Examples:

The mass function must sum to one over the whole lattice:

>>> import numpy as np, rcopula as rc
>>> from scipy import stats
>>> from rcopula.discrete import discrete_pmf
>>> margins = [stats.poisson(2.0), stats.poisson(3.0)]
>>> grid = np.array([[i, j] for i in range(40) for j in range(45)])
>>> total = discrete_pmf(rc.ClaytonCopula(2.0), grid, margins).sum()
>>> bool(abs(total - 1.0) < 1e-9)
True

And it must reproduce the margins when summed over the other coordinate:

>>> marginal = discrete_pmf(rc.ClaytonCopula(2.0), grid, margins).reshape(40, 45).sum(axis=1)
>>> bool(np.allclose(marginal, margins[0].pmf(np.arange(40)), atol=1e-9))
True

mixed_pdf ¶

mixed_pdf(copula: Copula, x: ArrayLike, margins: list[Any], discrete: ArrayLike, *, step: float = 1e-05) -> NDArray[np.float64]

Compute the likelihood of each row when some variables are discrete and one is continuous.

For example, a continuous measurement paired with a yes/no outcome. The result is a probability for the discrete coordinates and a density for the continuous one, which is what a likelihood for such data needs.

The joint density of a mixed vector is a partial derivative: differentiate the copula along the continuous coordinates and difference it along the discrete ones, then multiply by the continuous marginal densities.

With one discrete coordinate this reduces to a difference of h-functions,

.. math::

f(x_1, x_2) = f_1(x_1)\left[
    \frac{\partial C}{\partial u_1}(u_1, F_2(x_2))
    - \frac{\partial C}{\partial u_1}(u_1, F_2(x_2^-))
\right],

which is the form used in the transportation and biostatistics literature for joining a discrete choice to a continuous response.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Any `d`-dimensional copula from this package.	required
`x`	`array_like of float, shape (n, d) or (d,)`	Observed values on the margins' own scale.	required
`margins`	`list of frozen scipy.stats distributions, length d`	Discrete ones must provide `cdf` and `pmf`; continuous ones `cdf` and `pdf`.	required
`discrete`	`array_like of bool, shape (d,)`	Which coordinates are discrete (`True`) and which continuous (`False`). At most one may be `False` unless all are.	required
`step`	`float`	Finite-difference step, on the copula scale, for the derivative along the continuous coordinate. Used whenever that derivative has no closed form: for every copula with `dim > 2`, and for bivariate families other than the Archimedean, Gaussian and Student t ones (whose conditional CDFs are analytic, so `step` does not affect them). The step is shrunk near the edges of the unit interval so the difference stays inside it. Must be in `(0, 0.5)`.	`1e-5`

Returns:

Type	Description
`numpy.ndarray of float, shape (n,)`	The joint mass-density of each row, non-negative.

Raises:

Type	Description
`ValueError`	If `discrete`, `x` or `margins` does not match `copula.dim`, or `step` is not in `(0, 0.5)`.
`NotImplementedError`	If more than one coordinate is continuous while at least one is discrete.

Notes

With no discrete coordinates this is the ordinary copula density times the marginal densities, and with all of them it is :func:discrete_pmf; both limits are checked in the test suite. For a bivariate Archimedean, Gaussian or t copula the continuous derivative goes through the analytic :func:~rcopula.conditional_cdf; otherwise it is a central difference of the copula CDF with step step, which is only as accurate as that CDF (in d > 2 the elliptical CDFs are themselves numerical integrals, so a larger step than the default can be the more accurate choice there).

Examples:

A continuous margin paired with a Bernoulli one. The density integrates and sums to one:

>>> import numpy as np, rcopula as rc
>>> from scipy import stats
>>> from rcopula.discrete import mixed_pdf
>>> margins = [stats.norm(), stats.bernoulli(0.3)]
>>> grid = np.linspace(-8, 8, 2001)
>>> rows = np.concatenate([
...     np.column_stack([grid, np.zeros_like(grid)]),
...     np.column_stack([grid, np.ones_like(grid)]),
... ])
>>> values = mixed_pdf(rc.GaussianCopula(0.5), rows, margins, [False, True])
>>> total = np.trapezoid(values[: grid.size], grid) + np.trapezoid(values[grid.size :], grid)
>>> bool(abs(total - 1.0) < 1e-6)
True

discrete_loglik ¶

discrete_loglik(copula: Copula, x: ArrayLike, margins: list[Any]) -> float

Score how well a copula with discrete margins explains the data (higher is better).

Returns the log-likelihood: the sum over rows of the log of :func:discrete_pmf. Comparing it across parameter values (or, with care, families) shows which explains the data best; :func:fit_discrete maximises it.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	Any `d`-dimensional copula from this package.	required
`x`	`array_like of float or int, shape (n, d)`	Observed values on the margins' own scale.	required
`margins`	`list of DiscreteMargin, length d`	One frozen discrete distribution per column.	required

Returns:

Type	Description
`float`	The log-likelihood. Rows with probability below `1e-300` contribute `log(1e-300)` (about -691) rather than minus infinity.

Raises:

Type	Description
`ValueError`	As for :func:`discrete_pmf`.

Examples:

>>> import numpy as np, rcopula as rc
>>> from scipy import stats
>>> from rcopula.discrete import discrete_loglik
>>> margins = [stats.poisson(2.0), stats.poisson(2.0)]
>>> x = rc.CopulaDistribution(rc.GaussianCopula(0.7), margins).rvs(300, random_state=0)
>>> strong = discrete_loglik(rc.GaussianCopula(0.7), x, margins)
>>> weak = discrete_loglik(rc.GaussianCopula(0.0), x, margins)
>>> bool(strong > weak)
True

fit_discrete ¶

fit_discrete(x: ArrayLike, copula: Copula, margins: list[Any], *, start: ArrayLike | None = None) -> DiscreteFitResult

Estimate a copula's parameters from discrete data (counts, codes), by maximum likelihood.

The margins are taken as given -- fit them separately, which is the inference-functions-for-margins two-step and is what everyone does. The copula parameter is then the only unknown.

Parameters:

Name	Type	Description	Default
`x`	`array_like of float or int, shape (n, d)`	Observed counts or codes, on the margins' own scale.	required
`copula`	`Copula`	The family to fit, with `dim == d`. Its current parameters are the starting point, and any parameters it marks as fixed stay fixed.	required
`margins`	`list of DiscreteMargin, length d`	One already-fitted frozen discrete distribution per column.	required
`start`	`array_like of float or None`	Starting values for the free parameters only, overriding the copula's current ones: exactly one finite value per free parameter. `None` uses the copula's current values.	`None`

Returns:

Type	Description
`DiscreteFitResult`	The fitted copula, its parameters, log-likelihood, convergence flag and the log-likelihood of the same margins under independence. If the copula has no free parameters nothing is optimised: the result holds the copula as given and its log-likelihood, and `independent_loglik` is still the independence copula's, so the likelihood ratio compares the given copula with independence.

Raises:

Type	Description
`ValueError`	If `x` or `margins` does not match `copula.dim`, or `start` does not have one finite value per free parameter.

Notes

Rank-based estimation (method="itau") is not offered here on purpose: with ties, the sample Kendall tau does not estimate the copula's tau, and inverting it produces a parameter with no defensible interpretation. See :func:tau_upper_bound for the size of the distortion.

Examples:

>>> import numpy as np, rcopula as rc
>>> from scipy import stats
>>> from rcopula.discrete import fit_discrete
>>> margins = [stats.poisson(4.0), stats.poisson(4.0)]
>>> x = rc.CopulaDistribution(rc.GaussianCopula(0.6), margins).rvs(2000, random_state=0)
>>> result = fit_discrete(x, rc.GaussianCopula(0.0), margins)
>>> bool(abs(result.params[0] - 0.6) < 0.06)
True

distributional_transform ¶

distributional_transform(x: ArrayLike, margins: list[Any], *, random_state: Any = None, replicates: int = 1) -> NDArray[np.float64]

Turn discrete observations into values between 0 and 1 by adding controlled randomness.

Each count is replaced by a random point inside the probability interval it occupies, so the result is exactly uniform and can be fed to any method that assumes continuous data (fitting, plots, tests). The randomness is part of the answer, so repeat with different seeds and compare.

The distributional transform randomises within each atom,

.. math:: U = F(X^-) + V\,\bigl(F(X) - F(X^-)\bigr), \qquad V \sim U(0,1),

with :math:V independent of :math:X. The result is exactly uniform, not approximately -- which is what separates this from ad-hoc jittering -- and its copula is a copula of :math:X, chosen at random from the identified class.

Parameters:

Name	Type	Description	Default
`x`	`array_like of float or int, shape (n, d) or (d,)`	Observed values on the margins' own scale.	required
`margins`	`list of frozen scipy.stats distributions, length d`	One per column. A margin with a `pmf` method is treated as discrete; any other (continuous) margin is simply passed through its `cdf`.	required
`random_state`	`(int, Generator or None)`	Seed or generator for the randomisation.	`None`
`replicates`	`int`	Draw this many independent transforms (at least 1) and return their average on the copula scale. Averaging reduces the randomisation noise but biases the result towards the middle of each atom, so the default is 1 and anything else is a deliberate trade.	`1`

Returns:

Type	Description
`numpy.ndarray of float, shape (n, d)`	Pseudo-observations, clipped to `[1e-12, 1 - 1e-12]`.

Raises:

Type	Description
`ValueError`	If `margins` does not have one entry per column of `x`, or `replicates < 1`.

Examples:

>>> import numpy as np
>>> from scipy import stats
>>> from rcopula.discrete import distributional_transform
>>> x = stats.poisson(3.0).rvs(20000, random_state=0)[:, None]
>>> u = distributional_transform(x, [stats.poisson(3.0)], random_state=0)
>>> bool(abs(u.mean() - 0.5) < 0.01)          # exactly uniform, so mean 1/2
True
>>> bool(abs(u.var() - 1 / 12) < 0.005)
True

tau_upper_bound ¶

tau_upper_bound(margins: list[Any], *, support: int = 200) -> float

Compute the largest rank correlation (Kendall's tau-b) two discrete margins allow.

Use it to judge a measured correlation between discrete variables: a tau of 0.1 may already be the strongest dependence the margins permit.

Evaluates :math:\tau_b at the comonotone coupling, which is where it is maximised. Since :math:\tau_b divides by :math:\sqrt{(1-\sum_i p_i^2)(1-\sum_j q_j^2)} it already corrects for the ties within each margin, and so reaches 1 when the two margins are identical -- comonotone Poisson(1) pairs really do give :math:\tau_b = 1.

The ceiling bites when the margins differ, because then no coupling can align their atoms. Two Bernoullis with success probabilities 0.1 and 0.9 cannot exceed 0.111 however they are coupled, so reading a fitted 0.1 there as "weak dependence" inverts the truth: it is as strong as the margins permit.

Parameters:

Name	Type	Description	Default
`margins`	`list of DiscreteMargin, length 2`	Two frozen discrete distributions supported on the non-negative integers `0, 1, 2, ...`.	required
`support`	`int`	How far up the lattice to evaluate: values `0, ..., support` are used. The tail beyond this contributes less than its mass, so the default is generous for anything with a finite mean.	`200`

Returns:

Type	Description
`float`	The attainable maximum of Kendall's tau-b, in `[0, 1]`. Returns 0.0 if either margin puts all its mass on one value.

Raises:

Type	Description
`ValueError`	If `margins` does not have exactly two entries, `support` is not a non-negative integer, or a margin puts mass on negative or non-integer values (which the lattice `0, ..., support` would silently miss).

Examples:

Identical margins reach 1; mismatched ones do not come close:

>>> from scipy import stats
>>> from rcopula.discrete import tau_upper_bound
>>> round(tau_upper_bound([stats.poisson(1.0), stats.poisson(1.0)]), 4)
1.0
>>> round(tau_upper_bound([stats.bernoulli(0.1), stats.bernoulli(0.9)]), 4)
0.1111
>>> round(tau_upper_bound([stats.poisson(3.0), stats.nbinom(4, 0.5)]), 4)
0.9438

checkerboard ¶

checkerboard(copula: Copula, margins: list[Any], *, support: int = 60) -> NDArray[np.float64]

Tabulate the probability of every pair of count values under a bivariate copula model.

The table is the checkerboard copula's mass on the lattice induced by the margins, and is the natural thing to plot as a picture of a copula fitted to discrete data.

Since the copula is identified only on :math:\mathrm{Ran}\,F_1 \times \mathrm{Ran}\,F_2, an infinite family of copulas fits any discrete data equally well. The checkerboard member spreads each cell's mass uniformly over its rectangle: it is the canonical representative, it is the one the distributional transform targets on average, and it is the one to plot when a picture of "the" fitted copula is wanted.

Parameters:

Name	Type	Description	Default
`copula`	`Copula`	A bivariate copula (`dim == 2`).	required
`margins`	`list of DiscreteMargin, length 2`	Two frozen discrete distributions supported on the non-negative integers.	required
`support`	`int`	Lattice extent: values `0, ..., support` are tabulated in each coordinate.	`60`

Returns:

Type	Description
`numpy.ndarray of float, shape (support + 1, support + 1)`	Cell probabilities; entry `[i, j]` is :math:`P(X_1 = i, X_2 = j)`. They sum to one up to the tail truncation.

Raises:

Type	Description
`ValueError`	If the copula is not bivariate, `margins` does not have two entries, `support` is not a non-negative integer, or a margin puts mass on negative or non-integer values.

Examples:

>>> import numpy as np, rcopula as rc
>>> from scipy import stats
>>> from rcopula.discrete import checkerboard
>>> mass = checkerboard(rc.GaussianCopula(0.5), [stats.poisson(3.0)] * 2)
>>> bool(abs(mass.sum() - 1.0) < 1e-8)
True
>>> bool(np.all(mass >= 0))
True

rcopula.insurance ¶

Insurance, actuarial and operational-risk utilities.

Tools to simulate total claims, size reinsurance layers and catastrophe bonds, and compute operational-risk capital, with dependence between lines or cells described by a copula.

Insurance is where copulas are least optional. An insurer's aggregate loss is a compound distribution -- a random number of claims, each of random size -- and compound distributions have no closed form worth using, no elliptical structure, and margins that are wildly non-normal. There is no covariance matrix that summarises how two lines of business co-move; there is only the joint distribution, and a copula is how you write one down.

The regulatory question is also the tail, not the middle. Solvency II asks for the 99.5% one-year VaR and Basel's operational-risk framework asks for 99.9%. At those quantiles the difference between "correlated" and "tail-dependent" is the whole answer.

================================= =========================================== :func:aggregate_loss Compound frequency-severity distribution. :func:operational_risk_capital Loss-distribution approach across cells. :func:excess_of_loss Reinsurance layer recovery. :func:layer_statistics Expected loss and attachment probabilities. :func:reinsurance_premium Premium for a layer, with loading. :func:catastrophe_bond Cat bond / ILS payoff analysis. ================================= ===========================================

References

Klugman, S. A., Panjer, H. H. and Willmot, G. E. (2019). Loss Models: From Data to Decisions, 5th ed. Wiley. Frees, E. W. and Valdez, E. A. (1998). Understanding relationships using copulas. North American Actuarial Journal 2(1), 1-25. The loss-ALAE study that introduced copulas to actuarial work. Frachot, A., Roncalli, T. and Salomon, E. (2004). The correlation problem in operational risk. Groupe de Recherche Operationnelle, Credit Lyonnais. Embrechts, P., Kluppelberg, C. and Mikosch, T. (1997). Modelling Extremal Events for Insurance and Finance. Springer. Cummins, J. D. (2008). CAT bonds and other risk-linked securities. Risk Management and Insurance Review 11(1), 23-47.

Variate ¶

Bases: Protocol

Anything you can draw random numbers from: the type of a frequency or severity input.

You rarely implement this yourself. Any scipy.stats frozen distribution -- stats.poisson(50) for a claim count, stats.lognorm(1.5, scale=1000) for a claim size -- already qualifies. It is a structural type (a :class:typing.Protocol): any object with a matching rvs method is accepted, with no need to subclass.

Methods:

Name	Description
`rvs`	Draw `size` random values. `size` is an int or tuple of ints; `random_state` is a :class:`numpy.random.Generator`. Must return an array_like of numbers with shape `size`.

Notes

Deliberately weaker than :class:~rcopula.distribution.Margin, which needs cdf/pdf/ppf. A compound distribution is built by simulation, so all that is required here is rvs -- satisfied by every scipy.stats frozen distribution, discrete or continuous, and by anything else that draws.

LayerStatistics ¶

Bases: NamedTuple

Headline numbers for one reinsurance layer, as returned by :func:layer_statistics.

A named tuple, so fields can be read by name (stats.expected_loss) or unpacked like an ordinary tuple. All monetary fields are in the same currency units as the simulated losses they were computed from.

Attributes:

Name	Type	Description
`attachment`	`float`	Loss level at which the layer starts paying (the cedant's retention).
`limit`	`float`	Most the layer can pay; the layer covers losses from `attachment` to `attachment + limit`.
`expected_loss`	`float`	Average amount the layer pays per simulated period, in currency units.
`attachment_probability`	`float`	Probability (0 to 1) that the loss exceeds `attachment`, i.e. the layer pays anything at all.
`exhaustion_probability`	`float`	Probability (0 to 1) that the loss reaches `attachment + limit`, i.e. the layer pays out in full.
`expected_loss_ratio`	`float`	`expected_loss / limit`: expected payout as a fraction of the limit, the usual way layers of different size are compared.

Notes

The repr follows the market convention "<limit> xs <attachment>", so a layer paying 5m above a 10m retention prints as 5e+06 xs 1e+07.

aggregate_loss ¶

aggregate_loss(frequency: Variate, severity: Variate, n: int = 100000, random_state: Generator | int | None = None) -> NDArray[np.float64]

Simulate total yearly losses as a random number of claims of random size.

This is the compound (frequency-severity) aggregate loss: draw how many claims occur, draw how large each one is, add them up. Use it to turn a claim-count model and a claim-size model into a distribution of total loss for a policy, a line of business, or an operational-risk cell, which you can then feed to :func:~rcopula.risk.value_at_risk, :func:layer_statistics and friends.

Parameters:

Name	Type	Description	Default
`frequency`	`Variate (e.g. scipy.stats frozen discrete distribution)`	Number of claims per period, e.g. `scipy.stats.poisson(50)`. Must produce non-negative integers.	required
`severity`	`Variate (e.g. scipy.stats frozen continuous distribution)`	Size of each individual claim, in currency units, e.g. `scipy.stats.lognorm(1.5, scale=1000)`.	required
`n`	`int`	Number of simulated periods (scenarios).	`100_000`
`random_state`	`(int, Generator or None)`	Seed or generator for reproducibility. A Generator is used (and advanced) as is; an int or None seeds a fresh one.	`None`

Returns:

Type	Description
`numpy.ndarray of float, shape (n,)`	Total loss in each simulated period, in the units of `severity`. Periods with zero claims have a total of exactly 0.

Raises:

Type	Description
`ValueError`	If `frequency` produces a negative claim count.

Notes

:math:S = \sum_{i=1}^{N} X_i with :math:N from frequency and each :math:X_i from severity. This is the actuarial workhorse: no closed form in general, but its first two moments are exact and worth checking against,

.. math:: \mathbb{E}[S] = \mathbb{E}[N]\,\mathbb{E}[X], \qquad \mathrm{Var}[S] = \mathbb{E}[N]\,\mathrm{Var}[X] + \mathrm{Var}[N]\,\mathbb{E}[X]^2.

Examples:

The compound moments come out right:

>>> import numpy as np
>>> from scipy import stats
>>> from rcopula.insurance import aggregate_loss
>>> freq, sev = stats.poisson(40), stats.lognorm(0.8, scale=1000)
>>> s = aggregate_loss(freq, sev, n=60_000, random_state=0)
>>> expected = freq.mean() * sev.mean()
>>> bool(abs(s.mean() / expected - 1) < 0.02)
True

The aggregate is right-skewed even when the severity is only mildly so:

>>> bool(stats.skew(s) > 0)
True

operational_risk_capital ¶

operational_risk_capital(cells: list[tuple[Variate, Variate]], copula: Copula | None = None, alpha: float = 0.999, n: int = 100000, random_state: Generator | int | None = None) -> dict[str, float]

How much capital to hold for operational losses across several risk cells.

This is the loss-distribution approach (LDA). You describe each cell -- for example "retail banking x external fraud" -- by a claim-count and a loss-size distribution, choose a copula for how the cells' bad years coincide, and get back the 99.9% aggregate loss, the capital it implies, and how much diversification credit the dependence assumption allows.

Parameters:

Name	Type	Description	Default
`cells`	`list of (Variate, Variate) tuples, length d`	One `(frequency, severity)` pair of frozen distributions per cell, as for :func:`aggregate_loss`.	required
`copula`	`Copula or None`	Dependence across cells; its `dim` must equal `len(cells)`. `None` means independence, which is the assumption that produces the largest diversification credit.	`None`
`alpha`	`float`	Confidence level, between 0 and 1; 0.999 is the Basel operational-risk standard.	`0.999`
`n`	`int`	Number of simulated years.	`100_000`
`random_state`	`(int, Generator or None)`	Seed or generator for reproducibility.	`None`

Returns:

Type Description

dict of str to float

All amounts in the units of the severity distributions.

"var" Value-at-risk of the total loss at level alpha. "expected_shortfall" Average total loss in the worst 1 - alpha of years. "expected_loss" Mean total loss (provisioned, not capitalised). "capital" var - expected_loss: the capital requirement. "standalone_capital" Sum over cells of each cell's own VaR minus its own mean -- the capital if no diversification were recognised. "diversification_benefit" standalone_capital - capital.

Raises:

Type	Description
`ValueError`	If `cells` is empty, or `copula.dim` does not equal `len(cells)`.

Notes

Each cell -- a business line crossed with an event type -- gets its own compound distribution. The cells are then combined under copula, by rank reordering, so each cell's approved loss distribution is preserved exactly while the dependence becomes the copula's. That matters in practice: cell distributions are validated and signed off individually, and a joint model that quietly changes them is not usable.

Capital is the difference between the aggregate VaR and the expected loss, which is the Basel definition -- expected loss is provisioned, not capitalised.

Examples:

Dependence between cells consumes the diversification credit:

>>> from scipy import stats
>>> from rcopula import GaussianCopula
>>> from rcopula.insurance import operational_risk_capital
>>> cells = [(stats.poisson(30), stats.lognorm(1.2, scale=500))] * 3
>>> free = operational_risk_capital(cells, None, n=40_000, random_state=0)
>>> tied = operational_risk_capital(
...     cells, GaussianCopula(0.8, dim=3), n=40_000, random_state=0)
>>> bool(tied["capital"] > free["capital"])
True
>>> bool(free["diversification_benefit"] > tied["diversification_benefit"])
True

excess_of_loss ¶

excess_of_loss(losses: ArrayLike, attachment: float, limit: float) -> NDArray[np.float64]

How much a reinsurance layer pays back for each loss.

An excess-of-loss layer pays nothing until the loss passes the attachment point, then pays the excess, up to at most limit. Use it to turn simulated gross losses into the reinsurer's recoveries (and, by subtraction, the insurer's net losses).

Parameters:

Name	Type	Description	Default
`losses`	`array_like of float, any shape`	Gross losses, in currency units (positive numbers are losses).	required
`attachment`	`float`	Retention: the loss level at which the layer starts paying. Must be `>= 0`.	required
`limit`	`float`	Maximum amount the layer pays. Must be `> 0`.	required

Returns:

Type	Description
numpy.ndarray of float, same shape as ``losses``	Recovery for each loss, between 0 and `limit`.

Raises:

Type	Description
`ValueError`	If `attachment < 0` or `limit <= 0`.

Notes

A "limit excess of attachment" layer pays :math:\min(\max(L - a, 0), \ell): nothing until the cedant's retention is exhausted, then pound for pound, then nothing once the limit is used up. Structurally identical to a CDO tranche, which is not a coincidence -- both are call spreads on an aggregate loss.

Examples:

>>> import numpy as np
>>> from rcopula.insurance import excess_of_loss
>>> losses = np.array([0.0, 50.0, 150.0, 400.0, 1000.0])
>>> excess_of_loss(losses, attachment=100.0, limit=200.0)
array([  0.,   0.,  50., 200., 200.])

layer_statistics ¶

layer_statistics(losses: ArrayLike, attachment: float, limit: float) -> LayerStatistics

Expected payout and the chance of being hit, for one reinsurance layer.

Summarises a limit excess of attachment layer against a set of simulated losses: how much it pays on average, how often it pays anything, and how often it pays out in full.

Parameters:

Name	Type	Description	Default
`losses`	`array_like of float, shape (n,)`	Simulated gross losses, in currency units (positive numbers are losses), e.g. from :func:`aggregate_loss`.	required
`attachment`	`float`	Loss level at which the layer starts paying. Must be `>= 0`.	required
`limit`	`float`	Maximum amount the layer pays. Must be `> 0`.	required

Returns:

Type	Description
`LayerStatistics`	Named tuple with fields `attachment`, `limit`, `expected_loss`, `attachment_probability`, `exhaustion_probability` and `expected_loss_ratio`; see :class:`LayerStatistics`.

Raises:

Type	Description
`ValueError`	If `attachment < 0` or `limit <= 0`.

Examples:

>>> from scipy import stats
>>> from rcopula.insurance import aggregate_loss, layer_statistics
>>> s = aggregate_loss(stats.poisson(40), stats.lognorm(1.0, scale=1000),
...                    n=40_000, random_state=0)
>>> low = layer_statistics(s, 60_000, 20_000)
>>> high = layer_statistics(s, 150_000, 20_000)
>>> bool(low.attachment_probability > high.attachment_probability)
True
>>> bool(low.expected_loss > high.expected_loss)
True

reinsurance_premium ¶

reinsurance_premium(losses: ArrayLike, attachment: float, limit: float, loading: float = 0.25, method: str = 'expected_value', alpha: float = 0.99) -> float

What to charge for a reinsurance layer: expected payout plus a risk loading.

Computes the reinsurer's recoveries on the layer from simulated losses (see :func:excess_of_loss) and applies one of three classical premium principles. Use the tail-based "expected_shortfall" method for high, remote layers, where the mean alone understates the risk being sold.

Parameters:

Name	Type	Description	Default
`losses`	`array_like of float, shape (n,)`	Simulated gross losses, in currency units (positive numbers are losses).	required
`attachment`	`float`	Loss level at which the layer starts paying. Must be `>= 0`.	required
`limit`	`float`	Maximum amount the layer pays. Must be `> 0`.	required
`loading`	`float`	Risk loading :math:`\theta`. Its meaning depends on `method` (a markup on the mean, a multiple of the standard deviation, or a weight on the tail); see below.	`0.25`
`method`	`str`	Premium principle to apply: one of `"expected_value"`, `"standard_deviation"` or `"expected_shortfall"`.	`"expected_value"`
`alpha`	`float`	Confidence level for the expected shortfall, between 0 and 1. Only used when `method="expected_shortfall"`.	`0.99`

Returns:

Type	Description
`float`	Premium, in the same currency units as `losses`.

Raises:

Type	Description
`ValueError`	If `method` is not one of the three names above, or if `attachment < 0` or `limit <= 0`.

Notes

Three classical principles:

expected_value :math:(1 + \theta)\,\mathbb{E}[R]. Simple, and the usual starting point, but it charges the same for two layers with equal mean and wildly different volatility. standard_deviation :math:\mathbb{E}[R] + \theta\,\mathrm{sd}(R). Prices the uncertainty. expected_shortfall :math:(1-\theta)\mathbb{E}[R] + \theta\,\mathrm{ES}_\alpha(R). Prices the tail, which is what a high layer is actually selling.

Examples:

Tail-sensitive principles charge more for a volatile layer:

>>> from scipy import stats
>>> from rcopula.insurance import aggregate_loss, reinsurance_premium
>>> s = aggregate_loss(stats.poisson(40), stats.lognorm(1.2, scale=1000),
...                    n=40_000, random_state=0)
>>> ev = reinsurance_premium(s, 100_000, 50_000, method="expected_value")
>>> es = reinsurance_premium(s, 100_000, 50_000, method="expected_shortfall")
>>> bool(es > ev)
True

catastrophe_bond ¶

catastrophe_bond(losses: ArrayLike, attachment: float, exhaustion: float, coupon: float = 0.06, risk_free: float = 0.03, maturity: float = 1.0) -> dict[str, float]

Expected loss, hit probability and relative value of a catastrophe bond.

A cat bond pays a coupon but loses principal if the sponsor's catastrophe losses are large. This function takes simulated sponsor losses and an attachment/exhaustion range and returns the figures the insurance-linked securities (ILS) market quotes, so different deals can be compared on risk and price.

Parameters:

Name	Type	Description	Default
`losses`	`array_like of float, shape (n,)`	Simulated sponsor losses over the whole term of the bond (one number per scenario covering `maturity` years), in currency units (positive numbers are losses).	required
`attachment`	`float`	Loss level at which principal starts to be written down. Must be `>= 0`.	required
`exhaustion`	`float`	Loss level at which all principal is lost. Must be greater than `attachment`.	required
`coupon`	`float`	Annual coupon rate as a decimal (0.06 = 6%).	`0.06`
`risk_free`	`float`	Annual risk-free rate as a decimal, continuously compounded for discounting.	`0.03`
`maturity`	`float`	Term in years, `> 0`. The whole-term expected loss is divided by `maturity` to annualise it before it is compared with the annual spread.	`1.0`

Returns:

Type Description

dict of str to float

Rates and losses are fractions of principal (0.02 = 2%).

"expected_loss" Average fraction of principal lost over the whole term. "annual_expected_loss" expected_loss / maturity: the per-year expected loss the ILS market quotes. "attachment_probability" Probability that the loss exceeds attachment (any principal lost). "exhaustion_probability" Probability that the loss reaches exhaustion (all principal lost). "conditional_severity" Average fraction of principal lost, given that some is lost (0.0 if no scenario loses principal). "spread" coupon - risk_free. "multiple" spread / annual_expected_loss (inf if the expected loss is 0). Both numbers are per year, so the ratio does not depend on the term. "expected_return" spread - annual_expected_loss: annual excess return net of expected losses. "fair_price" exp(-risk_free * maturity) * (1 + coupon * maturity - expected_loss), per unit of principal.

Raises:

Type	Description
`ValueError`	Unless `0 <= attachment < exhaustion` and `maturity > 0`.

Notes

Principal is written down linearly between attachment and exhaustion, so the investor is short a call spread on the sponsor's loss -- the mirror image of :func:excess_of_loss.

Returns the three numbers the ILS market quotes: expected loss (the fraction of principal expected to be lost), attachment probability, and the multiple -- spread divided by expected loss -- which is how relative value is judged across deals. The spread is an annual rate while losses cover the whole term, so the expected loss is annualised (divided by maturity) before the multiple and expected return are formed; with the default maturity=1 the two coincide.

Examples:

A remote layer has a low expected loss and therefore a high multiple:

>>> from scipy import stats
>>> from rcopula.insurance import aggregate_loss, catastrophe_bond
>>> s = aggregate_loss(stats.poisson(20), stats.lognorm(1.5, scale=2000),
...                    n=60_000, random_state=0)
>>> near = catastrophe_bond(s, 60_000, 120_000)
>>> far = catastrophe_bond(s, 200_000, 300_000)
>>> bool(far["expected_loss"] < near["expected_loss"])
True
>>> bool(far["multiple"] > near["multiple"])
True

Applications¶

rcopula.risk ¶

RiskSummary ¶

value_at_risk ¶

expected_shortfall ¶

rank_reorder ¶

simulate_losses ¶

diversification_benefit ¶

risk_contributions ¶

covar ¶

delta_covar ¶

marginal_expected_shortfall ¶

stress_scenario ¶

rcopula.credit ¶

default_indicators ¶

default_times ¶

portfolio_loss ¶

tranche_loss ¶

tranche_expected_loss ¶

tranche_spread ¶

nth_to_default_probability ¶

vasicek_loss_cdf ¶

implied_correlation ¶

rcopula.derivatives ¶

MonteCarloPrice ¶

SmileMargin ¶

cdf ¶

ppf ¶

pdf ¶

CmsLeg ¶

lognormal_terminal ¶

black76 ¶

margrabe ¶

kirk_spread ¶

implied_volatility ¶

basket_option ¶

rainbow_option ¶

spread_option ¶

cms_convexity_adjustment ¶

cms_margin ¶

cms_spread_option ¶

basket_implied_vol ¶

rcopula.portfolio ¶

BacktestResult ¶

EfficientFrontier ¶

mispricing_index ¶

pairs_signal ¶

backtest_pairs ¶

simulate_returns ¶

mean_cvar_weights ¶

min_variance_weights ¶

efficient_frontier ¶

rcopula.statarb ¶

multivariate_spearman ¶

diagonal_distance ¶

tail_concentration ¶

select_pairs ¶

select_partners ¶

Copula-GARCH margins¶

rcopula.garch ¶

GarchResult dataclass ¶

persistence property ¶

unconditional_mean property ¶

unconditional_vol property ¶

half_life property ¶

n_params property ¶

aic property ¶

bic property ¶

innovation ¶

forecast_variance ¶

forecast_mean ¶

forecast_vol ¶

CopulaGarch ¶

dim property ¶

names property ¶

fit classmethod ¶

simulate ¶

forecast ¶

forecast_risk ¶

summary ¶

GarchResult `dataclass` ¶

persistence `property` ¶

unconditional_mean `property` ¶

unconditional_vol `property` ¶

half_life `property` ¶

n_params `property` ¶

aic `property` ¶

bic `property` ¶

dim `property` ¶

names `property` ¶

fit `classmethod` ¶

family `property` ¶

market `property` ¶

group_loadings `property` ¶

groups `property` ¶

group_labels `property` ¶

n_groups `property` ¶

n_factors `property` ¶

df `property` ¶

idiosyncratic `property` ¶

param_bounds `property` ¶

unfitted `classmethod` ¶

DynamicPath `dataclass` ¶

loglik `property` ¶

DynamicFitResult `dataclass` ¶

n_params `property` ¶

aic `property` ¶

bic `property` ¶

persistence `property` ¶

DccResult `dataclass` ¶

persistence `property` ¶

DiscreteFitResult `dataclass` ¶

n_params `property` ¶

aic `property` ¶

bic `property` ¶