Copulas for Multi-Asset Modeling and Portfolio Risk

TL;DR – Copulas model dependence between assets beyond simple correlation, capturing tail dependence crucial for risk management. This guide covers Gaussian, t-copulas, Archimedean families, and portfolio simulation with production code.

1. Why Copulas?#

Correlation fails to capture:

Tail dependence – Assets crash together more than correlation suggests
Non-linear dependence – Asymmetric relationships
Non-normal margins – Fat tails, skewness

Copulas separate:

Marginal distributions – Individual asset behavior
Dependence structure – How assets move together

Sklar's Theorem: Any joint distribution can be decomposed into margins and a copula.

2. Copula Basics #

2.1 Definition #

A copula $C$ is a multivariate CDF with uniform margins: $C(u_1, \ldots, u_n) = P(U_1 \leq u_1, \ldots, U_n \leq u_n)$

where $U_i \sim \text{Uniform}(0, 1)$ .

Key property: For any joint distribution $F$ with margins $F_1, \ldots, F_n$ : $F(x_1, \ldots, x_n) = C(F_1(x_1), \ldots, F_n(x_n))$

2.2 Practical Workflow #

Fit marginal distributions to each asset
Transform to uniform via probability integral transform
Fit copula to uniform data
Simulate from copula and transform back

python

1import numpy as np
2from scipy import stats
3import matplotlib.pyplot as plt
4
5def fit_marginals(returns):
6    """
7    Fit marginal distributions to each asset.
8    
9    Args:
10        returns: Array of shape (n_samples, n_assets)
11        
12    Returns:
13        List of fitted distributions
14    """
15    n_assets = returns.shape[1]
16    marginals = []
17    
18    for i in range(n_assets):
19        # Fit t-distribution (captures fat tails)
20        params = stats.t.fit(returns[:, i])
21        marginals.append(('t', params))
22    
23    return marginals
24
25def transform_to_uniform(returns, marginals):
26    """Transform returns to uniform via CDF."""
27    n_samples, n_assets = returns.shape
28    U = np.zeros_like(returns)
29    
30    for i in range(n_assets):
31        dist_name, params = marginals[i]
32        if dist_name == 't':
33            U[:, i] = stats.t.cdf(returns[:, i], *params)
34        elif dist_name == 'norm':
35            U[:, i] = stats.norm.cdf(returns[:, i], *params)
36    
37    return U
38
39# Example
40np.random.seed(42)
41n_samples, n_assets = 1000, 3
42
43# Generate correlated returns
44mean = [0.05, 0.07, 0.06]
45cov = [[0.04, 0.02, 0.01],
46       [0.02, 0.06, 0.015],
47       [0.01, 0.015, 0.05]]
48returns = np.random.multivariate_normal(mean, cov, n_samples)
49
50# Fit and transform
51marginals = fit_marginals(returns)
52U = transform_to_uniform(returns, marginals)
53
54print("Uniform data shape:", U.shape)
55print("Uniform data range:", U.min(), "to", U.max())
56

3. Gaussian Copula #

3.1 Definition #

The Gaussian copula with correlation matrix $\Sigma$ : $C^{Gauss}(u_1, \ldots, u_n; \Sigma) = \Phi_\Sigma(\Phi^{-1}(u_1), \ldots, \Phi^{-1}(u_n))$

where:

$\Phi$ = standard normal CDF
$\Phi_\Sigma$ = multivariate normal CDF with correlation $\Sigma$

Properties:

No tail dependence (independence in extremes)
Easy to fit and simulate
Used in CDO pricing (infamous in 2008 crisis)

3.2 Implementation #

python

1class GaussianCopula:
2    """Gaussian copula for multivariate dependence."""
3    
4    def __init__(self, correlation_matrix):
5        self.corr = np.array(correlation_matrix)
6        self.n_dim = len(correlation_matrix)
7    
8    def fit(self, U):
9        """
10        Fit correlation matrix from uniform data.
11        
12        Args:
13            U: Uniform data of shape (n_samples, n_dim)
14        """
15        # Transform to normal
16        Z = stats.norm.ppf(U)
17        
18        # Compute correlation
19        self.corr = np.corrcoef(Z.T)
20        
21        return self
22    
23    def sample(self, n_samples):
24        """
25        Generate samples from Gaussian copula.
26        
27        Returns:
28            Uniform samples of shape (n_samples, n_dim)
29        """
30        # Generate correlated normals
31        Z = np.random.multivariate_normal(
32            np.zeros(self.n_dim),
33            self.corr,
34            n_samples
35        )
36        
37        # Transform to uniform
38        U = stats.norm.cdf(Z)
39        
40        return U
41    
42    def pdf(self, U):
43        """Copula density."""
44        Z = stats.norm.ppf(U)
45        
46        # Multivariate normal density / product of marginal densities
47        det_corr = np.linalg.det(self.corr)
48        inv_corr = np.linalg.inv(self.corr)
49        
50        exponent = -0.5 * np.sum(Z @ (inv_corr - np.eye(self.n_dim)) * Z, axis=1)
51        
52        return np.exp(exponent) / np.sqrt(det_corr)
53
54# Example usage
55copula = GaussianCopula(correlation_matrix=[[1, 0.5], [0.5, 1]])
56U_samples = copula.sample(n_samples=1000)
57
58# Visualize
59plt.figure(figsize=(10, 5))
60plt.subplot(1, 2, 1)
61plt.scatter(U_samples[:, 0], U_samples[:, 1], alpha=0.5, s=10)
62plt.xlabel('U1')
63plt.ylabel('U2')
64plt.title('Gaussian Copula Samples (rho=0.5)')
65plt.grid(True, alpha=0.3)
66
67plt.subplot(1, 2, 2)
68plt.hist2d(U_samples[:, 0], U_samples[:, 1], bins=30, cmap='Blues')
69plt.xlabel('U1')
70plt.ylabel('U2')
71plt.title('Density')
72plt.colorbar()
73plt.tight_layout()
74plt.show()
75

