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.

Modelado de datos de cola con la distribución de Pareto generalizado

En este ejemplo se muestra cómo ajustar los datos de cola a la distribución de Pareto generalizado mediante la estimación de máxima verosimilitud.

Ajustar una distribución paramétrica a los datos a veces resulta en un modelo que concuerda bien con los datos en regiones de alta densidad, pero mal en áreas de baja densidad. Para las distribuciones unimodales, como la normal o Student t, estas regiones de baja densidad se conocen como las "colas" de la distribución. Una razón por la que un modelo puede caber mal en las colas es que por definición, hay menos datos en las colas en las que basar una elección de modelo, y por lo que los modelos a menudo se eligen en función de su capacidad para ajustar los datos cerca del modo. Otra razón podría ser que la distribución de datos reales es a menudo más complicada que los modelos paramétricos habituales.

Sin embargo, en muchas aplicaciones, el ajuste de los datos en la cola es la principal preocupación. La distribución de Pareto generalizado (GP) se desarrolló como una distribución que puede modelar colas de una amplia variedad de distribuciones, basadas en argumentos teóricos. Un enfoque de la adaptación de distribución que involucra al GP es utilizar un ajuste no paramétrico (la función de distribución acumulativa empírica, por ejemplo) en las regiones donde hay muchas observaciones, y para ajustar el GP a la cola (s) de los datos.

La distribución generalizada de Pareto

El Pareto generalizado (GP) es una distribución sesgada a la derecha, parametrizada con un parámetro de forma, k, y un parámetro de escala, Sigma. k también se conoce como el parámetro "índice de cola" y puede ser positivo, cero o negativo.

x = linspace(0,10,1000); plot(x,gppdf(x,-.4,1),'-', x,gppdf(x,0,1),'-', x,gppdf(x,2,1),'-'); xlabel('x / sigma'); ylabel('Probability density'); legend({'k < 0' 'k = 0' 'k > 0'}); 

Tenga en cuenta que para k < 0, el GP tiene cero probabilidad por encima de un límite superior de-(1/k). Para k > = 0, el GP no tiene límite superior. También, el GP se utiliza a menudo en conjunción con un tercer parámetro, umbral que desplaza el límite inferior lejos de cero. No necesitaremos esa generalidad.

La distribución de GP es una generalización tanto de la distribución exponencial (k = 0) como de la distribución de Pareto (k > 0). El GP incluye esas dos distribuciones en una familia más grande para que sea posible un rango continuo de formas.

Simulación de datos de excedencia

La distribución de GP se puede definir constructivamente en términos de excedencias. Comenzando con una distribución de probabilidad cuya cola derecha desciende a cero, como la normal, podemos muestrear valores aleatorios independientemente de esa distribución. Si arreglamos un valor de umbral, desechamos todos los valores que están por debajo del umbral y restamos el umbral de los valores que no se desechan, el resultado se conoce como excedencias. La distribución de las excedencias es aproximadamente un GP. De forma similar, podemos establecer un umbral en la cola izquierda de una distribución e ignorar todos los valores por encima de ese umbral. El umbral debe estar lo suficientemente lejos en la cola de la distribución original para que la aproximación sea razonable.

La distribución original determina el parámetro de forma, k, de la distribución de GP resultante. Las distribuciones cuyas colas se caen como un polinomio, como la t de Student, conducen a un parámetro de forma positiva. Las distribuciones cuyas colas disminuyen exponencialmente, como la normal, corresponden a un parámetro de forma cero. Las distribuciones con colas finitas, como la beta, corresponden a un parámetro de forma negativa.

Las aplicaciones del mundo real para la distribución de GP incluyen modelar extremos de los retornos del mercado de valores y modelar inundaciones extremas. Para este ejemplo, usaremos datos simulados, generados a partir de la distribución t de Student con 5 grados de libertad. Tomaremos el 5% más grande de 2000 observaciones de la distribución t, y luego Restaremos el 95% cuantil para obtener excedencias.

rng(3,'twister'); x = trnd(5,2000,1); q = quantile(x,.95); y = x(x>q) - q; n = numel(y) 
 n =     100  

Ajuste de la distribución utilizando la máxima verosimilitud

La distribución de GP se define para 0 < Sigma e-inf < k < INF. Sin embargo, la interpretación de los resultados de la estimación de máxima verosimilitud es problemática cuando k <-1/2. Afortunadamente, esos casos corresponden a las colas de ajuste de distribuciones como la beta o triangular, y por lo tanto no presentará un problema aquí.

paramEsts = gpfit(y); kHat      = paramEsts(1)   % Tail index parameter sigmaHat  = paramEsts(2)   % Scale parameter 
 kHat =      0.0987   sigmaHat =      0.7156  

Como es de esperar, ya que los datos simulados se generaron utilizando una distribución t, la estimación de k es positiva.

Comprobando el ajuste visualmente

Para evaluar visualmente cuán bueno es el ajuste, trazaremos un histograma escalado de los datos de la cola, superados con la función de densidad del GP que hemos calculado. El histograma se escala de modo que la barra de alturas veces su ancho suma a 1.

