Metadata-Version: 2.5
Name: rcopula
Version: 0.2.0
Summary: Copula modelling in Python: a full-featured replication of R's `copula` package.
Project-URL: Homepage, https://github.com/Chrisebell24/python-copula
Project-URL: Documentation, https://chrisebell24.github.io/python-copula/
Project-URL: Repository, https://github.com/Chrisebell24/python-copula
Project-URL: Issues, https://github.com/Chrisebell24/python-copula/issues
Project-URL: Changelog, https://github.com/Chrisebell24/python-copula/blob/main/CHANGELOG.md
Author-email: Chris Bell <chris.e.bell24@gmail.com>
License: MIT License
        
        Copyright (c) 2026 Chris Bell
        
        Permission is hereby granted, free of charge, to any person obtaining a copy
        of this software and associated documentation files (the "Software"), to deal
        in the Software without restriction, including without limitation the rights
        to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
        copies of the Software, and to permit persons to whom the Software is
        furnished to do so, subject to the following conditions:
        
        The above copyright notice and this permission notice shall be included in all
        copies or substantial portions of the Software.
        
        THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
        IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
        FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
        AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
        LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
        OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
        SOFTWARE.
License-File: LICENSE
License-File: NOTICE
Keywords: archimedean,copula,dependence,risk,sklar,statistics,tail dependence,vine copula
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Financial and Insurance Industry
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Mathematics
Classifier: Typing :: Typed
Requires-Python: >=3.10
Requires-Dist: matplotlib>=3.7
Requires-Dist: numpy>=1.24
Requires-Dist: pandas>=2.0
Requires-Dist: scipy>=1.10
Provides-Extra: dev
Requires-Dist: mpmath>=1.3; extra == 'dev'
Requires-Dist: mypy>=1.11; extra == 'dev'
Requires-Dist: pytest-cov>=5.0; extra == 'dev'
Requires-Dist: pytest>=8.0; extra == 'dev'
Requires-Dist: ruff>=0.6; extra == 'dev'
Provides-Extra: docs
Requires-Dist: mkdocs-material>=9.5; extra == 'docs'
Requires-Dist: mkdocs>=1.6; extra == 'docs'
Requires-Dist: mkdocstrings[python]>=0.26; extra == 'docs'
Provides-Extra: garch
Requires-Dist: arch>=6.0; extra == 'garch'
Description-Content-Type: text/markdown

# rcopula

**Copula modelling in Python — a full-featured replication of R's [`copula`](https://cran.r-project.org/package=copula) package, verified against it numerically.**

