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 generalizada de Pareto

Este ejemplo muestra cómo ajustar los datos de cola a la distribución de Pareto generalizada por estimación de máxima verosimilitud.

Ajustar una distribución paramétrica a los datos a veces da como resultado un modelo que está bien de acuerdo con los datos en regiones de alta densidad, pero mal en áreas de baja densidad. Para distribuciones unimodales, como la normal o la t de Student, estas regiones de baja densidad se conocen como las "colas" de la distribución. Una razón por la que un modelo podría encajar mal en las colas es que, por definición, hay menos datos en las colas en los que basar una elección de modelo, 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, ajustar los datos en la cola es la principal preocupación. La distribución generalizada de Pareto (GP) se desarrolló como una distribución que puede modelar colas de una amplia variedad de distribuciones, basándose en argumentos teóricos. Un enfoque para el ajuste de distribución que implica el GP es utilizar un ajuste no paramétrico (la función de distribución acumulativa empírica, por ejemplo) en regiones donde hay muchas observaciones, y ajustar el GP a la cola 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'}); 

Observe 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. Además, el GP se utiliza a menudo junto con un tercer parámetro de umbral que desplaza el límite inferior lejos de cero. No necesitaremos esa generalidad aquí.

La distribución 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 superación

La distribución de GP se puede definir de forma constructiva en términos de excedencias. Comenzando con una distribución de probabilidad cuya cola derecha cae a cero, como la normal, podemos muestrear valores aleatorios independientemente de esa distribución. Si fijamos un valor de umbral, tiramos todos los valores que están por debajo del umbral y restamos el umbral de los valores que no se tiran, el resultado se conoce como superaciones. La distribución de las superaciones es aproximadamente un GP. Del mismo modo, 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 shape, k, de la distribución 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 el modelado de extremos de los retornos del mercado de valores y el modelado de inundaciones extremas. En este ejemplo, usaremos datos simulados, generados a partir de la distribución t de un estudiante con 5 grados de libertad. Tomaremos el 5% más grande de 2000 observaciones de la distribución t, y luego restamos el cuantil del 95% para obtener superaciones.

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 con la máxima verosimilitud

La distribución GP se define para 0 < sigma y -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 colas de ajuste de distribuciones como la beta o triangular, por lo que no presentarán ningún problema aquí.

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

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

Comprobación del ajuste visualmente

Para evaluar visualmente lo bueno que es el ajuste, trazaremos un histograma escalado de los datos de cola, superpuesto con la función de densidad del GP que hemos estimado. El histograma se escala de modo que las alturas de la barra multiplican su anchura 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 bandeja bastante pequeño, por lo que hay una buena oferta de ruido en el histograma. Aun 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 instalado.

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

Cálculo de errores estándar para las estimaciones de parámetros

Para cuantificar la precisión de las estimaciones, usaremos 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 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 menor que la de sigma -- su error estándar está en el orden de la estimación en sí. Los parámetros de forma suelen ser difíciles de estimar. Es importante tener en cuenta que el cálculo de estos errores estándar supone que el modelo GP es correcto y que tenemos suficientes datos para la aproximación asintótica a la matriz de covarianza que se va a contener.

Comprobación de la asunción de normalidad asintótica

La interpretación de los errores estándar suele implicar asumir que, si el mismo ajuste pudiera repetirse muchas veces en los datos procedentes de la misma fuente, las estimaciones de máxima probabilidad de los parámetros seguirían aproximadamente una distribución normal. Por ejemplo, los intervalos de confianza a menudo se basan 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. Generaremos 1000 conjuntos de datos de réplica mediante el remuestreo a partir de los datos, ajustaremos una distribución GP a cada uno y guardaremos 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 arranque.

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

Uso de una transformación de parámetros

El histograma de las estimaciones de arranque para k parece ser sólo un poco asimétrico, mientras que para las estimaciones de sigma definitivamente aparece sesgado a la derecha. Un remedio común para esa asimetría es estimar el parámetro y su error estándar en la escala de registro, 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 aparece 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 adecuada.

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 arranque para k y log(sigma) aparecen aceptablemente cerca de 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 la suposición de normalidad, y luego exponente para transformar ese intervalo de nuevo a la escala original para sigma.

De hecho, eso es exactamente lo que hace la función entre bastidores.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 del 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).