4. Student-t Copula #

4.1 Motivation #

Key advantage: Exhibits tail dependence – assets crash together.

Student-t copula with correlation $\Sigma$ and degrees of freedom $\nu$ : $C^t(u_1, \ldots, u_n; \Sigma, \nu) = t_{\Sigma, \nu}(t_\nu^{-1}(u_1), \ldots, t_\nu^{-1}(u_n))$

Tail dependence coefficient: $\lambda = 2 \bar{t}_{\nu+1}\left(-\sqrt{\frac{(\nu+1)(1-\rho)}{1+\rho}}\right)$

where $\rho$ is correlation and $\bar{t}$ is survival function.

4.2 Implementation #

python

1class StudentTCopula:
2    """Student-t copula with tail dependence."""
3    
4    def __init__(self, correlation_matrix, df):
5        self.corr = np.array(correlation_matrix)
6        self.df = df  # Degrees of freedom
7        self.n_dim = len(correlation_matrix)
8    
9    def fit(self, U, method='mle'):
10        """
11        Fit correlation and degrees of freedom.
12        
13        Uses maximum likelihood estimation.
14        """
15        from scipy.optimize import minimize
16        
17        def neg_log_likelihood(params):
18            df = params[0]
19            # Transform to t-distributed
20            X = stats.t.ppf(U, df)
21            
22            # Compute correlation
23            corr = np.corrcoef(X.T)
24            
25            # Log-likelihood (simplified)
26            try:
27                mvt = stats.multivariate_t(loc=np.zeros(self.n_dim),
28                                          shape=corr, df=df)
29                return -np.sum(mvt.logpdf(X))
30            except:
31                return 1e10
32        
33        # Optimize
34        result = minimize(neg_log_likelihood, x0=[10], bounds=[(2.1, 30)])
35        self.df = result.x[0]
36        
37        # Recompute correlation with optimal df
38        X = stats.t.ppf(U, self.df)
39        self.corr = np.corrcoef(X.T)
40        
41        return self
42    
43    def sample(self, n_samples):
44        """Generate samples from t-copula."""
45        # Generate from multivariate t
46        from scipy.stats import multivariate_t
47        
48        mvt = multivariate_t(loc=np.zeros(self.n_dim),
49                            shape=self.corr, df=self.df)
50        X = mvt.rvs(n_samples)
51        
52        # Transform to uniform
53        U = stats.t.cdf(X, self.df)
54        
55        return U
56    
57    def tail_dependence(self, rho):
58        """
59        Compute tail dependence coefficient.
60        
61        Args:
62            rho: Correlation between two assets
63        """
64        from scipy.stats import t as t_dist
65        
66        arg = -np.sqrt((self.df + 1) * (1 - rho) / (1 + rho))
67        lambda_tail = 2 * t_dist.sf(arg, self.df + 1)
68        
69        return lambda_tail
70
71# Example
72t_copula = StudentTCopula(correlation_matrix=[[1, 0.5], [0.5, 1]], df=5)
73U_t = t_copula.sample(n_samples=1000)
74
75# Compare tail behavior
76plt.figure(figsize=(14, 5))
77
78plt.subplot(1, 3, 1)
79plt.scatter(U_samples[:, 0], U_samples[:, 1], alpha=0.3, s=10, label='Gaussian')
80plt.xlabel('U1')
81plt.ylabel('U2')
82plt.title('Gaussian Copula')
83plt.grid(True, alpha=0.3)
84
85plt.subplot(1, 3, 2)
86plt.scatter(U_t[:, 0], U_t[:, 1], alpha=0.3, s=10, label='Student-t', color='red')
87plt.xlabel('U1')
88plt.ylabel('U2')
89plt.title('Student-t Copula (df=5)')
90plt.grid(True, alpha=0.3)
91
92plt.subplot(1, 3, 3)
93# Focus on lower tail
94mask_gauss = (U_samples[:, 0] < 0.1) & (U_samples[:, 1] < 0.1)
95mask_t = (U_t[:, 0] < 0.1) & (U_t[:, 1] < 0.1)
96plt.scatter(U_samples[mask_gauss, 0], U_samples[mask_gauss, 1], 
97           alpha=0.5, s=20, label='Gaussian')
98plt.scatter(U_t[mask_t, 0], U_t[mask_t, 1], 
99           alpha=0.5, s=20, label='Student-t', color='red')
100plt.xlabel('U1')
101plt.ylabel('U2')
102plt.title('Lower Tail (U1, U2 < 0.1)')
103plt.legend()
104plt.grid(True, alpha=0.3)
105
106plt.tight_layout()
107plt.show()
108
109print("Tail dependence (t-copula):", t_copula.tail_dependence(0.5))
110

