Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

Simulación de variables aleatorias dependientes utilizando cópulas

Este ejemplo muestra cómo utilizar cópulas para generar datos de distribuciones multivariadas cuando hay relaciones complicadas entre las variables, o cuando las variables individuales provienen de distribuciones diferentes.

MATLAB® es una herramienta ideal para ejecutar simulaciones que incorporan entradas aleatorias o ruido. Estadísticas y machine learning Toolbox™ proporciona funciones para crear secuencias de datos aleatorios de acuerdo con muchas distribuciones univariadas comunes. El cuadro de herramientas también incluye algunas funciones para generar datos aleatorios de distribuciones multivariadas, como la t multivariada normal y multivariada. Sin embargo, no existe una forma incorporada de generar distribuciones multivariadas para todas las distribuciones marginales, o en los casos en que las variables individuales proceden de distribuciones diferentes.

Recientemente, los cópulas se han popularizado en los modelos de simulación. Las copulas son funciones que describen las dependencias entre las variables y proporcionan una forma de crear distribuciones para modelar datos multivariados correlacionados. Mediante una cópula, un analista de datos puede construir una distribución multivariada especificando distribuciones univariadas marginales y eligiendo una cópula concreta para proporcionar una estructura de correlación entre variables. Las distribuciones bivariadas, así como las distribuciones en dimensiones más altas, son posibles. En este ejemplo, discutimos cómo utilizar cópulas para generar datos aleatorios multivariados dependientes en MATLAB, utilizando Statistics y machine learning Toolbox.

Dependencia entre entradas de simulación

Una de las decisiones de diseño para una simulación Monte-Carlo es una elección de distribuciones de probabilidad para las entradas aleatorias. La selección de una distribución para cada variable individual es a menudo sencilla, pero decidir qué dependencias deben existir entre las entradas puede no serlo. Idealmente, los datos de entrada a una simulación deben reflejar lo que se sabe acerca de la dependencia entre las cantidades reales que se modelan. Sin embargo, puede haber poca o ninguna información sobre la cual basar cualquier dependencia en la simulación, y en tales casos, es una buena idea experimentar con diferentes posibilidades, con el fin de determinar la sensibilidad del modelo.

Sin embargo, puede ser difícil generar realmente entradas aleatorias con dependencia cuando tienen distribuciones que no son de una distribución multivariada estándar. Además, algunas de las distribuciones multivariadas estándar pueden modelar solo tipos muy limitados de dependencia. Siempre es posible hacer que las entradas sean independientes, y si bien eso es una elección simple, no siempre es sensato y puede conducir a conclusiones equivocadas.

Por ejemplo, una simulación Monte-Carlo de riesgo financiero podría tener insumos aleatorios que representen diferentes fuentes de pérdidas de seguro. Estas entradas pueden ser modeladas como variables aleatorias lognormal. Una pregunta razonable es cómo la dependencia entre estas dos entradas afecta a los resultados de la simulación. De hecho, podría ser sabido a partir de datos reales que las mismas condiciones aleatorias afectan a ambas fuentes, e ignorando que en la simulación podría conducir a conclusiones equivocadas.

La simulación de variables aleatorias lognormal independientes es trivial. La forma más sencilla sería utilizar la función.lognrnd Aquí, usaremos la función para generar n pares de variables aleatorias normales independientes, y luego exponerlos.mvnrnd Observe que la matriz de covarianza utilizada aquí es diagonal, es decir, independencia entre las columnas de Z.

n = 1000; sigma = .5; SigmaInd = sigma.^2 .* [1 0; 0 1] 
 SigmaInd =      0.2500         0          0    0.2500  
ZInd = mvnrnd([0 0], SigmaInd, n); XInd = exp(ZInd); plot(XInd(:,1),XInd(:,2),'.'); axis equal; axis([0 5 0 5]); xlabel('X1'); ylabel('X2'); 

Los r.v. de lognormal dependientes bivariables también son fáciles de generar, utilizando una matriz de covarianza con términos fuera de la diagonal no-cero.