bins = 0:.25:7; h = bar(bins,histc(y,bins)/(length(y)*.25),'histc'); h.FaceColor = [.9 .9 .9]; ygrid = linspace(0,1.1*max(y),100); line(ygrid,gppdf(ygrid,kHat,sigmaHat)); xlim([0,6]); xlabel('Exceedance'); ylabel('Probability Density'); 

Hemos utilizado un ancho de contenedor bastante pequeño, por lo que hay una buena cantidad de ruido en el histograma. Aún así, la densidad ajustada sigue la forma de los datos, por lo que el modelo GP parece ser una buena opción.

También podemos comparar el CDF empírico con el CDF ajustado.

[F,yi] = ecdf(y); plot(yi,gpcdf(yi,kHat,sigmaHat),'-'); hold on; stairs(yi,F,'r'); hold off; legend('Fitted Generalized Pareto CDF','Empirical CDF','location','southeast'); 

Errores estándar de cálculo para las estimaciones de parámetros

Para cuantificar la precisión de las estimaciones, usaremos los errores estándar calculados a partir de la matriz de covarianza asintótica de los estimadores de máxima verosimilitud. La función calcula, como su segunda salida, una aproximación numérica a esa matriz de covarianza.gplike Alternativamente, podríamos haber llamado con dos argumentos de salida, y habría devuelto los intervalos de confianza para los parámetros.gpfit

[nll,acov] = gplike(paramEsts, y); stdErr = sqrt(diag(acov)) 
 stdErr =      0.1158     0.1093  

Estos errores estándar indican que la precisión relativa de la estimación para k es bastante inferior a la de Sigma: su error estándar está en el orden de la estimación. Los parámetros de forma suelen ser difíciles de estimar. Es importante tener en cuenta que el cómputo de estos errores estándar supuso que el modelo GP es correcto, y que tenemos suficientes datos para la aproximación asintótica a la matriz de covarianzas para sostener.

Comprobando la asunción de normalidad asintótica

Por lo general, la interpretación de los errores estándar implica suponer que, si el mismo ajuste se repite muchas veces en los datos procedentes de la misma fuente, las estimaciones de máxima verosimilitud de los parámetros seguirán aproximadamente una distribución normal. Por ejemplo, los intervalos de confianza se basan a menudo en esta suposición.

Sin embargo, esa aproximación normal puede o no ser buena. Para evaluar lo bueno que es en este ejemplo, podemos usar una simulación de arranque. Vamos a generar 1000 replicar datasets remuestreando de los datos, ajustar una distribución GP a cada uno, y guardar todas las estimaciones de réplica.

replEsts = bootstrp(1000,@gpfit,y); 

Como una comprobación aproximada de la distribución de muestreo de los estimadores de parámetros, podemos ver histogramas de las réplicas de bootstrap.

subplot(2,1,1); hist(replEsts(:,1)); title('Bootstrap estimates of k'); subplot(2,1,2); hist(replEsts(:,2)); title('Bootstrap estimates of sigma'); 

Mediante una transformación de parámetros

El histograma de las estimaciones de bootstrap para k parece ser sólo un poco asimétrico, mientras que para las estimaciones de Sigma definitivamente aparece sesgada a la derecha. Un remedio común para esa asimetría es estimar el parámetro y su error estándar en la escala logaritmo, donde una aproximación normal puede ser más razonable. Una gráfica Q-Q es una mejor manera de evaluar la normalidad que un histograma, porque la no normalidad se muestra como puntos que no siguen aproximadamente una línea recta. Vamos a comprobar que para ver si la transformación de registro para Sigma es apropiada.

subplot(1,2,1); qqplot(replEsts(:,1)); title('Bootstrap estimates of k'); subplot(1,2,2); qqplot(log(replEsts(:,2))); title('Bootstrap estimates of log(sigma)'); 

Las estimaciones de bootstrap para k y log (Sigma) parecen aceptablemente cercanas a la normalidad. Una gráfica Q-Q para las estimaciones de Sigma, en la escala sin registrar, confirmaría la asimetría que ya hemos visto en el histograma. Por lo tanto, sería más razonable construir un intervalo de confianza para Sigma calculando primero uno para log (Sigma) bajo el supuesto de normalidad, y luego exponenciación para transformar ese intervalo de nuevo a la escala original para Sigma.

De hecho, eso es exactamente lo que hace la función detrás de las escenas.gpfit

[paramEsts,paramCI] = gpfit(y); 
kHat kCI  = paramCI(:,1) 
 kHat =      0.0987   kCI =     -0.1283     0.3258  
sigmaHat sigmaCI  = paramCI(:,2) 
 sigmaHat =      0.7156   sigmaCI =      0.5305     0.9654  

Observe que mientras que el intervalo de confianza de 95% para k es simétrico sobre la estimación de máxima verosimilitud, el intervalo de confianza para Sigma no lo es. Esto se debe a que se creó transformando un CI simétrico para log (Sigma).