5. Archimedean Copulas #

5.1 Families #

Archimedean copulas are defined via generator function $\phi$ : $C(u_1, u_2) = \phi^{-1}(\phi(u_1) + \phi(u_2))$

Clayton Copula (lower tail dependence): $C(u_1, u_2) = (u_1^{-\theta} + u_2^{-\theta} - 1)^{-1/\theta}$

Gumbel Copula (upper tail dependence): $C(u_1, u_2) = \exp\left(-\left[(-\log u_1)^\theta + (-\log u_2)^\theta\right]^{1/\theta}\right)$

Frank Copula (symmetric, no tail dependence): $C(u_1, u_2) = -\frac{1}{\theta} \log\left(1 + \frac{(e^{-\theta u_1} - 1)(e^{-\theta u_2} - 1)}{e^{-\theta} - 1}\right)$

5.2 Implementation #

python

1class ClaytonCopula:
2    """Clayton copula with lower tail dependence."""
3    
4    def __init__(self, theta):
5        self.theta = theta  # theta > 0
6    
7    def cdf(self, u1, u2):
8        """Copula CDF."""
9        return (u1**(-self.theta) + u2**(-self.theta) - 1)**(-1/self.theta)
10    
11    def sample(self, n_samples):
12        """Generate samples using conditional sampling."""
13        U1 = np.random.uniform(0, 1, n_samples)
14        V = np.random.uniform(0, 1, n_samples)
15        
16        # Conditional distribution
17        U2 = (U1**(-self.theta) * (V**(-self.theta/(1+self.theta)) - 1) + 1)**(-1/self.theta)
18        
19        return np.column_stack([U1, U2])
20    
21    def tail_dependence_lower(self):
22        """Lower tail dependence coefficient."""
23        return 2**(-1/self.theta)
24
25class GumbelCopula:
26    """Gumbel copula with upper tail dependence."""
27    
28    def __init__(self, theta):
29        self.theta = theta  # theta >= 1
30    
31    def cdf(self, u1, u2):
32        """Copula CDF."""
33        return np.exp(-((-np.log(u1))**self.theta + 
34                        (-np.log(u2))**self.theta)**(1/self.theta))
35    
36    def sample(self, n_samples):
37        """Generate samples."""
38        # Use Laplace transform method
39        from scipy.stats import levy_stable
40        
41        # Generate from stable distribution
42        alpha = 1 / self.theta
43        stable = levy_stable(alpha=alpha, beta=1)
44        V = stable.rvs(n_samples)
45        
46        E1 = np.random.exponential(1, n_samples)
47        E2 = np.random.exponential(1, n_samples)
48        
49        U1 = np.exp(-(E1 / V)**(1/self.theta))
50        U2 = np.exp(-(E2 / V)**(1/self.theta))
51        
52        return np.column_stack([U1, U2])
53    
54    def tail_dependence_upper(self):
55        """Upper tail dependence coefficient."""
56        return 2 - 2**(1/self.theta)
57
58# Compare copulas
59clayton = ClaytonCopula(theta=2)
60gumbel = GumbelCopula(theta=2)
61
62U_clayton = clayton.sample(1000)
63U_gumbel = gumbel.sample(1000)
64
65plt.figure(figsize=(14, 5))
66
67plt.subplot(1, 3, 1)
68plt.scatter(U_clayton[:, 0], U_clayton[:, 1], alpha=0.3, s=10)
69plt.title('Clayton (Lower Tail Dep)')
70plt.xlabel('U1')
71plt.ylabel('U2')
72plt.grid(True, alpha=0.3)
73
74plt.subplot(1, 3, 2)
75plt.scatter(U_gumbel[:, 0], U_gumbel[:, 1], alpha=0.3, s=10, color='red')
76plt.title('Gumbel (Upper Tail Dep)')
77plt.xlabel('U1')
78plt.ylabel('U2')
79plt.grid(True, alpha=0.3)
80
81plt.subplot(1, 3, 3)
82print("Clayton lower tail dep:", clayton.tail_dependence_lower())
83print("Gumbel upper tail dep:", gumbel.tail_dependence_upper())
84plt.text(0.1, 0.5, 'Clayton lower tail:\n{:.3f}'.format(clayton.tail_dependence_lower()),
85         fontsize=12)
86plt.text(0.1, 0.3, 'Gumbel upper tail:\n{:.3f}'.format(gumbel.tail_dependence_upper()),
87         fontsize=12, color='red')
88plt.xlim(0, 1)
89plt.ylim(0, 1)
90plt.title('Tail Dependence Comparison')
91plt.axis('off')
92
93plt.tight_layout()
94plt.show()
95

6. Portfolio Risk Simulation #

6.1 Complete Workflow #

python