rho = .7; SigmaDep = sigma.^2 .* [1 rho; rho 1] 
 SigmaDep =      0.2500    0.1750     0.1750    0.2500  
ZDep = mvnrnd([0 0], SigmaDep, n); XDep = exp(ZDep); 

Un segundo gráfico de dispersión ilustra la diferencia entre estas dos distribuciones bivariadas.

plot(XDep(:,1),XDep(:,2),'.'); axis equal; axis([0 5 0 5]); xlabel('X1'); ylabel('X2'); 

Está claro que hay más tendencia en el segundo DataSet para que los valores grandes de x1 se asocien con valores grandes de x2, y de manera similar para valores pequeños. Esta dependencia está determinada por el parámetro de correlación, Rho, de la normal bivariada subyacente. Las conclusiones extraídas de la simulación bien podrían depender de si x1 y x2 fueron generados con dependencia o no.

La distribución lognormal bivariada es una solución simple en este caso, y por supuesto se generaliza fácilmente a dimensiones más altas y casos donde las distribuciones marginales son lognormals.different También existen otras distribuciones multivariadas, por ejemplo, la t multivariada y las distribuciones Dirichlet se utilizan para simular variables aleatorias t y beta dependientes, respectivamente. Pero la lista de distribuciones multivariadas simples no es larga, y sólo se aplican en casos donde los marginales están todos en la misma familia (o incluso las mismas distribuciones exactas). Esto puede ser una limitación real en muchas situaciones.

Un método más general para construir distribuciones bivariadas dependientes

Aunque la construcción anterior que crea un lognormal bivariado es simple, sirve para ilustrar un método que es más generalmente aplicable. En primer lugar, generamos pares de valores a partir de una distribución normal bivariada. Existe una dependencia estadística entre estas dos variables, y cada una tiene una distribución marginal normal. A continuación, una transformación (la función exponencial) se aplica por separado a cada variable, cambiando las distribuciones marginales en lognormals. Las variables transformadas todavía tienen una dependencia estadística.

Si se pudiera encontrar una transformación adecuada, este método podría generalizarse para crear vectores aleatorios bivariados dependientes con otras distribuciones marginales. De hecho, existe un método general de construcción de tal transformación, aunque no tan simple como la mera exponenciación.

Por definición, aplicar el CDF normal (denotado aquí por PHI) a una variable aleatoria normal estándar da como resultado un r.v. que es uniforme en el intervalo [0,1]. Para ver esto, si Z tiene una distribución normal estándar, entonces el CDF de U = PHI (Z) es

  Pr{U <= u0} = Pr{PHI(Z) <= u0} = Pr{Z <= PHI^(-1)(u0)} = u0,

y esa es la CDF de una U (0,1) r.v. histogramas de algunos valores normales y transformados simulados demuestran ese hecho.

n = 1000; z = normrnd(0,1,n,1); hist(z,-3.75:.5:3.75); xlim([-4 4]); title('1000 Simulated N(0,1) Random Values'); xlabel('Z'); ylabel('Frequency'); 

u = normcdf(z); hist(u,.05:.1:.95); title('1000 Simulated N(0,1) Values Transformed to U(0,1)'); xlabel('U'); ylabel('Frequency'); 

Ahora, tomando prestado de la teoría de la generación de números aleatorios univariados, aplicando el CDF inverso de cualquier distribución F a una variable aleatoria U (0,1) da como resultado un r.v. cuya distribución es exactamente F. Esto se conoce como el método de inversión. La prueba es esencialmente lo opuesto a la prueba anterior para el caso delantero. Otro histograma ilustra la transformación a una distribución gamma.

x = gaminv(u,2,1); hist(x,.25:.5:9.75); title('1000 Simulated N(0,1) Values Transformed to Gamma(2,1)'); xlabel('X'); ylabel('Frequency'); 

Esta transformación de dos pasos se puede aplicar a cada variable de una normal bivariada estándar, creando r.v. dependientes con distribuciones marginales arbitrarias. Dado que la transformación funciona en cada componente por separado, los dos r.v. resultantes no tienen ni siquiera las mismas distribuciones marginales. La transformación se define como

  Z = [Z1 Z2] ~ N([0 0],[1 rho; rho 1])   U = [PHI(Z1) PHI(Z2)]   X = [G1(U1) G2(U2)]

