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 mediante copulas

En este ejemplo se muestra cómo utilizar copulas para generar datos a partir de distribuciones multivariadas cuando hay relaciones complicadas entre las variables o cuando las variables individuales proceden de distribuciones diferentes.

MATLAB® es una herramienta ideal para ejecutar simulaciones que incorporan entradas aleatorias o ruido. Statistics and 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 a partir de distribuciones multivariadas, como la t multivariante normal y multivariante. Sin embargo, no hay ninguna forma integrada de generar distribuciones multivariadas para todas las distribuciones marginales, o en los casos en que las variables individuales son de distribuciones diferentes.

Recientemente, las copulas se han vuelto populares en los modelos de simulación. Las copulas son funciones que describen las dependencias entre variables y proporcionan una manera de crear distribuciones para modelar datos multivariantes correlacionados. Mediante una cópula, un analista de datos puede construir una distribución multivariante especificando distribuciones univariadas marginales y eligiendo una cópula determinada 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, se describe cómo utilizar copulas para generar datos aleatorios multivariantes dependientes en MATLAB, mediante Statistics and Machine Learning Toolbox.

Dependencia entre entradas de simulación

Una de las decisiones de diseño para una simulación de Montecarlo es una opción de distribuciones de probabilidad para las entradas aleatorias. Seleccionar una distribución para cada variable individual suele ser sencillo, pero decidir qué dependencias deben existir entre las entradas puede no serlo. Idealmente, los datos de entrada en una simulación deben reflejar lo que se conoce sobre la dependencia entre las cantidades reales que se modelan. Sin embargo, puede haber poca o ninguna información sobre la que 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 solo pueden modelar tipos de dependencia muy limitados. Siempre es posible hacer que las entradas sean independientes, y aunque es una elección sencilla, no siempre es sensata y puede llevar a conclusiones equivocadas.

Por ejemplo, una simulación de Riesgo Financiero de Montecarlo podría tener insumos aleatorios que representen diferentes fuentes de pérdidas de seguros. Estas entradas pueden modelarse como variables aleatorias lognormales. Una pregunta razonable es cómo la dependencia entre estas dos entradas afecta los resultados de la simulación. De hecho, podría saberse a partir de datos reales que las mismas condiciones aleatorias afectan a ambas fuentes, e ignorar que en la simulación podría conducir a conclusiones equivocadas.

La simulación de variables aleatorias lognormales 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 exponentilas.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 registro sortado dependientes también son fáciles de generar, utilizando una matriz de covarianza con términos fuera diagonales distintos de 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); 

Una segunda gráfica 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 conjunto de datos para que los valores grandes de X1 se asocien con valores grandes de X2, e incluso para valores pequeños. Esta dependencia viene 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 se generaron o no 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 en los que las distribuciones marginales son lognormales.diferente También existen otras distribuciones multivariadas, por ejemplo, las distribuciones t multivariada s y 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 los casos en que 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 la construcción de distribuciones bivariadas dependientes

Aunque la construcción anterior que crea un lognormal bivariante es simple, sirve para ilustrar un método que es más aplicable generalmente. En primer lugar, generamos pares de valores a partir de una distribución normal bivariada. Hay 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 lognormales. 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 para construir una transformación de este tipo, aunque no tan simple como la simple exponenciación.

Por definición, la aplicación del CDF normal (denotado aquí por PHI) a una variable aleatoria normal estándar da como resultado una 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 ese es el CDF de un 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, el préstamo de la teoría de la generación de números aleatorios univariados, la aplicación del 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 contrario de la prueba anterior para el caso de reenvío. 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 en 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 ni siquiera necesitan tener 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 CDF inversos de dos distribuciones posiblemente diferentes. Por ejemplo, podemos generar vectores aleatorios a partir de una distribución bivariada con 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 gráfica tiene histogramas junto a una gráfica 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]); 