1def simulate_portfolio_returns(
2    marginals,
3    copula,
4    weights,
5    n_scenarios=10000
6):
7    """
8    Simulate portfolio returns using copula.
9    
10    Args:
11        marginals: List of (dist_name, params) for each asset
12        copula: Fitted copula object
13        weights: Portfolio weights
14        n_scenarios: Number of scenarios
15        
16    Returns:
17        Portfolio returns
18    """
19    # Sample from copula
20    U = copula.sample(n_scenarios)
21    
22    # Transform to asset returns
23    n_assets = len(marginals)
24    returns = np.zeros((n_scenarios, n_assets))
25    
26    for i in range(n_assets):
27        dist_name, params = marginals[i]
28        if dist_name == 't':
29            returns[:, i] = stats.t.ppf(U[:, i], *params)
30        elif dist_name == 'norm':
31            returns[:, i] = stats.norm.ppf(U[:, i], *params)
32    
33    # Compute portfolio returns
34    portfolio_returns = returns @ weights
35    
36    return portfolio_returns, returns
37
38# Example: 3-asset portfolio
39weights = np.array([0.4, 0.3, 0.3])
40
41# Fit copula to historical data
42t_copula_3d = StudentTCopula(
43    correlation_matrix=[[1, 0.6, 0.4],
44                       [0.6, 1, 0.5],
45                       [0.4, 0.5, 1]],
46    df=6
47)
48
49# Simulate
50port_returns, asset_returns = simulate_portfolio_returns(
51    marginals,
52    t_copula_3d,
53    weights,
54    n_scenarios=10000
55)
56
57# Risk metrics
58VaR_95 = np.percentile(port_returns, 5)
59CVaR_95 = np.mean(port_returns[port_returns <= VaR_95])
60
61print("Portfolio VaR (95%): {:.2%}".format(VaR_95))
62print("Portfolio CVaR (95%): {:.2%}".format(CVaR_95))
63
64# Plot distribution
65plt.figure(figsize=(12, 5))
66
67plt.subplot(1, 2, 1)
68plt.hist(port_returns, bins=50, alpha=0.7, density=True)
69plt.axvline(VaR_95, color='r', linestyle='--', label='VaR 95%')
70plt.axvline(CVaR_95, color='darkred', linestyle='--', label='CVaR 95%')
71plt.xlabel('Portfolio Return')
72plt.ylabel('Density')
73plt.title('Portfolio Return Distribution')
74plt.legend()
75plt.grid(True, alpha=0.3)
76
77plt.subplot(1, 2, 2)
78plt.scatter(asset_returns[:, 0], asset_returns[:, 1], alpha=0.1, s=1)
79plt.xlabel('Asset 1 Return')
80plt.ylabel('Asset 2 Return')
81plt.title('Joint Return Distribution')
82plt.grid(True, alpha=0.3)
83
84plt.tight_layout()
85plt.show()
86

6.2 Stress Testing #

python

1def stress_test_portfolio(marginals, copula, weights, stress_scenarios):
2    """
3    Stress test portfolio under extreme scenarios.
4    
5    Args:
6        stress_scenarios: List of (asset_idx, percentile) tuples
7    """
8    results = []
9    
10    for scenario_name, conditions in stress_scenarios.items():
11        # Generate scenarios
12        U = copula.sample(10000)
13        returns = np.zeros((10000, len(marginals)))
14        
15        for i in range(len(marginals)):
16            dist_name, params = marginals[i]
17            returns[:, i] = stats.t.ppf(U[:, i], *params)
18        
19        # Apply stress conditions
20        mask = np.ones(10000, dtype=bool)
21        for asset_idx, percentile in conditions:
22            threshold = np.percentile(returns[:, asset_idx], percentile)
23            mask &= (returns[:, asset_idx] <= threshold)
24        
25        # Compute stressed portfolio returns
26        stressed_returns = returns[mask] @ weights
27        
28        results.append({
29            'scenario': scenario_name,
30            'mean': np.mean(stressed_returns),
31            'VaR_95': np.percentile(stressed_returns, 5),
32            'worst': np.min(stressed_returns)
33        })
34    
35    return results
36
37# Define stress scenarios
38stress_scenarios = {
39    'Market Crash': [(0, 5), (1, 5), (2, 5)],  # All assets in bottom 5%
40    'Asset 1 Crash': [(0, 1)],  # Asset 1 in bottom 1%
41    'Correlation Breakdown': [(0, 5), (1, 95)]  # Asset 1 down, Asset 2 up
42}
43
44stress_results = stress_test_portfolio(marginals, t_copula_3d, weights, stress_scenarios)
45
46for result in stress_results:
47    print("Scenario: {}".format(result['scenario']))
48    print("  Mean return: {:.2%}".format(result['mean']))
49    print("  VaR 95%: {:.2%}".format(result['VaR_95']))
50    print("  Worst case: {:.2%}".format(result['worst']))
51    print()
52

7. Production Checklist #

Validate copula choice – Test goodness-of-fit
Check tail dependence – Measure empirically
Use robust estimation – Handle outliers
Implement model selection – AIC/BIC comparison
Cache fitted copulas – Expensive to refit
Monitor correlation stability – Rolling window
Handle high dimensions – Vine copulas for 10+ assets
Validate simulations – Compare to historical
Document assumptions – Stationarity, ergodicity
Implement backtesting – VaR exceedances

References #

Nelsen, R. (2006). An Introduction to Copulas. Springer.
McNeil, A., Frey, R., & Embrechts, P. (2005). Quantitative Risk Management. Princeton.
Joe, H. (2014). Dependence Modeling with Copulas. CRC Press.
Cherubini, U., Luciano, E., & Vecchiato, W. (2004). Copula Methods in Finance. Wiley.

Copulas are essential for multi-asset risk management. Use Gaussian for simplicity, Student-t for tail dependence, and Archimedean for specific tail behavior. Always validate with historical stress tests.