donde G1 y G2 son CDFs inversas de dos distribuciones posiblemente diferentes. Por ejemplo, podemos generar vectores aleatorios a partir de una distribución bivariada con los marginales t (5) y gamma (2, 1).

n = 1000; rho = .7; Z = mvnrnd([0 0], [1 rho; rho 1], n); U = normcdf(Z); X = [gaminv(U(:,1),2,1) tinv(U(:,2),5)]; 

Esta trama tiene histogramas junto a un gráfico de dispersión para mostrar tanto las distribuciones marginales, como la dependencia.

[n1,ctr1] = hist(X(:,1),20); [n2,ctr2] = hist(X(:,2),20); subplot(2,2,2); plot(X(:,1),X(:,2),'.'); axis([0 12 -8 8]); h1 = gca; title('1000 Simulated Dependent t and Gamma Values'); xlabel('X1 ~ Gamma(2,1)'); ylabel('X2 ~ t(5)'); subplot(2,2,4); bar(ctr1,-n1,1); axis([0 12 -max(n1)*1.1 0]); axis('off'); h2 = gca; subplot(2,2,1); barh(ctr2,-n2,1); axis([-max(n2)*1.1 0 -8 8]); axis('off'); h3 = gca; h1.Position = [0.35 0.35 0.55 0.55]; h2.Position = [.35 .1 .55 .15]; h3.Position = [.1 .35 .15 .55]; colormap([.8 .8 1]); 

Rango coeficientes de correlación

La dependencia entre x1 y x2 en esta construcción está determinada por el parámetro de correlación, Rho, de la normal bivariada subyacente. Sin embargo, es cierto que la correlación lineal de x1 y x2 es Rho.not Por ejemplo, en el caso lognormal original, hay un formulario cerrado para esa correlación:

  cor(X1,X2) = (exp(rho.*sigma.^2) - 1) ./ (exp(sigma.^2) - 1)

que es estrictamente menos que Rho a menos que Rho es exactamente uno. En casos más generales, sin embargo, como la construcción gamma/t anterior, la correlación lineal entre x1 y x2 es difícil o imposible de expresar en términos de Rho, pero las simulaciones se pueden utilizar para mostrar que ocurre el mismo efecto.

Esto se debe a que el coeficiente de correlación lineal expresa la dependencia entre r.v. y cuando las transformaciones no lineales se aplican a los r.v., no se preserva la correlación lineal.linear En cambio, un coeficiente de correlación de rango, como el Tau de Kendall o el Rho de Spearman, es más apropiado.

En términos generales, estas correlaciones de rango miden el grado en que los valores grandes o pequeños de una r.v. asocian con valores grandes o pequeños de otro. Sin embargo, a diferencia del coeficiente de correlación lineal, miden la asociación sólo en términos de rangos. Como consecuencia, la correlación de rangos se preserva bajo cualquier transformación monotónica. En particular, el método de transformación que acabamos de describir conserva la correlación de rangos. Por lo tanto, conociendo la correlación de rango de la normal bivariada Z determina exactamente la correlación de rango de la r.v. transformada final de X. Mientras que Rho todavía se necesita para parametrizar la normal bivariada subyacente, Tau de Kendall o Rho de Spearman son más útiles en la descripción de la dependencia entre los r.v., porque son invariantes a la elección de la distribución marginal.

Resulta que para la normalidad bivariada, hay un simple 1-1 mapeo entre Tau de Kendall o Rho de Spearman, y el coeficiente de correlación lineal Rho:

  tau   = (2/pi)*arcsin(rho)     or   rho = sin(tau*pi/2)   rho_s = (6/pi)*arcsin(rho/2)   or   rho = 2*sin(rho_s*pi/6)
