¿Cómo tratar las cópulas?

Question

Quiero modelar un par de variables $(X,Y)$ utilizando cópulas. Mi idea es modelar las distribuciones marginales y luego utilizar cópulas para combinar las distribuciones marginales y la cópula para obtener una distribución conjunta.

Sin embargo, no estoy seguro de cómo funciona.

Digamos que $X$ y $Y$ siguen distribuciones Gamma y que he encontrado los parámetros óptimos para ajustar mis curvas. El siguiente es mi código (Python) para ajustar la distribución Gamma a $X$ , utilizo casi lo mismo para $Y$ .

# x = [80, 100, 285, 186, 134, 89, 78, 45, 98, ...]
sns.distplot(x)
fit_alpha, fit_loc, fit_beta = stats.gamma.fit(x)
rv = stats.gamma(a=fit_alpha, loc=fit_loc, scale=fit_beta)
absc = np.linspace(0, 500, 100)
sns.lineplot(absc, rv.pdf(absc))

Ahora puedo trazar las funciones de distribución acumulativa de mis observaciones (las fdc son las fdc de las distribuciones Gamma ajustadas), y compruebo que las gráficas obtenidas son (casi) las de las distribuciones uniformes (lo que creo que es una condición para usar cópulas). Igual que antes, aquí está el código para $X$ pero hice lo mismo para $Y$ .

x_prime = stats.gamma.cdf(x, a=fit_alpha, loc=fit_loc, scale=fit_beta)
plt.hist(x_prime);

Entonces puedo mostrar el gráfico de dispersión de x_prime y y_prime (que son las cdf ajustadas de mis observaciones). Creo que lo que obtengo en esta fase es un cópula .

plt.scatter(x_prime, y_prime, alpha=0.22)

¿Fue un procedimiento adecuado? Y lo que es más importante, ¿qué debo hacer con eso para obtener una distribución conjunta?

Answer 1

Al hacerlo, la modelización de la función de cópula dependerá del modelo paramétrico que haya utilizado para sus marginales. Incluso si conoce el verdadero modelo que generó sus datos (lo que rara vez ocurre en los problemas del mundo real), el método que utilizó para estimar los parámetros no es exacto y la cópula que obtenga seguirá dependiendo de los marginales.

Para evitarlo, hay que transformar los datos en datos del rango . Por ejemplo, supongamos que tenemos un vector aleatorio bidimensional $\mathbf{X} = (X_1, X_2)$ y una muestra de 4 observaciones $\mathcal{S} = \left\{(-1,3), (0.5,1), (0, -2), (1,0)\right\}$ . La variable de rango $\mathbf{U}_i[j]$ viene dado por el rango de $\mathbf{X}_i[j]$ para todos $1\leq j \leq 4$ . Así que de mi ejemplo tienes $\mathcal{R} = \left\{(1,4), (3,3), (2, 1), (4,2)\right\}$ .

Utilizando esta transformación, se pueden renormalizar los datos para tener valores en $[0, 1]$ y obtener la cópula empírica. A continuación, eligiendo un modelo paramétrico para la cópula, se pueden estimar los mejores parámetros utilizando, por ejemplo, el método MLE.

Termino con un ejemplo utilizando la biblioteca OpenTURNS ( https://openturns.github.io/www/index.html ) que ya tiene un montón de métodos implementados para tratar con cópulas. Digamos que tengo datos generados a partir de una distribución que tiene una cópula gaussiana con un parámetro de correlación $\rho = 0.8$ y marginales Gamma con parámetro de tasa $k = 1$ y el parámetro de forma $\lambda = 1$ :

import openturns as ot
import matplotlib.pyplot as plt
from openturns.viewer import View

ot.RandomGenerator.SetSeed(42)

marginals = [ot.Gamma()] * 2 # default parameters are [1, 1]
copula = ot.NormalCopula(ot.CorrelationMatrix([[1, 0.8], [0.8, 1]]))

distribution = ot.ComposedDistribution(marginals, copula)

sample = distribution.getSample(1000)}

Si transformo mis datos en datos de rango, obtengo la cópula (la normalización elegida es para que los datos estén en $(0,1)$ ):

rank_sample = (sample.rank() + 0.5) / (sample.getSize() + 1)

Entonces, puedo estimar la función de cópula utilizando el modelo que quiera. Como conozco el verdadero modelo de cópula, lo utilizaré:

copula_fit = ot.NormalCopulaFactory().build(sample)
print(copula_fit.getParameter())
# returns 0.796686

Y entonces, puedo estimar los marginales:

marginals_fit = [ot.ExponentialFactory().build(sample.getMarginal(0)),
                 ot.ExponentialFactory().build(sample.getMarginal(1))]
print([mf.getParameter() for mf in marginals_fit])
# returns [1.11684,1.08329] and [1.10042,1.07048]

Nótese que he elegido estimar la cópula y luego los marginales pero podría haber hecho lo contrario ya que los dos problemas son independientes.

Por último, puedes fusionar los marginales y la cópula aprendida en una distribución y trazar su densidad:

fig, ax = plt.subplots()
distribution_fit = ot.ComposedDistribution(marginals_fit, copula_fit)
View(distribution_fit.drawPDF(), figure=fig, axes=[ax], add_legend=False)
plt.show()

O puede evaluar la densidad en cualquier punto:

print(distribution_fit.computePDF([0.5, 0.5])) # Returns 0.6708