TL;DR – Copulas model dependence between assets beyond simple correlation, capturing tail dependence crucial for risk management. This guide covers Gaussian, t-copulas, Archimedean families, and portfolio simulation with production code.

1. Why Copulas?#

Correlation fails to capture:

Tail dependence – Assets crash together more than correlation suggests
Non-linear dependence – Asymmetric relationships
Non-normal margins – Fat tails, skewness

Copulas separate:

Marginal distributions – Individual asset behavior
Dependence structure – How assets move together

Sklar's Theorem: Any joint distribution can be decomposed into margins and a copula.

2. Copula Basics #

2.1 Definition #

A copula $C$ is a multivariate CDF with uniform margins: $C(u_1, \ldots, u_n) = P(U_1 \leq u_1, \ldots, U_n \leq u_n)$

where $U_i \sim \text{Uniform}(0, 1)$ .

Key property: For any joint distribution $F$ with margins $F_1, \ldots, F_n$ : $F(x_1, \ldots, x_n) = C(F_1(x_1), \ldots, F_n(x_n))$

2.2 Practical Workflow #

Fit marginal distributions to each asset
Transform to uniform via probability integral transform
Fit copula to uniform data
Simulate from copula and transform back

python

1import numpy as np
2from scipy import stats
3import matplotlib.pyplot as plt
4
5def fit_marginals(returns):
6    """
7    Fit marginal distributions to each asset.
8    
9    Args:
10        returns: Array of shape (n_samples, n_assets)
11        
12    Returns:
13        List of fitted distributions
14    """
15    n_assets = returns.shape[1]
16    marginals = []
17    
18    for i in range(n_assets):
19        # Fit t-distribution (captures fat tails)
20        params = stats.t.fit(returns[:, i])
21        marginals.append(('t', params))
22    
23    return marginals
24
25def transform_to_uniform(returns, marginals):
26    """Transform returns to uniform via CDF."""
27    n_samples, n_assets = returns.shape
28    U = np.zeros_like(returns)
29    
30    for i in range(n_assets):
31        dist_name, params = marginals[i]
32        if dist_name == 't':
33            U[:, i] = stats.t.cdf(returns[:, i], *params)
34        elif dist_name == 'norm':
35            U[:, i] = stats.norm.cdf(returns[:, i], *params)
36    
37    return U
38
39# Example
40np.random.seed(42)
41n_samples, n_assets = 1000, 3
42
43# Generate correlated returns
44mean = [0.05, 0.07, 0.06]
45cov = [[0.04, 0.02, 0.01],
46       [0.02, 0.06, 0.015],
47       [0.01, 0.015, 0.05]]
48returns = np.random.multivariate_normal(mean, cov, n_samples)
49
50# Fit and transform
51marginals = fit_marginals(returns)
52U = transform_to_uniform(returns, marginals)
53
54print("Uniform data shape:", U.shape)
55print("Uniform data range:", U.min(), "to", U.max())
56

3. Gaussian Copula #

3.1 Definition #

The Gaussian copula with correlation matrix $\Sigma$ : $C^{Gauss}(u_1, \ldots, u_n; \Sigma) = \Phi_\Sigma(\Phi^{-1}(u_1), \ldots, \Phi^{-1}(u_n))$

where:

$\Phi$ = standard normal CDF
$\Phi_\Sigma$ = multivariate normal CDF with correlation $\Sigma$

Properties:

No tail dependence (independence in extremes)
Easy to fit and simulate
Used in CDO pricing (infamous in 2008 crisis)

3.2 Implementation #

python

1class GaussianCopula:
2    """Gaussian copula for multivariate dependence."""
3    
4    def __init__(self, correlation_matrix):
5        self.corr = np.array(correlation_matrix)
6        self.n_dim = len(correlation_matrix)
7    
8    def fit(self, U):
9        """
10        Fit correlation matrix from uniform data.
11        
12        Args:
13            U: Uniform data of shape (n_samples, n_dim)
14        """
15        # Transform to normal
16        Z = stats.norm.ppf(U)
17        
18        # Compute correlation
19        self.corr = np.corrcoef(Z.T)
20        
21        return self
22    
23    def sample(self, n_samples):
24        """
25        Generate samples from Gaussian copula.
26        
27        Returns:
28            Uniform samples of shape (n_samples, n_dim)
29        """
30        # Generate correlated normals
31        Z = np.random.multivariate_normal(
32            np.zeros(self.n_dim),
33            self.corr,
34            n_samples
35        )
36        
37        # Transform to uniform
38        U = stats.norm.cdf(Z)
39        
40        return U
41    
42    def pdf(self, U):
43        """Copula density."""
44        Z = stats.norm.ppf(U)
45        
46        # Multivariate normal density / product of marginal densities
47        det_corr = np.linalg.det(self.corr)
48        inv_corr = np.linalg.inv(self.corr)
49        
50        exponent = -0.5 * np.sum(Z @ (inv_corr - np.eye(self.n_dim)) * Z, axis=1)
51        
52        return np.exp(exponent) / np.sqrt(det_corr)
53
54# Example usage
55copula = GaussianCopula(correlation_matrix=[[1, 0.5], [0.5, 1]])
56U_samples = copula.sample(n_samples=1000)
57
58# Visualize
59plt.figure(figsize=(10, 5))
60plt.subplot(1, 2, 1)
61plt.scatter(U_samples[:, 0], U_samples[:, 1], alpha=0.5, s=10)
62plt.xlabel('U1')
63plt.ylabel('U2')
64plt.title('Gaussian Copula Samples (rho=0.5)')
65plt.grid(True, alpha=0.3)
66
67plt.subplot(1, 2, 2)
68plt.hist2d(U_samples[:, 0], U_samples[:, 1], bins=30, cmap='Blues')
69plt.xlabel('U1')
70plt.ylabel('U2')
71plt.title('Density')
72plt.colorbar()
73plt.tight_layout()
74plt.show()
75