subplot(1,1,1); rho = -1:.01:1; tau = 2.*asin(rho)./pi; rho_s = 6.*asin(rho./2)./pi; plot(rho,tau,'-', rho,rho_s,'-', [-1 1],[-1 1],'k:'); axis([-1 1 -1 1]); xlabel('rho'); ylabel('Rank correlation coefficient'); legend('Kendall''s tau', 'Spearman''s rho_s', 'location','northwest'); 

Por lo tanto, es fácil crear la correlación de rango deseada entre x1 y x2, independientemente de sus distribuciones marginales, eligiendo el valor de parámetro Rho correcto para la correlación lineal entre Z1 y Z2.

Observe que para la distribución normal multivariada, la correlación de rango de Spearman es casi idéntica a la correlación lineal. Sin embargo, esto no es cierto una vez que transformamos a las variables aleatorias finales.

Cópulas

El primer paso de la construcción descrita anteriormente define lo que se conoce como una cópula, concretamente, una cópula Gaussiana. Una cópula bivariada es simplemente una distribución de probabilidad en dos variables aleatorias, cada una de cuyas distribuciones marginales es uniforme. Estas dos variables pueden ser completamente independientes, relacionadas de forma determinista (por ejemplo, U2 = U1), o cualquier cosa en el medio. La familia de cópulas Gaussianas bivariadas es parametrizada por Rho = [1 Rho; Rho 1], la matriz de correlación lineal. U1 y U2 se acercan a la dependencia lineal como los enfoques Rho +/-1, y se acercan la independencia completa como Rho enfoques cero.

Los diagramas de dispersión de algunos valores aleatorios simulados para varios niveles de Rho ilustran la gama de diferentes posibilidades para los cópulas Gaussianos:

n = 500; Z = mvnrnd([0 0], [1 .8; .8 1], n); U = normcdf(Z,0,1); subplot(2,2,1); plot(U(:,1),U(:,2),'.'); title('rho = 0.8'); xlabel('U1'); ylabel('U2'); Z = mvnrnd([0 0], [1 .1; .1 1], n); U = normcdf(Z,0,1); subplot(2,2,2); plot(U(:,1),U(:,2),'.'); title('rho = 0.1'); xlabel('U1'); ylabel('U2'); Z = mvnrnd([0 0], [1 -.1; -.1 1], n); U = normcdf(Z,0,1); subplot(2,2,3); plot(U(:,1),U(:,2),'.'); title('rho = -0.1'); xlabel('U1'); ylabel('U2'); Z = mvnrnd([0 0], [1 -.8; -.8 1], n); U = normcdf(Z,0,1); subplot(2,2,4); plot(U(:,1),U(:,2),'.'); title('rho = -0.8'); xlabel('U1'); ylabel('U2'); 

La dependencia entre U1 y U2 está completamente separada de las distribuciones marginales de x1 = G (U1) y x2 = G (U2). X1 y x2 se pueden dar distribuciones marginales, y todavía tienen la misma correlación de rango.any Este es uno de los principales atractivos de los cópulas--permiten esta especificación separada de dependencia y distribución marginal.

t copulas

Se puede construir una familia diferente de cópulas partiendo de una distribución de t bivariada, y transformándola utilizando el CDF t correspondiente. La distribución de la t bivariada se parametriza con Rho, la matriz de correlación lineal, y Nu, los grados de libertad. Así, por ejemplo, podemos hablar de una t (1) o una t (5) cópula, basada en la t multivariada con uno y cinco grados de libertad, respectivamente.

Los diagramas de dispersión de algunos valores aleatorios simulados para varios niveles de Rho ilustran el rango de diferentes posibilidades para t (1) cópulas:

n = 500; nu = 1; T = mvtrnd([1 .8; .8 1], nu, n); U = tcdf(T,nu); subplot(2,2,1); plot(U(:,1),U(:,2),'.'); title('rho = 0.8'); xlabel('U1'); ylabel('U2'); T = mvtrnd([1 .1; .1 1], nu, n); U = tcdf(T,nu); subplot(2,2,2); plot(U(:,1),U(:,2),'.'); title('rho = 0.1'); xlabel('U1'); ylabel('U2'); T = mvtrnd([1 -.1; -.1 1], nu, n); U = tcdf(T,nu); subplot(2,2,3); plot(U(:,1),U(:,2),'.'); title('rho = -0.1'); xlabel('U1'); ylabel('U2'); T = mvtrnd([1 -.8; -.8 1], nu, n); U = tcdf(T,nu); subplot(2,2,4); plot(U(:,1),U(:,2),'.'); title('rho = -0.8'); xlabel('U1'); ylabel('U2'); 