Coeficientes de correlación de rango

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.No 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 se aplican transformaciones no lineales a esas r.v., la correlación lineal no se conserva.Lineal En su lugar, 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 un r.v. se 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 rango sin reservas bajo cualquier transformación monotónica. En particular, el método de transformación que acaba de describir conserva la correlación de clasificación. Por lo tanto, conocer la correlación de rango de la Z normal bivariante determina exactamente la correlación de rango de la X de r.v. transformada final. Mientras que rho todavía es necesario para parametrizar la normal bivariada subyacente, el tau de Kendall o el rho de Spearman son más útiles para describir la dependencia entre r.v., porque son invariables a la elección de la distribución marginal.

Resulta que para la normal bivariada, hay un simple mapeo 1-1 entre el tau de Kendall o el 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 correcto del parámetro rho para la correlación lineal entre Z1 y Z2.

Observe que para la distribución normal multivariante, 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.

Copulas

El primer paso de la construcción descrita anteriormente define lo que se conoce como una cópula, específicamente, 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 determinísticamente (por ejemplo, U2 a U1), o cualquier cosa en el medio. La familia de copulas gaussianas bivariadas está parametrizada por Rho [1 rho; rho 1], la matriz de correlación lineal. U1 y U2 se acercan a la dependencia lineal a medida que rho se acerca +/- 1, y se acercan a la independencia completa a medida que rho se acerca a cero.

Las gráficas de dispersión de algunos valores aleatorios simulados para varios niveles de rho ilustran el rango de diferentes posibilidades para las copulas gaussas:

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 es completamente independiente de las distribuciones marginales de X1 a G(U1) y X2 a G(U2). X1 y X2 pueden recibir distribuciones marginales y todavía tienen la misma correlación de rango.Cualquier Este es uno de los principales atractivos de las copulas -- permiten esta especificación separada de dependencia y distribución marginal.

t Copulas

Se puede construir una familia diferente de copulas partiendo de una distribución t bivariada y transformándola utilizando el CDF t correspondiente. La distribución 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 copula t(1) o t(5), basada en la t multivariada con uno y cinco grados de libertad, respectivamente.

Las gráficas de dispersión de algunos valores aleatorios simulados para varios niveles de rho ilustran el rango de diferentes posibilidades para t(1) copulas:

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 copula 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 copula t es también la misma función de rho que para un gaussiano. Sin embargo, como demuestran estas gráficas, una copula t(1) difiere bastante de una cópula gaussiana, incluso cuando sus componentes tienen la misma correlación de rango. La diferencia está en su estructura de dependencia. No es de extrañar que, a medida que se hace más grande el parámetro nu de grados de libertad, una copula t(nu) se acerca a la cópula gaussiana correspondiente.

Al igual que con una cópula gaussiana, cualquier distribución marginal se puede imponer a través de una copula t. Por ejemplo, utilizando una copula t con 1 grado de libertad, podemos volver a generar vectores aleatorios a partir de una distribución bivariada con 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 cópula gaussiana, la distribución construida aquí, basada en una copula t(1), tiene las mismas distribuciones marginales y la misma correlación de rango entre variables, pero una estructura de dependencia muy 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 determinada en una aplicación puede basarse en datos reales observados, o diferentes copulas pueden utilizarse como una forma de determinar la sensibilidad de los resultados de la simulación a la distribución de entrada.

Copulas de orden superior

Las copulas gaussianas y t se conocen como copulas elípticas. Es fácil generalizar las copulas elípticas a un mayor número de dimensiones. Por ejemplo, podemos simular datos de una distribución trivariada con marginales Gamma(2,1), Beta(2,2) y t(5) utilizando una cópula gaussiana de la siguiente manera.

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'); 

Observe que la relación entre el parámetro de correlación lineal rho y, por ejemplo, el tau de Kendall, se mantiene para cada entrada en la matriz de correlación que Rho utilizó aquí. Podemos verificar que las correlaciones de rango de muestra 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 multivariadas dependientes utilizando una cópula, hemos visto que necesitamos 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 siguen 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 ajustar un modelo paramétrico por separado a cada conjunto de datos 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 los pasos 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 deseemos experimentar con diferentes copulas y correlaciones. Aquí, usaremos una copula 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 volverían idénticos a medida que simulamos más pares de valores. Observe que los valores se extraen 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 normalmente distribuida, a los valores simulados finales. Esto equivale a utilizar una versión suavizada del CDF inverso empírico.