4. Student-t Copula #

4.1 Motivation #

Key advantage: Exhibits tail dependence – assets crash together.

Student-t copula with correlation $\Sigma$ and degrees of freedom $\nu$ : $C^t(u_1, \ldots, u_n; \Sigma, \nu) = t_{\Sigma, \nu}(t_\nu^{-1}(u_1), \ldots, t_\nu^{-1}(u_n))$

Tail dependence coefficient: $\lambda = 2 \bar{t}_{\nu+1}\left(-\sqrt{\frac{(\nu+1)(1-\rho)}{1+\rho}}\right)$

where $\rho$ is correlation and $\bar{t}$ is survival function.

4.2 Implementation #

python

1class StudentTCopula:
2    """Student-t copula with tail dependence."""
3    
4    def __init__(self, correlation_matrix, df):
5        self.corr = np.array(correlation_matrix)
6        self.df = df  # Degrees of freedom
7        self.n_dim = len(correlation_matrix)
8    
9    def fit(self, U, method='mle'):
10        """
11        Fit correlation and degrees of freedom.
12        
13        Uses maximum likelihood estimation.
14        """
15        from scipy.optimize import minimize
16        
17        def neg_log_likelihood(params):
18            df = params[0]
19            # Transform to t-distributed
20            X = stats.t.ppf(U, df)
21            
22            # Compute correlation
23            corr = np.corrcoef(X.T)
24            
25            # Log-likelihood (simplified)
26            try:
27                mvt = stats.multivariate_t(loc=np.zeros(self.n_dim),
28                                          shape=corr, df=df)
29                return -np.sum(mvt.logpdf(X))
30            except:
31                return 1e10
32        
33        # Optimize
34        result = minimize(neg_log_likelihood, x0=[10], bounds=[(2.1, 30)])
35        self.df = result.x[0]
36        
37        # Recompute correlation with optimal df
38        X = stats.t.ppf(U, self.df)
39        self.corr = np.corrcoef(X.T)
40        
41        return self
42    
43    def sample(self, n_samples):
44        """Generate samples from t-copula."""
45        # Generate from multivariate t
46        from scipy.stats import multivariate_t
47        
48        mvt = multivariate_t(loc=np.zeros(self.n_dim),
49                            shape=self.corr, df=self.df)
50        X = mvt.rvs(n_samples)
51        
52        # Transform to uniform
53        U = stats.t.cdf(X, self.df)
54        
55        return U
56    
57    def tail_dependence(self, rho):
58        """
59        Compute tail dependence coefficient.
60        
61        Args:
62            rho: Correlation between two assets
63        """
64        from scipy.stats import t as t_dist
65        
66        arg = -np.sqrt((self.df + 1) * (1 - rho) / (1 + rho))
67        lambda_tail = 2 * t_dist.sf(arg, self.df + 1)
68        
69        return lambda_tail
70
71# Example
72t_copula = StudentTCopula(correlation_matrix=[[1, 0.5], [0.5, 1]], df=5)
73U_t = t_copula.sample(n_samples=1000)
74
75# Compare tail behavior
76plt.figure(figsize=(14, 5))
77
78plt.subplot(1, 3, 1)
79plt.scatter(U_samples[:, 0], U_samples[:, 1], alpha=0.3, s=10, label='Gaussian')
80plt.xlabel('U1')
81plt.ylabel('U2')
82plt.title('Gaussian Copula')
83plt.grid(True, alpha=0.3)
84
85plt.subplot(1, 3, 2)
86plt.scatter(U_t[:, 0], U_t[:, 1], alpha=0.3, s=10, label='Student-t', color='red')
87plt.xlabel('U1')
88plt.ylabel('U2')
89plt.title('Student-t Copula (df=5)')
90plt.grid(True, alpha=0.3)
91
92plt.subplot(1, 3, 3)
93# Focus on lower tail
94mask_gauss = (U_samples[:, 0] < 0.1) & (U_samples[:, 1] < 0.1)
95mask_t = (U_t[:, 0] < 0.1) & (U_t[:, 1] < 0.1)
96plt.scatter(U_samples[mask_gauss, 0], U_samples[mask_gauss, 1], 
97           alpha=0.5, s=20, label='Gaussian')
98plt.scatter(U_t[mask_t, 0], U_t[mask_t, 1], 
99           alpha=0.5, s=20, label='Student-t', color='red')
100plt.xlabel('U1')
101plt.ylabel('U2')
102plt.title('Lower Tail (U1, U2 < 0.1)')
103plt.legend()
104plt.grid(True, alpha=0.3)
105
106plt.tight_layout()
107plt.show()
108
109print("Tail dependence (t-copula):", t_copula.tail_dependence(0.5))
110