Una cópula t tiene distribuciones marginales uniformes para U1 y U2, al igual que una cópula Gaussiana. La correlación de rango Tau o rho_s entre los componentes en una t cópula es también la misma función de Rho que para un gaussiano. Sin embargo, como demuestran estas parcelas, una cópula de t (1) difiere bastante de una cópula gaussiana, incluso cuando sus componentes tienen la misma correlación de rangos. La diferencia está en su estructura de dependencia. No es sorprendente, ya que el parámetro de grados de libertad nu se hace más grande, una t (nu) cópula se aproxima a la cópula gaussiana correspondiente.

Al igual que con una cópula gaussiana, cualquier distribución marginal se puede imponer sobre una copula. Por ejemplo, utilizando una cópula t con 1 grado de libertad, podemos generar de nuevo vectores aleatorios a partir de una distribución bivariada con los marginales GAM (2, 1) y t (5):

subplot(1,1,1); n = 1000; rho = .7; nu = 1; T = mvtrnd([1 rho; rho 1], nu, n); U = tcdf(T,nu); X = [gaminv(U(:,1),2,1) tinv(U(:,2),5)];  [n1,ctr1] = hist(X(:,1),20); [n2,ctr2] = hist(X(:,2),20); subplot(2,2,2); plot(X(:,1),X(:,2),'.'); axis([0 15 -10 10]); h1 = gca; title('1000 Simulated Dependent t and Gamma Values'); xlabel('X1 ~ Gamma(2,1)'); ylabel('X2 ~ t(5)'); subplot(2,2,4); bar(ctr1,-n1,1); axis([0 15 -max(n1)*1.1 0]); axis('off'); h2 = gca; subplot(2,2,1); barh(ctr2,-n2,1); axis([-max(n2)*1.1 0 -10 10]); axis('off'); h3 = gca; h1.Position = [0.35 0.35 0.55 0.55]; h2.Position = [.35 .1 .55 .15]; h3.Position = [.1 .35 .15 .55]; colormap([.8 .8 1]); 

En comparación con la distribución bivariada gamma/t construida anteriormente, que se basaba en una copula gaussiana, la distribución construida aquí, basada en una copula t (1), tiene las mismas distribuciones marginales y la misma correlación de rangos entre variables, pero una muy estructura de dependencia diferente. Esto ilustra el hecho de que las distribuciones multivariadas no se definen de forma exclusiva por sus distribuciones marginales o por sus correlaciones. La elección de una cópula en particular en una aplicación puede basarse en datos reales observados, o diferentes cópulas pueden utilizarse como una forma de determinar la sensibilidad de los resultados de la simulación a la distribución de entrada.

Las copulas de orden superior

Las cópulas gaussiana y t son conocidas como cópulas elípticas. Es fácil generalizar cópulas elípticas a un mayor número de dimensiones. Por ejemplo, podemos simular datos de una distribución trivariada con los marginales gamma (2, 1), beta (2, 2) y t (5) utilizando una cópula gaussiana como se indica a continuación.

subplot(1,1,1); n = 1000; Rho = [1 .4 .2; .4 1 -.8; .2 -.8 1]; Z = mvnrnd([0 0 0], Rho, n); U = normcdf(Z,0,1); X = [gaminv(U(:,1),2,1) betainv(U(:,2),2,2) tinv(U(:,3),5)]; plot3(X(:,1),X(:,2),X(:,3),'.'); grid on; view([-55, 15]); xlabel('U1'); ylabel('U2'); zlabel('U3'); 

Tenga en cuenta que la relación entre el parámetro de correlación lineal Rho y, por ejemplo, Tau de Kendall, contiene para cada entrada en la matriz de correlación rho utilizado aquí. Podemos verificar que las correlaciones de rango de muestreo de los datos son aproximadamente iguales a los valores teóricos.

