rcopula

Copula modelling in Python — a full-featured replication of R's copula package, verified against it numerically.

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

Install

pip install rcopula

Quick start

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.

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

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".

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.

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).

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.

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.

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.

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.

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 has the theory, and rc.statarb.select_partners picks which assets to group together.

More finance tutorials

Four more in the same step-by-step style, each with a runnable script that checks its own numbers, are on the finance tutorials page:

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.

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-two scripts in examples/, each of which runs and asserts its own claims — so they cannot drift out of date without failing:

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.

There is a vine copula tutorial 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 (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 and CONTRIBUTING.md.

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