5. Archimedean Copulas #

5.1 Families #

Archimedean copulas are defined via generator function $\phi$ : $C(u_1, u_2) = \phi^{-1}(\phi(u_1) + \phi(u_2))$

Clayton Copula (lower tail dependence): $C(u_1, u_2) = (u_1^{-\theta} + u_2^{-\theta} - 1)^{-1/\theta}$

Gumbel Copula (upper tail dependence): $C(u_1, u_2) = \exp\left(-\left[(-\log u_1)^\theta + (-\log u_2)^\theta\right]^{1/\theta}\right)$

Frank Copula (symmetric, no tail dependence): $C(u_1, u_2) = -\frac{1}{\theta} \log\left(1 + \frac{(e^{-\theta u_1} - 1)(e^{-\theta u_2} - 1)}{e^{-\theta} - 1}\right)$

5.2 Implementation #

python

1class ClaytonCopula:
2    """Clayton copula with lower tail dependence."""
3    
4    def __init__(self, theta):
5        self.theta = theta  # theta > 0
6    
7    def cdf(self, u1, u2):
8        """Copula CDF."""
9        return (u1**(-self.theta) + u2**(-self.theta) - 1)**(-1/self.theta)
10    
11    def sample(self, n_samples):
12        """Generate samples using conditional sampling."""
13        U1 = np.random.uniform(0, 1, n_samples)
14        V = np.random.uniform(0, 1, n_samples)
15        
16        # Conditional distribution
17        U2 = (U1**(-self.theta) * (V**(-self.theta/(1+self.theta)) - 1) + 1)**(-1/self.theta)
18        
19        return np.column_stack([U1, U2])
20    
21    def tail_dependence_lower(self):
22        """Lower tail dependence coefficient."""
23        return 2**(-1/self.theta)
24
25class GumbelCopula:
26    """Gumbel copula with upper tail dependence."""
27    
28    def __init__(self, theta):
29        self.theta = theta  # theta >= 1
30    
31    def cdf(self, u1, u2):
32        """Copula CDF."""
33        return np.exp(-((-np.log(u1))**self.theta + 
34                        (-np.log(u2))**self.theta)**(1/self.theta))
35    
36    def sample(self, n_samples):
37        """Generate samples."""
38        # Use Laplace transform method
39        from scipy.stats import levy_stable
40        
41        # Generate from stable distribution
42        alpha = 1 / self.theta
43        stable = levy_stable(alpha=alpha, beta=1)
44        V = stable.rvs(n_samples)
45        
46        E1 = np.random.exponential(1, n_samples)
47        E2 = np.random.exponential(1, n_samples)
48        
49        U1 = np.exp(-(E1 / V)**(1/self.theta))
50        U2 = np.exp(-(E2 / V)**(1/self.theta))
51        
52        return np.column_stack([U1, U2])
53    
54    def tail_dependence_upper(self):
55        """Upper tail dependence coefficient."""
56        return 2 - 2**(1/self.theta)
57
58# Compare copulas
59clayton = ClaytonCopula(theta=2)
60gumbel = GumbelCopula(theta=2)
61
62U_clayton = clayton.sample(1000)
63U_gumbel = gumbel.sample(1000)
64
65plt.figure(figsize=(14, 5))
66
67plt.subplot(1, 3, 1)
68plt.scatter(U_clayton[:, 0], U_clayton[:, 1], alpha=0.3, s=10)
69plt.title('Clayton (Lower Tail Dep)')
70plt.xlabel('U1')
71plt.ylabel('U2')
72plt.grid(True, alpha=0.3)
73
74plt.subplot(1, 3, 2)
75plt.scatter(U_gumbel[:, 0], U_gumbel[:, 1], alpha=0.3, s=10, color='red')
76plt.title('Gumbel (Upper Tail Dep)')
77plt.xlabel('U1')
78plt.ylabel('U2')
79plt.grid(True, alpha=0.3)
80
81plt.subplot(1, 3, 3)
82print("Clayton lower tail dep:", clayton.tail_dependence_lower())
83print("Gumbel upper tail dep:", gumbel.tail_dependence_upper())
84plt.text(0.1, 0.5, 'Clayton lower tail:\n{:.3f}'.format(clayton.tail_dependence_lower()),
85         fontsize=12)
86plt.text(0.1, 0.3, 'Gumbel upper tail:\n{:.3f}'.format(gumbel.tail_dependence_upper()),
87         fontsize=12, color='red')
88plt.xlim(0, 1)
89plt.ylim(0, 1)
90plt.title('Tail Dependence Comparison')
91plt.axis('off')
92
93plt.tight_layout()
94plt.show()
95

6. Portfolio Risk Simulation #

6.1 Complete Workflow #

python