tauTheoretical = 2.*asin(Rho)./pi 
 tauTheoretical =      1.0000    0.2620    0.1282     0.2620    1.0000   -0.5903     0.1282   -0.5903    1.0000  
tauSample = corr(X, 'type','Kendall') 
 tauSample =      1.0000    0.2655    0.1060     0.2655    1.0000   -0.6076     0.1060   -0.6076    1.0000  

Copulas y distribuciones marginales empíricas

Para simular datos multivariados dependientes usando una cópula, hemos visto que tenemos que especificar

  1) the copula family (and any shape parameters),   2) the rank correlations among variables, and   3) the marginal distributions for each variable

Supongamos que tenemos dos conjuntos de datos de devolución de stock, y nos gustaría ejecutar una simulación de Monte Carlo con entradas que sigan las mismas distribuciones que nuestros datos.

load stockreturns nobs = size(stocks,1); subplot(2,1,1); hist(stocks(:,1),10); xlabel('X1'); ylabel('Frequency'); subplot(2,1,2); hist(stocks(:,2),10); xlabel('X2'); ylabel('Frequency'); 

(Estos dos vectores de datos tienen la misma longitud, pero eso no es crucial.)

Podríamos encajar un modelo paramétrico por separado para cada DataSet, y usar esas estimaciones como nuestras distribuciones marginales. Sin embargo, un modelo paramétrico puede no ser lo suficientemente flexible. En su lugar, podemos usar un modelo empírico para las distribuciones marginales. Sólo necesitamos una forma de calcular el CDF inverso.

El CDF inverso empírico para estos datasets es sólo una función de escalera, con pasos en los valores 1/Nobs, 2/Nobs,... 1. Las alturas de paso son simplemente los datos ordenados.

invCDF1 = sort(stocks(:,1)); n1 = length(stocks(:,1)); invCDF2 = sort(stocks(:,2)); n2 = length(stocks(:,2)); subplot(1,1,1); stairs((1:nobs)/nobs, invCDF1,'b'); hold on; stairs((1:nobs)/nobs, invCDF2,'r'); hold off; legend('X1','X2'); xlabel('Cumulative Probability'); ylabel('X'); 

Para la simulación, es posible que queramos experimentar con diferentes cópulas y correlaciones. Aquí, usaremos una cópula bivariada t (5) con un parámetro de correlación negativa bastante grande.

n = 1000; rho = -.8; nu = 5; T = mvtrnd([1 rho; rho 1], nu, n); U = tcdf(T,nu); X = [invCDF1(ceil(n1*U(:,1))) invCDF2(ceil(n2*U(:,2)))];  [n1,ctr1] = hist(X(:,1),10); [n2,ctr2] = hist(X(:,2),10); subplot(2,2,2); plot(X(:,1),X(:,2),'.'); axis([-3.5 3.5 -3.5 3.5]); h1 = gca; title('1000 Simulated Dependent Values'); xlabel('X1'); ylabel('X2'); subplot(2,2,4); bar(ctr1,-n1,1); axis([-3.5 3.5 -max(n1)*1.1 0]); axis('off'); h2 = gca; subplot(2,2,1); barh(ctr2,-n2,1); axis([-max(n2)*1.1 0 -3.5 3.5]); axis('off'); h3 = gca; h1.Position = [0.35 0.35 0.55 0.55]; h2.Position = [.35 .1 .55 .15]; h3.Position = [.1 .35 .15 .55]; colormap([.8 .8 1]); 

Los histogramas marginales de los datos simulados coinciden estrechamente con los de los datos originales y se convertirían en idénticos a medida que simulamos más pares de valores. Tenga en cuenta que los valores se dibujan a partir de los datos originales y, dado que solo hay 100 observaciones en cada conjunto de datos, los datos simulados son algo "discretos". Una manera de superar esto sería agregar una pequeña cantidad de variación aleatoria, posiblemente distribuida, a los valores simulados finales. Esto equivale a usar una versión suavizada del CDF inverso empírico.