[![CI](https://github.com/Chrisebell24/python-copula/actions/workflows/ci.yml/badge.svg)](https://github.com/Chrisebell24/python-copula/actions/workflows/ci.yml)
[![License: MIT](https://img.shields.io/badge/License-MIT-blue.svg)](https://github.com/Chrisebell24/python-copula/blob/main/LICENSE)

> **Status: pre-1.0, under active development.** The API may change before 1.0.

---

## Why this exists

R's `copula` package is the most complete copula toolkit in any language. Python's options
are fragmented, and roughly **60–70% of R `copula`'s statistical surface has no Python
equivalent at all**:

- no multiplier-bootstrap goodness-of-fit anywhere;
- no nested Archimedean copulas anywhere, because they need an exponentially
  tilted stable sampler that does not exist in NumPy, SciPy or any copula package;
- no Khoudraji device, general rotation wrapper, or arbitrary-copula mixtures;
- no Kendall distribution function `K`, no general Rosenblatt transform;
- **not one Python package returns standard errors for a fitted copula**;
- families missing under one API: Joe, AMH, Plackett, FGM, Marshall–Olkin, Galambos,
  Hüsler–Reiss, Tawn, t-EV, Fréchet bounds.

`rcopula` closes that gap, and proves it: every deterministic quantity is checked against
fixtures generated by running the real R package — **and separately against the
definitions themselves**, so agreeing with R cannot be mistaken for being right.

It also goes past R where the gap was worth closing:

- **time-varying copulas** — Patton (2006) and score-driven (GAS) recursions, plus
  Engle's DCC for a correlation matrix that moves;
- **discrete and mixed margins** — the exact inclusion-exclusion likelihood, with the
  identifiability caveat stated rather than buried;
- **vines** — R/C/D, in pure NumPy;
- **automatic family selection**, **JSON serialization**, and **bootstrap confidence
  intervals** with measured coverage.

## What copulas are for

Sklar's theorem says any joint distribution splits cleanly into **its margins** and **a
copula** that carries all the dependence. That separation is the whole point:

- **Model dependence and margins independently.** Fit each margin however you like — a
  GARCH residual distribution, a kernel density, a fitted Gamma — then fit the dependence
  structure separately.
- **Dependence measures that don't lie.** Kendall's τ and Spearman's ρ are invariant under
  monotone marginal transforms; Pearson correlation is not. Pearson also reports ~0 for
  strongly dependent pairs (e.g. `y = x²`).
- **Tail dependence.** Gaussian dependence says extreme joint events are essentially
  impossible. Real markets disagree. A t or Clayton copula prices joint crashes correctly;
  a Gaussian one does not. This distinction is not cosmetic — it changes VaR materially.

### Choosing a family

| If your data look like… | Use |
|---|---|
| Symmetric dependence, no joint extremes | `GaussianCopula` |
| Symmetric dependence **with** joint extremes in both tails | `StudentCopula` |
| Joint crashes but not joint booms (lower-tail dependence) | `ClaytonCopula` |
| Joint booms but not joint crashes (upper-tail dependence) | `GumbelCopula`, `JoeCopula` |
| Flexible middle, negative dependence allowed, no tail dependence | `FrankCopula` |
| Weak dependence only (τ ∈ [−0.18, 0.33]) | `AMHCopula` |
| Asymmetric tails | `KhoudrajiCopula`, `RotatedCopula` |
| Block / hierarchical group structure | `NestedArchimedean` |
| High-dimensional, pair-specific dependence | `VineCopula` |

## Install

```bash
pip install rcopula
```

## Quick start

```python
import rcopula as rc
from scipy import stats

# Construct and sample
cop = rc.ClaytonCopula(theta=2.0, dim=3)
u = cop.rvs(1000, random_state=42)

cop.pdf(u)  # density
cop.cdf(u)  # distribution function
cop.tau()  # Kendall's tau  -> 0.5
cop.lambda_()  # TailDependence(lower=0.7071, upper=0.0)

# Calibrate from a target dependence
rc.GumbelCopula.from_tau(0.5)  # theta = 2.0

# Fit, with standard errors
res = rc.fit(rc.GumbelCopula(dim=2), u[:, :2], method="mpl")
print(res.summary())

# Attach margins (Sklar, second direction)
mv = rc.CopulaDistribution(
    cop, margins=[stats.norm(1, 2), stats.expon(scale=1 / 3), stats.gamma(2)]
)
x = mv.rvs(1000, random_state=0)

# Goodness of fit
rc.gof_test(rc.ClaytonCopula(dim=3), x, simulation="mult")
```

## Tutorial: vine copulas for financial markets

A correlation matrix says how assets move together *on average*. It cannot say
that small caps crash with the S&P far more often than they rally with it, or
that Treasuries rally in a sell-off. A **vine copula** can: it links assets in
pairs, and every pair gets its own shape. Eight steps, from a table of returns to
a crash-aware portfolio. The full runnable version, with the data simulated so it
works offline, is
[`examples/28_vine_markets_howto.py`](https://github.com/Chrisebell24/python-copula/blob/main/examples/28_vine_markets_howto.py).

**1. Get daily returns.** Any `DataFrame` of returns, one column per asset.

```python
import numpy as np
import yfinance as yf  # pip install yfinance -- or use your own data

tickers = ["SPY", "QQQ", "IWM", "HYG", "TLT"]  # large caps, tech, small caps, credit, bonds
prices = yf.download(tickers, start="2015-01-01")["Close"][tickers]
returns = np.log(prices).diff().dropna()
```

**2. Strip out volatility.** Calm and stormy periods hit every asset at once, and
a copula fitted to raw returns mistakes that for dependence. Fit a GARCH to each
asset and keep the *standardised residuals*: how surprising each day was, given
how volatile things already were. Then convert them to ranks in (0, 1), where
0.01 means "a 1-in-100 bad day for this asset".

```python
import pandas as pd
import rcopula as rc
from rcopula.garch import fit_garch

margins = [fit_garch(returns[t], dist="t", name=t) for t in tickers]
u = rc.pseudo_obs(pd.DataFrame({m.name: m.resid for m in margins}))
```

**3. Fit the vine.** One line. A C-vine puts the most connected asset at the
centre and links every other asset to it; higher trees capture what is left once
that asset is accounted for. Each link picks its own family by AIC.

```python
vine = rc.fit_vine(u, structure="C")
print(vine.describe())
# C-vine copula, dim 5, order [0, 1, 2, 3, 4]       <- SPY (column 0) is the centre
#   tree 1  0,1      Student copula, rho=0.877, df=4.06
#   tree 1  0,2      Clayton copula, theta=2.59
#   tree 1  0,3      Clayton copula, theta=1.32
#   tree 1  0,4      Frank copula, theta=-1.94
#   tree 2  1,2|0    Gaussian copula, rho=0.208      <- QQQ-IWM once SPY is known
#   ...
```

**4. Read what it found.** Tree 1 is the part to read. Kendall's τ is the overall
co-movement (−1 to 1); tail dependence is the chance that two assets have an
extreme day *together* (0 to 1).

```python
for i, pair in enumerate(vine.pair_copulas[0]):
    asset = tickers[vine.order[i + 1]]
    tail = pair.lambda_()
    print(
        f"SPY-{asset}  {pair.name:<8} tau={pair.tau():+.2f}  "
        f"crash together={tail.lower:.2f}  rally together={tail.upper:.2f}"
    )
# SPY-QQQ  Student  tau=+0.68  crash together=0.59  rally together=0.59
# SPY-IWM  Clayton  tau=+0.56  crash together=0.77  rally together=0.00
# SPY-HYG  Clayton  tau=+0.40  crash together=0.59  rally together=0.00
# SPY-TLT  Frank    tau=-0.21  crash together=0.00  rally together=0.00
```

Small caps and credit fall with the S&P but do not rally with it; bonds move the
other way. A Gaussian copula, which is what a correlation matrix implies, sets
every one of those tail numbers to zero.

**5. Simulate tomorrow and measure risk.** Put the GARCH margins and the vine
back together. The forecast starts from *today's* volatility, so the risk number
moves with the market.

```python
from rcopula.garch import CopulaGarch

weights = [0.30, 0.15, 0.15, 0.10, 0.30]
model = CopulaGarch(margins, vine, innovations="parametric")
model.forecast_risk(weights, alpha=0.99, n=200_000)
# {'var': 0.0203, 'expected_shortfall': 0.0271, ...}   Gaussian copula: 0.0186 and 0.0240
```

**6. Ask how often they crash together.** Simulate from the vine and count the
days when SPY, QQQ and IWM all have a 1-in-50 bad day at once.

```python
sims = vine.rvs(400_000)
np.mean(np.all(sims[:, :3] < 0.02, axis=1))
# vine 1.13%   data 1.04%   Gaussian copula 0.48%   independent 0.0008%
```

**7. Stress test.** Keep only the simulated days when SPY has a 1-in-100 loss,
and look at everything else on those days.

```python
crash_days = sims[sims[:, 0] < 0.01]
(crash_days < 0.05).mean(axis=0)  # chance each asset has its own worst-5% day
# QQQ 93%   IWM 99%   HYG 92%   TLT 2%
```

**8. Build a crash-aware portfolio.** Feed simulated days to an optimiser that
minimises expected shortfall (the average loss on the worst 1% of days), capped at
40% per fund. Because the scenarios come from the vine, the optimiser knows which
assets fail together.

```python
from rcopula.portfolio import mean_cvar_weights

scenarios = model.forecast(horizon=1, n=5_000)
w = mean_cvar_weights(scenarios, alpha=0.99, bounds=(0.0, 0.40))
```

Optimised on Gaussian-copula scenarios instead, the portfolio holds less of the
TLT hedge, and its real bad-day loss comes out a little higher. That gain is
modest. Most of what a vine adds is in steps 5–7: knowing *before* the bad day
that the risk number is too low, and by how much.

**Which vine?** Use a **C-vine** when one asset or factor drives the rest (an
index, a sector ETF, oil for energy names). Use a **D-vine** (`structure="D"`) for a
chain, such as points on a yield curve or futures expiries. For 20+ assets, add
`truncate=2` to model only the first two trees and treat the rest as independent.
The [vine tutorial](https://chrisebell24.github.io/python-copula/vines/) has the
theory, and `rc.statarb.select_partners` picks which assets to group together.

### More finance tutorials

Three more in the same step-by-step style, each with a runnable script that checks
its own numbers, are on the [finance tutorials](https://chrisebell24.github.io/python-copula/finance-tutorials/) page:

| | |
|---|---|
| [Risk management](https://chrisebell24.github.io/python-copula/finance-tutorials/#risk-management-for-a-trading-book) | VaR and expected shortfall in dollars, a VaR backtest, which desk the risk comes from, how much diversification really buys, contagion between desks, reverse stress testing — [`examples/29`](https://github.com/Chrisebell24/python-copula/blob/main/examples/29_risk_management_howto.py) |
| [Trading strategies](https://chrisebell24.github.io/python-copula/finance-tutorials/#trading-strategies-with-copulas) | Picking pairs, a copula mispricing signal, an out-of-sample backtest with costs against the z-score trade, a vine basket trade, and an alarm for when a pair breaks — [`examples/30`](https://github.com/Chrisebell24/python-copula/blob/main/examples/30_trading_strategies_howto.py) |
| [Valuing odd assets](https://chrisebell24.github.io/python-copula/finance-tutorials/#valuing-odd-assets) | Worst-of notes, basket puts, first-to-default baskets, catastrophe bonds and private stakes, priced at the same correlation under different tails — [`examples/31`](https://github.com/Chrisebell24/python-copula/blob/main/examples/31_valuing_odd_assets_howto.py) |

## Datasets

Nothing is bundled: R `copula`'s datasets are GPL-3 and this package is MIT, so
copying them in would be a licence conflict. Instead the loader ships URLs,
licences and SHA-256 digests, and fetches permissively-licensed sources on first
use.

```python
flood = rc.datasets.load("nwis_peaks")  # USGS, public domain
wine = rc.datasets.load("uci_wine_quality_red")  # UCI, CC BY 4.0
flood.attrs["citation"]
```

Five sources: USGS annual peak flows and NOAA GHCN-Daily (both live endpoints, so
unpinnable), and three static UCI tables — abalone, red wine quality and NASA
aerofoil self-noise — whose SHA-256 is checked on every fetch. A test asserts no
data file is ever committed.

## Examples

Thirty-one scripts in [`examples/`](https://github.com/Chrisebell24/python-copula/tree/main/examples/), each of which **runs and asserts
its own claims** — so they cannot drift out of date without failing:

```bash
python examples/04_tail_dependence.py     # what tau does not tell you
python examples/13_hydrology_return_periods.py
pytest tests/test_examples.py             # run them all
```

They cover the R book chapters; applied finance (CDO tranches, basket and CMS
spread options, pairs trading, copula-GARCH, operational risk); hydrology,
climate compound extremes, structural reliability and survival analysis; the
statistics and machine-learning uses (mutual information, robust covariance,
synthetic data, conditional independence); vines; the nine diagnostic plots;
time-varying dependence; count data; quasi-random sampling; Archimedean
internals at d = 100; nested trees; and what rounding does to a rank statistic.

## Coming from R?

The API is Pythonic (scipy/statsmodels verbs), but R's option strings are preserved so you
can paste them across.

| R | `rcopula` |
|---|---|
| `claytonCopula(2, dim = 3)` | `rc.ClaytonCopula(theta=2, dim=3)` |
| `rCopula(n, cop)` | `cop.rvs(n)` |
| `dCopula(u, cop)` / `pCopula(u, cop)` | `cop.pdf(u)` / `cop.cdf(u)` |
| `tau(cop)` / `rho(cop)` / `lambda(cop)` | `cop.tau()` / `cop.rho()` / `cop.lambda_()` |
| `iTau(claytonCopula(), 0.5)` | `rc.ClaytonCopula.from_tau(0.5)` |
| `fitCopula(cop, u, method = "mpl")` | `rc.fit(cop, u, method="mpl")` |
| `pobs(x)` | `rc.pseudo_obs(x)` |
| `mvdc(cop, margins, paramMargins)` | `rc.CopulaDistribution(cop, margins=[...])` |
| `fitMvdc(x, mvdc, start)` | `rc.fit_joint(mv, x)` — margins and copula together |
| `cCopula(u, cop)` | `rc.rosenblatt(cop, u)` |
| `gofCopula(cop, x, simulation = "mult")` | `rc.gof_test(cop, x, simulation="mult")` |
| `gofT2stat(u1, u2)` | `rc.gof_two_sample(x, y)` — did the dependence change? |
| `dependogram(x, d)` | `rc.dependogram(x)` — which subsets are dependent |
| `rotCopula(cop, flip)` | `rc.RotatedCopula(cop, flip)` / `rc.survival(cop)` |
| `khoudrajiCopula(c1, c2, shapes)` | `rc.KhoudrajiCopula(c1, c2, shapes)` |
| `mixCopula(list(c1, c2), w)` | `rc.MixtureCopula([c1, c2], w)` |
| `opower(cop, alpha)` | `rc.opower(cop, alpha)` — a second dial on the tails |
| `margCopula(cop, keep)` | `rc.marginal_copula(cop, [0, 3])` |
| `fitLambda(x)` | `rc.fit_lambda(x)` — tail dependence without a family |
| `onacopula("G", C(1.5, , list(C(4, 1:3))))` | `rc.NestedArchimedean(rc.GumbelCopula(1.5), ...)` |
| `enacopula(u, cop, method = "etau")` | `rc.fit_nested(structure, u)` |
| *(VineCopula::RVineStructureSelect)* | `rc.fit_vine(u, structure="D")` |
| `retstable(n, V0, h, alpha)` | `rcopula.special.stable.retstable(...)` |
| `pK(t, cop, d)` / `qK` / `dK` / `rK` | `rc.kendall_cdf/ppf/pdf/rvs(cop, ...)` |
| `Kn(u, x)` | `rc.kendall_empirical(x, u)` |
| `xvCopula(cop, x, k)` | `rc.cross_validate(cop, x, k)` |
| `contour(cop, dCopula)` / `persp(...)` | `rc.plots.contour(cop)` / `rc.plots.surface(cop)` |
| `pairsRosenblatt(...)` | `rc.plots.pairs_rosenblatt(cop, u)` — *where* a fit fails |
| *(no equivalent)* | `rc.select_copula(u)` — rank every admissible family |
| *(no equivalent)* | `rc.plots.tail_concentration(u, [c1, c2])` — where families disagree |
| *(no equivalent)* | `rc.dynamic.fit_dynamic(u, cop)` — Patton / GAS time-varying dependence |
| *(no equivalent)* | `rc.dynamic.fit_dcc(u)` — Engle's dynamic conditional correlation |
| *(no equivalent)* | `rc.discrete.fit_discrete(x, cop, margins)` — exact likelihood for counts |
| *(no equivalent)* | `rc.bootstrap.bootstrap_measure(u, "tau")` — BCa intervals, `n_jobs=` |
| *(no equivalent)* | `rc.serialize.to_json(cop)` — a model file you can read and diff |
| *(no equivalent)* | `rc.statarb.select_pairs(returns)` — rank every pair, six criteria |
| *(no equivalent)* | `rc.statarb.select_partners(returns, "AAPL")` — vine quadruples |
| `rtrafo` / inverse `cCopula` | `rc.inverse_rosenblatt(cop, z)` |
| `rAntitheticVariates` | `rc.sampling.antithetic_rvs(cop, n)` |
| `rLatinHypercube` | `rc.sampling.latin_hypercube_rvs(cop, n)` |
| *(qrng vignette)* | `rc.sampling.quasi_rvs(cop, n)` — Sobol/Halton through the copula |
| `htrafo(u, cop)` | `rc.htrafo(cop, u)` |

There is a [vine copula tutorial](https://chrisebell24.github.io/python-copula/vines/)
that reproduces Aas et al. (2009) equation by equation — the D-vine and C-vine
factorisations agree with the printed formulas to 2e-15, and the paper's count of
"12 D-vines and 12 C-vines" in four dimensions is reproduced by enumeration.

The full translation table, an API reference, and the notes on where R is the
less accurate side are in the [documentation](https://chrisebell24.github.io/python-copula/)
(`make docs` to read it locally).

## How parity is verified

R's RNG cannot be matched bit-for-bit by NumPy, so verification splits in two:

- **Deterministic quantities** — densities, CDFs, generators, Pickands functions, τ/ρ/λ,
  `iTau`/`iRho`, Rosenblatt transforms, the Kendall function, and every GoF *statistic* —
  must match R to ~1e-10.
- **RNG-dependent quantities** — R's random draws are **exported to CSV and fed into
  Python**, making everything downstream deterministic. Fitted parameters must match to
  ~1e-6 and standard errors to ~1e-5. Only bootstrap p-values remain stochastic, and those
  are checked statistically.

Fixtures live in `tests/golden/` and are committed, so CI needs no R. Regenerate with
`make golden`. A published parity report lists every replicated function and its maximum
observed deviation.

## License and relationship to R `copula`

`rcopula` is MIT. R's `copula` is GPL-3. **No R source code was translated into this
project** — every algorithm is implemented from its originating published paper, and the R
package is used solely as a black-box test oracle. See [`NOTICE`](https://github.com/Chrisebell24/python-copula/blob/main/NOTICE) and
[`CONTRIBUTING.md`](https://github.com/Chrisebell24/python-copula/blob/main/CONTRIBUTING.md).

`rcopula` is not affiliated with or endorsed by the authors of the R `copula` package.