1def simulate_portfolio_returns(
2    marginals,
3    copula,
4    weights,
5    n_scenarios=10000
6):
7    """
8    Simulate portfolio returns using copula.
9    
10    Args:
11        marginals: List of (dist_name, params) for each asset
12        copula: Fitted copula object
13        weights: Portfolio weights
14        n_scenarios: Number of scenarios
15        
16    Returns:
17        Portfolio returns
18    """
19    # Sample from copula
20    U = copula.sample(n_scenarios)
21    
22    # Transform to asset returns
23    n_assets = len(marginals)
24    returns = np.zeros((n_scenarios, n_assets))
25    
26    for i in range(n_assets):
27        dist_name, params = marginals[i]
28        if dist_name == 't':
29            returns[:, i] = stats.t.ppf(U[:, i], *params)
30        elif dist_name == 'norm':
31            returns[:, i] = stats.norm.ppf(U[:, i], *params)
32    
33    # Compute portfolio returns
34    portfolio_returns = returns @ weights
35    
36    return portfolio_returns, returns
37
38# Example: 3-asset portfolio
39weights = np.array([0.4, 0.3, 0.3])
40
41# Fit copula to historical data
42t_copula_3d = StudentTCopula(
43    correlation_matrix=[[1, 0.6, 0.4],
44                       [0.6, 1, 0.5],
45                       [0.4, 0.5, 1]],
46    df=6
47)
48
49# Simulate
50port_returns, asset_returns = simulate_portfolio_returns(
51    marginals,
52    t_copula_3d,
53    weights,
54    n_scenarios=10000
55)
56
57# Risk metrics
58VaR_95 = np.percentile(port_returns, 5)
59CVaR_95 = np.mean(port_returns[port_returns <= VaR_95])
60
61print("Portfolio VaR (95%): {:.2%}".format(VaR_95))
62print("Portfolio CVaR (95%): {:.2%}".format(CVaR_95))
63
64# Plot distribution
65plt.figure(figsize=(12, 5))
66
67plt.subplot(1, 2, 1)
68plt.hist(port_returns, bins=50, alpha=0.7, density=True)
69plt.axvline(VaR_95, color='r', linestyle='--', label='VaR 95%')
70plt.axvline(CVaR_95, color='darkred', linestyle='--', label='CVaR 95%')
71plt.xlabel('Portfolio Return')
72plt.ylabel('Density')
73plt.title('Portfolio Return Distribution')
74plt.legend()
75plt.grid(True, alpha=0.3)
76
77plt.subplot(1, 2, 2)
78plt.scatter(asset_returns[:, 0], asset_returns[:, 1], alpha=0.1, s=1)
79plt.xlabel('Asset 1 Return')
80plt.ylabel('Asset 2 Return')
81plt.title('Joint Return Distribution')
82plt.grid(True, alpha=0.3)
83
84plt.tight_layout()
85plt.show()
86

6.2 Stress Testing #

python

1def stress_test_portfolio(marginals, copula, weights, stress_scenarios):
2    """
3    Stress test portfolio under extreme scenarios.
4    
5    Args:
6        stress_scenarios: List of (asset_idx, percentile) tuples
7    """
8    results = []
9    
10    for scenario_name, conditions in stress_scenarios.items():
11        # Generate scenarios
12        U = copula.sample(10000)
13        returns = np.zeros((10000, len(marginals)))
14        
15        for i in range(len(marginals)):
16            dist_name, params = marginals[i]
17            returns[:, i] = stats.t.ppf(U[:, i], *params)
18        
19        # Apply stress conditions
20        mask = np.ones(10000, dtype=bool)
21        for asset_idx, percentile in conditions:
22            threshold = np.percentile(returns[:, asset_idx], percentile)
23            mask &= (returns[:, asset_idx] <= threshold)
24        
25        # Compute stressed portfolio returns
26        stressed_returns = returns[mask] @ weights
27        
28        results.append({
29            'scenario': scenario_name,
30            'mean': np.mean(stressed_returns),
31            'VaR_95': np.percentile(stressed_returns, 5),
32            'worst': np.min(stressed_returns)
33        })
34    
35    return results
36
37# Define stress scenarios
38stress_scenarios = {
39    'Market Crash': [(0, 5), (1, 5), (2, 5)],  # All assets in bottom 5%
40    'Asset 1 Crash': [(0, 1)],  # Asset 1 in bottom 1%
41    'Correlation Breakdown': [(0, 5), (1, 95)]  # Asset 1 down, Asset 2 up
42}
43
44stress_results = stress_test_portfolio(marginals, t_copula_3d, weights, stress_scenarios)
45
46for result in stress_results:
47    print("Scenario: {}".format(result['scenario']))
48    print("  Mean return: {:.2%}".format(result['mean']))
49    print("  VaR 95%: {:.2%}".format(result['VaR_95']))
50    print("  Worst case: {:.2%}".format(result['worst']))
51    print()
52

7. Production Checklist #

Validate copula choice – Test goodness-of-fit
Check tail dependence – Measure empirically
Use robust estimation – Handle outliers
Implement model selection – AIC/BIC comparison
Cache fitted copulas – Expensive to refit
Monitor correlation stability – Rolling window
Handle high dimensions – Vine copulas for 10+ assets
Validate simulations – Compare to historical
Document assumptions – Stationarity, ergodicity
Implement backtesting – VaR exceedances

References #

Nelsen, R. (2006). An Introduction to Copulas. Springer.
McNeil, A., Frey, R., & Embrechts, P. (2005). Quantitative Risk Management. Princeton.
Joe, H. (2014). Dependence Modeling with Copulas. CRC Press.
Cherubini, U., Luciano, E., & Vecchiato, W. (2004). Copula Methods in Finance. Wiley.

NordVarg Team

Join 1,000+ Engineers

Related Posts

NordVarg Team

Join 1,000+ Engineers

Related Posts