MedialdeaVillanueva TFM

Descargar como pdf o txt
Descargar como pdf o txt
Está en la página 1de 65

Análisis de Valores

Extremos
Modelización Espacial

Autora: Adriana Medialdea Villanueva


Dirigido por: José Miguel Angulo Ibáñez

Máster en Estadı́stica Aplicada


Facultad de Ciencias. Universidad de Granada.

12 de Septiembre de 2016
Índice general

Resumen 5

1. Teorı́a de valores extremos univariante 9


1.1. Distribución de valores extremos generalizada . . . . . . . . . . . 9
1.2. Modelos de máximos por bloques . . . . . . . . . . . . . . . . . . 11
1.2.1. Estimación de los parámetros de la distribución GEV . . 12
1.2.2. Niveles de retorno . . . . . . . . . . . . . . . . . . . . . . 12
1.3. Modelos de excedencias de umbrales . . . . . . . . . . . . . . . . 13
1.3.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3.2. Distribución de Pareto generalizada . . . . . . . . . . . . 13
1.3.3. Estimación de los parámetros de la distribución de Pareto
generalizada . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.3.4. Selección del umbral . . . . . . . . . . . . . . . . . . . . . 17
1.3.5. Niveles de retorno . . . . . . . . . . . . . . . . . . . . . . 19
1.4. Modelización de extremos en secuencias de datos dependientes . 20

2. Teorı́a de valores extremos multivariante 23


2.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.2. Cópulas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3. Distribución de valores extremos multivariante . . . . . . . . . . 24
2.3.1. Caso bivariante . . . . . . . . . . . . . . . . . . . . . . . . 24
2.3.2. Extensión para n variables . . . . . . . . . . . . . . . . . . 26
2.3.3. Transformaciones . . . . . . . . . . . . . . . . . . . . . . . 26
2.4. Convergencia del vector de máximos . . . . . . . . . . . . . . . . 27
2.5. Funciones de dependencia . . . . . . . . . . . . . . . . . . . . . . 28
2.6. Modelización y estimación de parámetros . . . . . . . . . . . . . 30
2.6.1. Caso bivariante . . . . . . . . . . . . . . . . . . . . . . . . 30

3
4 ÍNDICE GENERAL

2.6.2. Extensión para n variables. Estimación mediante cópulas 31


2.7. Niveles de retorno . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.7.1. Caso bivariante . . . . . . . . . . . . . . . . . . . . . . . . 31
2.7.2. Extensión para n variables. Medida de intensidad . . . . . 31
2.8. Modelos de excedencias de umbrales . . . . . . . . . . . . . . . . 32
2.8.1. Caso bivariante . . . . . . . . . . . . . . . . . . . . . . . 32
2.8.2. Extensión a n variables . . . . . . . . . . . . . . . . . . . 34

3. Extremos espaciales 37
3.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2. Procesos max-estables espaciales . . . . . . . . . . . . . . . . . . 37
3.2.1. Representación espectral de los procesos max-estables . . 38
3.2.2. Modelos max-estables . . . . . . . . . . . . . . . . . . . . 39
3.3. Medidas de dependencia para extremos espaciales . . . . . . . . . 40
3.3.1. Función del coeficiente extremal . . . . . . . . . . . . . . 40
3.3.2. F-madograma . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.3.3. λ-madograma . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.4. Simulación de procesos espaciales max-estables . . . . . . . . . . 42
3.4.1. Simulaciones no condicionadas . . . . . . . . . . . . . . . 42
3.4.2. Simulaciones condicionadas . . . . . . . . . . . . . . . . . 42
3.5. Ajuste de un proceso max-estable Fréchet a los datos . . . . . . . 43
3.5.1. Mı́nimos cuadrados . . . . . . . . . . . . . . . . . . . . . . 43
3.5.2. Máxima verosimilitud . . . . . . . . . . . . . . . . . . . . 44
3.6. Selección del modelo . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.6.1. Criterio de información de Takeuchi . . . . . . . . . . . . 44
3.6.2. Estadı́stico de la tasa de verosimilitud . . . . . . . . . . . 45

4. Aplicación con R. Extremos espaciales 47


4.1. Ajuste de la distribución de valores extremos generalizada a los
datos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.2. Transformación de los datos a la distribución Fréchet. . . . . . . 52
4.3. Ajuste del proceso max-estable a los datos transformados. . . . . 53
4.3.1. Modelo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.3.2. Modelo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.3.3. Diagnosis y elección del modelo . . . . . . . . . . . . . . . 55
4.3.4. Predicciones. . . . . . . . . . . . . . . . . . . . . . . . . . 59
ÍNDICE GENERAL 5

5. Conclusión 61
6 ÍNDICE GENERAL
Resumen

La Teorı́a de Valores Extremos es una displicina estadśtica cuyos avances fun-


damentales son relativamente recientes. Está enfocada al análisis del comporta-
miento estocástico de los valores extremos de un proceso.

Con amplias aplicaciones en hidrologı́a, investigación medioambiental y me-


teorologı́a, seguros, aplicaciones financieras y geologı́a entre otras, tiene como
objetivo cuantificar el comportamiento de una serie de valores de un proceso,
siendo estos valores más extremos que los que usualmente se observan en él.

Esta teorı́a se desarrolla a partir de los resultados de Fréchet (1927), Fisher y


Tippett (1928) y von Mises (1936), que sientan las bases de la teorı́a probabilis-
tica de valores extremos unidimensional; sin embargo, su desarrollo teórico no
tendrá lugar hasta principios de los años 70 con la tesis dotoral de de Haan (1970)
y los resultados de inferencia estadı́stica desarrollados por Pickands (1975); estos
resultados supusieron las primeras contribuciones a la Teorı́a de Valores Extre-
mos Multivariante y motivaron el desarrollo de modelos alternativos basados en
excendencias de umbrales.

En la actualidad, los avances en esta teorı́a se está centrando principalmente en


el desarrollo de modelos y métodos para valores extremos de fenómenos espa-
ciales y otras estructuras más complejas.

En el presente trabajo se realiza una revisión de los principales enfoques existen-


tes en el análisis de valores extremos; prestando especial atención a la dimensión
espacial, junto con una aplicación práctica a datos geoestadı́sticos.

El primer capı́tulo está dedicado a la teorı́a univariante, donde se presenta el


punto de partida de esta disciplina a partir de la condición de valores extremos;
además de los modelos clásicos de máximos por bloques y excedencias de umbra-
les junto con las distintas distribuciones a las que dan lugar cada uno de ellos.
Finalmente, se introduce la modelización para el caso de datos dependientes.

El tercer capı́tulo está enfocado a las distribuciones de valores extremos multiva-


riantes y su aproximación para maximos por bloques y excedecias de umbrales,
prestando especial atención a su estructura de dependencia a través de cópulas.
En el cuarto capı́tulo se trata la modelización de la dependencia tras la inclu-
sión de una dimensión espacial, para lo que serán fundamentales los procesos
max-estables y las distribuciones finito-dimensionales que se derivarán a través
de ellos.

Finalmente, en el capı́tulo cinco se realiza una aplicación práctica con el software


R y la herramienta SpatialExtremes.
Capı́tulo 1

Teorı́a de valores extremos


univariante

1.1. Distribución de valores extremos generali-


zada
El modelo para el que se desarrolla la teorı́a de valores extremos está enfocado
a describir el comportamiento estadı́stico de

Mn = máx{X1 , ..., Xn },

donde X1 , ..., Xn es una secuencia de variables aleatorias independientes con


distribución común F y Mn representa el máximo del proceso sobre n unidades
de tiempos de observación.
La distribución de Mn podrı́a obtenerse de manera exacta a partir de la distri-
bución de las n variables, teniendo en cuenta las propiedades de independencia:

FMn (z) = P r[Mn 6 z] = P r[X1 6 z, ..., Xn 6 z]


(1.1)
= P r[X1 6 z] × ... × P r[Xn 6 z] = [F (z)]n .

y derivando, se obtendrı́a su función de densidad:

n−1
fMn (z) = n (FMn (z)) f (z) . (1.2)

La función de distribución anteriormente calculada converge a cero cuando n →


∞ para z > z ∗ y a uno para z ≤ z ∗ , con z ∗ = sup{z : F (z) < 1}. Por lo que,
para obtener una distribución lı́mite no degenerada, será necesario llevar a cabo
una normalización.
Esta solución, basada en el Teorema Central del Lı́mite, consiste en la búsqueda
de secuencias de constantes {bn ; n ≥ 1} y {an ; n ≥ 1} tales que la distribución
de

9
Mn − bn
Mn∗ = (1.3)
an

converge a una distribución no degenerada cuando n → ∞, es decir,

lı́m F n (an z + bn ) = G (z) . (1.4)


n→∞

El rango completo de distribuciones lı́mite que podrá seguir Mn∗ vendrá dado
por el Teorema de Valores Extremos:

Teorema 1.1.1 Si existen sucesiones de constantes {an > 0} y {bn } tales que
 
M n − bn
P ≥z → G(z), cuando n → ∞,
an

siendo G una función de distribución no degenerada, entonces G debe pertenecer


a una de las siguientes familias:

i. Gumbel: G(z) = exp − exp − z−b


  
a , −∞ < z < ∞
(
0, n z6b
ii. Fréchet: G(z) = z−b −α
 o
exp − a , z>b
( n  α o
z−b
exp − − a , z<b
iii. Weibull: G(z) =
1, z≥b

con parámetros a > 0, b y, en el caso de las familias ii y iii, α > 0.

Las anteriores distribuciones se conocen como distribuciones de valores ex-


tremos y serán las únicas a las que pueda converger la variable Mn∗ , indepen-
dientemente de como se distribuya F .
Las familias de distribuciones Fréchet, Gumbel y Weibull se pueden combinar
en una única familia con función de distribución
(   −1/ξ )
z−µ
G(z) = exp − 1 + ξ (1.5)
σ

definida en {z : 1 + ξ (z − µ) /σ > 0}, donde los parámetros de localización,


escala y forma satisfacen, respectivamente, −∞ < µ < ∞, σ > 0 y −∞ < ξ <
∞. Esta familia de distribuciones se conoce como familia de distribuciones
de valores extremos generalizada (GEV) (o familia de distribuciones de
Von Mises).
La especificación de ξ determinará el comportamiento de la cola de la distribu-
ción, de forma que según el valor que tome este parámetro se tendrá una de las
siguientes distribuciones:
Gumbel si ξ = 0
Fréchet si ξ > 0
Weibull si ξ < 0

De esta forma, el teorema 1.1.1 se replanteará como sigue:

Teorema 1.1.2 Si existen sucesiones de constantes {an > 0} y {bn } tales que

P r{(Mn − bn ) /an 6 z} → G(z), cuando n → ∞, (1.6)

para una distribución G no degenerada, entonces G pertenece a la familia de


distribuciones GEV,
(   −1/ξ )
z−µ
G(z) = exp − 1 + ξ . (1.7)
σ

El concepto de max-estabilidad, que se introduce a continuación, está estrecha-


mente relacionado con el teorema anterior.

Definición 1.1.1 Una distribución G se denomina max-estable si, para todo


n = 2, 3..., existen constantes αn > 0 y βn tales que

Gn (αn z + βn ) = G (z) . (1.8)

Es decir, la propiedad de max-estabilidad la sastisfarán aquellas distribuciones


en las que la operación de tomar máximos muestrales conduzca a una distribu-
ción idéntica, aunque con distintos parámetros de localización y escala.
Esta propiedad se relaciona con el Teorema de Valores Extremos a partir del
siguiente resultado:

Teorema 1.1.3 Una distribución es max-estable si y solo si es una distribución


de valores extremos generalizada.

1.2. Modelos de máximos por bloques


La familia de distribuciones GEV será útil para modelizar la distribución de los
máximos por bloques.
El procedimiento consiste en agrupar los datos en bloques de igual tamaño y,
a continuación, ajustar la distribución GEV al conjunto de los máximos corres-
pondientes a cada uno de los bloques.
El principal problema que presenta este método reside en la elección del tamaño
de los bloques, para la cual habrá que encontrar un equilibrio entre el sesgo y
la varianza.
La elección de bloques muy pequeños conducirá a una pobre aproximación del
modelo, con lo que se aumentará el sesgo al estimar y extrapolar. Por el con-
trario, la elección de bloques demasiado grandes aumentará la varianza de las
estimaciones.
Por cuestiones prácticas, en secuencias de datos temporales mensuales, se suelen
tomar bloques de longitud anual, de esta manera los máximos se distribuyen de
manera similar en cada uno de los bloques.

1.2.1. Estimación de los parámetros de la distribución GEV


Sean Z1 , ..., Zk v.a.i.i.d. con distribución GEV. La función de log-verosimilitud
para ξ 6= 0 es

k   
X zi − µ
log L(µ, σ, ξ) = − k log σ − (1 + 1/ξ) log 1 + ξ
i=1
σ
k   −1/ξ
X zi − µ
− 1+ξ ,
i=1
σ

con
 
zi − µ
1+ξ > 0, para i = 1, ..., k, (1.9)
σ

mientras que para el caso ξ = 0,

k   k   
X zi − µ X zi − µ
log L (µ, σ) = −k log σ − − exp − . (1.10)
i=1
σ i=1
σ

La maximización de estas ecuaciones no tiene solución analı́tica; sin embargo,


para un conjunto dado de datos, la maximización se obtiene de manera directa
usando algoritmos numéricos de optimización.

1.2.2. Niveles de retorno


Se denomina niveles de retorno a los cuantiles de la distribución de valores
extremos generalizada; se notarán como zp y se obtendrán como la inversa de
la función de distribución GEV:

 µ − σ 1 − {− log (1 − p)}−ξ ,
 h i
para ξ 6= 0
zp = ξ (1.11)
µ − σ log {− log (1 − p)} , para ξ = 0,

de forma que zp es el nivel de retorno asociado al periodo de retorno 1/p, es decir,
se espera que el nivel zp sea excedido en promedio una vez cada 1/p unidades
de tiempo. Equivalentemente, zp será excedido en una unidad de tiempo con
probabilidad p.
La inferencia de los niveles de retorno, para una probabilidad p dada, se llevará a
cabo por sustitución directa de las estimaciones de los parámetros del modelo.

1.3. Modelos de excedencias de umbrales

1.3.1. Introducción
Modelizar únicamente máximos de bloques es una aproximación poco eficiente
en el ánalisis de valores extremos si algunos de los bloques contienen eventos
más extremos que el resto. En este caso, será adecuado el uso de modelos de
excedencias de umbrales.
Sea X1 , X2 , ... una secuencia de variables aleatorias independientes e idéntica-
mente distribuidas con función de distribución marginal F . Se considerará un
evento extremo aquel que sobrepase el valor de un umbral u y se notarán las
excendencias de este umbral como Y = X − u.
De esta forma, el comportamiento de las excedencias de X sobre el umbral u
vendrá dado por la probabilidad condicionada:

1 − F (u + y)
P r{X > u + y|X > u} = , y > 0. (1.12)
1 − F (u)

Si se conociera la distribución F , la distribución de las excedencias de umbral


serı́a igualmente conocida; sin embargo, esto no sucede en la práctica, por lo
que será necesario el uso de la distribución GEV como aproximación.

1.3.2. Distribución de Pareto generalizada


La distribución de Pareto generalizada se utilizará como aproximación para la
distribución lı́mite de las excedencias de umbrales.

Teorema 1.3.1 Sea X1 , X2 , ... una secuencia de variables aleatorias indepen-


dientes con distribución común F, y sea

Mn = máx{X1 , ..., Xn }. (1.13)

Notaremos por X a un término arbitrario Xi de esta secuencia y se supondrá que


F satisface el teorema 1.1.2, es decir,

P r{Mn 6 z} ≈ G (z) , cuando n → ∞, (1.14)


donde
(   −1/ξ )
z−µ
G (z) = exp − 1 + ξ (1.15)
σ
para µ, σ > 0 y ξ. Entonces, para un umbral u suficientemente grande, la función
de distribución de (X − u), condicionada a X > u, tendrá la forma
  −1/ξ
 ξy
1 − 1 +
 , y>0 si ξ 6= 0
H (y) = σ̃ (1.16)
  y
1 − exp − , y>0 si ξ = 0,


σ̃
  ξ1
ξy
estando definida en {y : y > 0 y 1+ σ̃ > 0}, donde

σ̃ = σ + ξ (u − µ) . (1.17)

Esta familia se denomina familia de distribuciones de Pareto generaliza-


da, y viene caracterizada por los parámetros de escala σ̃ y de forma −∞ < ξ <
+∞.
Según el teorema anterior, si los máximos por bloques siguen una distribución
G, entonces la distribución de las excedencias de umbral se encuentra dentro de
la familia de distribuciones de Pareto generalizada.
El comportamiento de la distribución de Pareto generalizada está determinado
por el parámetro ξ:

Si ξ < 0, la distribución de los excesos tiene como lı́mite superior u − σ̃/ξ.


Si ξ ≥ 0, la distribución no tiene lı́mite superior. En concreto, para el caso
ξ > 0 se tiene la distribución de Pareto ordinaria.

En el caso en el que ξ = 0, la distribución se corresponde con una exponencial


de parámetro 1/σ̃.
A partir de este punto, por simplicidad en la notación, el parámetro de escala
de la distribución de Pareto Generalizada, σ̃, pasará a notarse como σ.

1.3.3. Estimación de los parámetros de la distribución de


Pareto generalizada
Dado un valor del umbral u y el número de datos, k, de la muestra original
X1 , ..., Xn que sobrepasan el umbral, la estimación de los parámetros σ y ξ se
puede llevar a cabo mediante diferentes métodos.

Método de máxima verosimilitud

El logaritmo de la función de verosimilitud para una muestra y1 , ..., yk de varia-


bles aleatorias independientes e idénticamente distribuidas con distribución de
Pareto generalizada viene dado por
 Xk  
1 ξyi
log L (σ, ξ) = −k log σ − +1 log 1 + , si ξ 6= 0 (1.18)
ξ i=1
σ

k
1X
log L (σ, ξ) = −k log σ − yi , si ξ = 0. (1.19)
σ i=1

Para maximizar esta función; en el caso en que ξ 6= 0, será necesario llevar a


cabo la siguiente reparametrización:

σ
(σ, ξ) → (τ, ξ) con τ = ,
ξ

con lo que se tiene

 Xk
1
log L (τ, ξ) = −k log ξ + k log τ − +1 log (1 + τ yi ) .
ξ i=1

ML
Los estimadores máximo-verosı́miles τ̂k,n y ξˆk,n
ML
se obtienen a partir de la ecua-
ción
! k
1 1 1X yi
ML
− +1 M Ly
= 0,
τ̂k,n ξˆk,n
ML k i=1 1 + τ̂k,n i

donde

k
1X
ξˆk,n
ML ML

= log 1 + τ̂k,n yi ,
k i=1

Para la obtención de las estimaciones de los parámetros será necesario utili-


zar técnicas numéricas, debido a que la maximización de la función de log-
verosimilitud no se puede llevar a cabo de forma analı́tica.

Método de los momentos ponderados

Los momentos ponderados de una variable aleatoria Y con función de distribu-


ción F se definen como

r s
Mp,r,s = E {Y p [F (Y )] [1 − F (Y )] } ,

para p, r y s números reales. Para el caso de la distribución de Pareto generali-


zada se considerará p = 1, r = 0 y s = 0, 1, 2, ..., obteniéndose
σ
M1,0,s = , para ξ < 1, (1.20)
(s + 1) (s + 1 + ξ)

siendo el equivalente muestral

k  s
1X j
M̂1,0,s = 1− yj,k .
k j=1 k+1

Los estimadores de los parámetros se obtienen resolviendo (1.20) para s = 0 y


s = 1:

M̂1,0,0
ξˆP W M = 2 − ,
M̂1,0,0 − 2M̂1,0,1
2M̂1,0,0 M̂1,0,1
σ̂ = .
M̂1,0,0 − 2M̂1,0,1

La aplicación de este método presenta algunos problemas. En el caso en que


ξ ≥ 1, los momentos no existen; y para el caso en el que ξ < 0, las estimaciones
pueden ser inconsistentes con los datos observados.

Método del percentil básico

Este método es válido en el caso en que ξ 6= 0; en el caso contrario, el parámetro


σ se podrá estimar de forma eficiente mediante el método de máxima verosimi-
litud.
Se considerarán dos estadı́sticos de orden, Yi,k y Yj,k , correspondientes a la
muestra Y1 , ..., Yk . La función de distribución acumulativa de una distribución
de Pareto generalizada, evaluada en estos estadı́sticos para sus correspondien-
tes valores percentuales, da como resultado un sistema de ecuaciones con dos
incógnitas:

− 1
1 − (1 + τ̂i,j yi,k ) = pi,n , (1.21)
ξˆi,j

− 1
1 − (1 + τ̂i,j yj,k ) = pj,n , (1.22)
ˆ
ξi,j

i
donde τ = ξ/σ y pi,n = n+1 .
Resolviendo este sistema, se obtiene

log (1 + τ̂i,j yi,k )


ξˆi,j = ,
− log (1 − pi,n )

ξˆi,j
σ̂i,j = .
τ̂i,j
Los estimadores finales se obtienen calculando ξˆi,j y σ̂i,j para todos los pares de
estadı́sticos de orden Yi,k < Yj,k :
n o
ξˆEP M = mediana ξˆi,j ; i<j ,

σ̂EP M = mediana {σ̂i,j ; i < j} .

1.3.4. Selección del umbral


Para identificar eventos extremos se definirá un umbral u; los eventos que su-
peren este umbral, {xi : xi > u}, se etiquetarán como x(1) , ..., x(k) , y se definirán
las excedencias de umbrales yj = x(j) − u, para j = 1, ..., k.
Teniendo en cuenta el teorema (1.3.1), los valores yj se considerarán realizacio-
nes independientes de una variable aletoria cuya distribución se podrá aproximar
mediante la familia de Pareto generalizada.

Para la elección del umbral habrá que tener en cuenta que una elección de
un umbral muy bajo aumentarı́a el sesgo del modelo, mientras que la elección
de un umbral demasiado alto generarı́a un número de excedencias muy pequeño,
con lo que aumentarı́a su varianza.

Se dispone de dos métodos para la elección del umbral:

Método 1: Este método se basa en la esperanza de la distribución de Pareto


generalizada, de forma que si Y sigue esta distribución con parámetros σ y ξ,

σ
E(Y ) = , con ξ > 1. (1.23)
1−ξ
Si suponemos que la distribución de Pareto generalizada es válida para modeli-
zar las excedencias de un cierto umbral al que notaremos u0 , para un término
arbitrario X de la serie se tendrá:

σu0
E (X − u0 |X > u0 ) = , con ξ > 1,
1−ξ
donde σu0 será el parámetro de escala de la distribución de las excedencias del
umbral u0 . Es decir, si la distribución de Pareto es válida para las excedencias
del umbral u0 , será igualmente válida para cualquier umbral u > u0 mediante
una sustitución del paramétro de escala por σu . De esta forma, y teniendo en
cuenta (1.17) , se tiene

σu σu + ξu
E (X − u|X > u) = = 0 . (1.24)
1−ξ 1−ξ
Se espera que los valores proporcionados por esta expresión varı́en linealmen-
te al variar u en los niveles en los que el ajuste de la distribución de Pareto
generalizada es apropiada, lo que se traduce en el siguiente procedimiento.
Se representarán los puntos
( nu
! )
1 X 
u, x(i) − u : u < xmax ,
nu i=1

donde x(1) , ..., x(nu ) son las observaciones que exceden u, y xmax es el máximo
Xi . Este gráfico se denomina gráfico de vida media residual.
Se tomará como umbral el valor máximo para el que el gráfico muestre una
tendencia lineal.

Método 2 Consiste en ajustar la distribución de Pareto generalizada en un


rango de umbrales y elegir el umbral adecuado según la estabilidad de los
parámetros estimados.

Según la caracterización de la familia de distribuciones de Pareto generalizada,


si una distribución de Pareto generalizada es útil para modelizar las excedencias
de un umbral u0 , las excedencias de un umbral u más alto que el anterior tam-
bién seguirán una distribución de Pareto generalizada, siendo sus parámetros de
forma idénticos, mientras que sus parámetros de escala satisfacen la siguiente
relación:

σu = σu0 + ξ (u − u0 ) si ξ 6= 0. (1.25)

En el caso en que ξ = 0 se realizará la siguiente reparametrización:

σ ∗ = σu − ξu.

Se representarán σ̂ ∗ y ξˆ frente a u. Las estimaciones de ambos parámetros


deberı́an ser constantes por encima de u0 si u0 es un umbral válido, para excen-
dencias que siguen la distribución de Pareto generalizada.
Los intervalos de confianza para ξˆ se obtendrán a través de la matriz de varianzas-
covarianzas V , mientras que los intervalos de confianza para σ̂ ∗ requieren el uso
del siguiente resultado:

Teorema 1.3.2 (Método Delta). Sea θ̂ el estimador máximo-verosı́mil de θ


con matriz de varianzas-covarianzas Σθ̂ . Sea ψ = h (θ) un vector de parámetros
de dimensión r definido como función de θ. Entonces,

ψ̂ ∼ Nr ψ; 5Tθ ψΣθ̂ 5θ ψ ,


donde 5θ ψ es la matriz de orden k × r de las derivadas parciales de ψ respecto


θ, dada por

∂ψ1 ∂ψ2 ∂ψr


 
∂θ1 ∂θ1 ··· ∂θ1
∂ψ1 ∂ψ2 ∂ψr
···
 
 ∂θ2 ∂θ2 ∂θ2 
5θ ψ = 
 .. .. .. .
.. 
 . . . . 
∂ψ1 ∂ψ2 ∂ψr
∂θk ∂θk ··· ∂θk
Por tanto, se tendrı́a que

V ar (σ ∗ ) ≈ 5σ ∗T V σ ∗ ,

donde

∂σ ∗ ∂σ ∗
 
∗T
5σ = , = [1, −u] .
∂σu ∂ξ

1.3.5. Niveles de retorno

Se denomina niveles de retorno, xk , a los cuantiles de la distribución de valores


extremos, con 1/k la probabilidad de que xk sea superado al menos una vez en
una unidad de tiempo.
Suponemos que las excedencias de un umbral u siguen una distribución de Pareto
generalizada; es decir, para x > u,

  −1/ξ
x−u
P r {X > x|X > u} = 1 + ξ . (1.26)
σ

Por lo que, notando ζu = P r [X > u], el nivel xk que se excede en promedio


cada k observaciones es la solución de la ecuación

  −1/ξ
xk − u 1
ζu 1 + ξ = . (1.27)
σ k

Despejando el término xk , se tiene que el nivel de retorno para excedencias de


umbrales viene dado por:

 u + σ (kζu )ξ − 1
 h i
para ξ 6= 0,
xk = ξ (1.28)
u + σ log (kζu ) para ξ = 0.

Para realizar la estimación de los niveles de retorno se sustituirán los parámetros


por sus estimaciones máximo verosı́miles.
Un concepto relacionado con los niveles de retorno es el tiempo de retorno, o
los tiempos en los que se superarı́a un cierto umbral; estos tiempos se definen
recursivamente como:

τ1 = mı́n {k : Xk > u} ,

τr = mı́n {k > τr−1 : Xk > u} , r > 1.


1.4. Modelización de extremos en secuencias de
datos dependientes
Hasta ahora se ha tenido en cuenta que los datos se habı́an obtenido a partir
de sucesiones de variables aleatorias independendientes; sin embargo, en la ma-
yorı́a de casos prácticos esta asunción no es realista. En estos casos, los usual
será considerar que la serie es estacionaria.
Una serie se considerará estacionaria si su distribuciones de probabilidad se man-
tienen estables a lo largo del tiempo. En este caso, las variables serán mutua-
mente dependientes, pero sus propiedades estocásticas permanecerán constantes
a lo largo del tiempo.
La estacionariedad se define de como:

Teorema 1.4.1 Sea X1 , X2 , ... una serie de tiempo. Esta serie se dice esta-
cionaria si para cada conjunto de ı́ndices 1 ≤ t1 < t2 < · · · < tm , la distri-
bución conjunta de (Xt1 , Xt2 , ..., Xtm ) coincide con la distribución conjunta de
(Xt1 +h , Xt2 +h , ..., Xtm +h ), para cualquier h.

Coles (2001) da la siguiente definición:

Definición 1.4.1 Una serie estacionaria X1 , X2 , ... satisface la condición de


los D (un ) si, para todo i1 < ... < ip < ji < ... < jq con ji − ip > l, se tiene que

|P [Xi1 ≤ un , ..., Xip ≤ un , Xj1 ≤ un , ..., Xjq ≤ un ]


− P [Xi1 ≤ un , ..., Xip ≤ un ] P [Xj1 ≤ un , ..., Xjq ≤ un ] | ≤ α (n, 1) ,
ln
donde α (n, ln ) → 0 para alguna sucesión ln de forma que n → 0 cuando
n → ∞.

Para secuencias de variables independientes la diferencia en probabilidad ante-


rior es igual a cero. En el caso en que las variables no son independientes, si
están suficientemente alejadas entre sı́ la diferencia en probabilidad es lo su-
ficientemente cercana a cero para no tener efecto en las leyes de lı́mites para
extremos, lo que se resume en el siguiente teorema:

Teorema 1.4.2 Sea X1 , X2 , ... un proceso estacionario y sea Mn = máx {X1 , ..., Xn }.
Entonces, si {an > 0} y {bn } son secuencias de constantes tales que

P r {(Mn − bn ) /an ≤ z} → G (z) ,


donde G es una función de distribución no degenerada y la condición D (un )
se satisface con un = an z + bn para cualquier número real z, G pertenece a la
familia de distribuciones de valores extremos generalizada.

Este resultado implica que dada una sucesión de variables suficientemente leja-
nas y dependientes entre sı́ (en el sentido especificado por la condición D (un )),
los valores máximos de estas series estacionarias se distribuyen según la misma
distribución lı́mite que en el caso de sucesiones de variables independientes. Sin
embargo, los parámetros de la distribución se verán afectados por la dependen-
cia.

Teorema 1.4.3 Sea X1 , X2 ... un proceso estacionario y sea X1∗ , X2∗ , ... una su-
cesión de variables independientes con la misma distribución marginal. Se define
Mn = máx {X1 , ..., Xn } y Mn∗ = máx {X1∗ , ..., Xn∗ }. Bajo condiciones adecuadas
de regularidad,
P r {(Mn∗ − bn ) /an ≤ z} → G1 (z)
cuando n ← ∞ para sucesiones {an > 0} y {bn }, donde G1 es una función de
distribución no degenerada, si y solo si

P r {(Mn − bn ) /an ≤ z} → G2 (z) ,

donde
G2 (z) = Gθ1 (z)
para una constante θ tal que 0 < θ ≤ 1.

De esta forma, si Mn∗ sigue una distribución GEV de parámetros (µ, σ, ξ), Mn
seguirá una distribución GEV con parámetros (µ∗ , σ ∗ , ξ) donde

σ
µ∗ = µ − 1 − θ−ξ ,

ξ
σ ∗ = σθξ .
Capı́tulo 2

Teorı́a de valores extremos


multivariante

2.1. Introducción
En el estudio de extremos para dos o más procesos, si bien cada uno de ellos pue-
de ser modelizado individualmente usando técnicas univariantes, sin embargo,
también será interesante estudiar las relaciones que puedan existir entre ellos.
En particular, puede ocurrir que alguna combinación de estos procesos sea de
mayor interés que los procesos individuales.
Sean Xi = (Xi,1 , ..., Xi,p ) , i ∈ {1, ..., n}, una secuencia de vectores aleatorios in-
dependientes e idénticamente distribuidos con función de distribución conjunta
F y marginales F1 , ..., Fp . El vector de máximos, que se notará como M , para
el caso multivariante se define como:

Mj = máx {Xi,j } ; para j = 1, ..., p,


i∈{1,...,n}

M = {M1 , ..., Mp } .

Este vector no se tiene que corresponder necesariamente con un vector observado


en la serie original.

2.2. Cópulas
En esta sección se describe uno de los procedimientos más usados para repre-
sentar la estructura de dependencia existente entre varias funciones de densidad
univariantes.
A través las cópulas, se podrán analizar y estimar de forma separada las compo-
nentes univariantes y la estructura de dependencia multivariante. La unión de

23
las diferentes estimaciones conducirá a una estimación de la función de densidad
conjunta.
Sea X = (X1 , ..., Xp ) un vector aleatorio con función de distribución F .
Supongamos que cada una de las variables Xj es continua con función de dis-
tribución marginal Fj ; entonces, aplicando una transformación de probabilidad,
Yj = Fj (Xj ) tendrá distribución uniforme (0, 1).
La función de distribución C del vector Y = (Y1 , ..., Yp ) es la cópula de F . Es
decir,

C (u) = P (Y1 ≤ u1 , ..., Yp ≤ up ) = F F1−1 (u1 ) , ..., Fp−1 (up ) ,




donde Fj−1 es la función cuantil de Fj .


Se podrá volver a la distribución original desde la cópula mediante la transfor-
mación cuantil:

F (X) = C (F1 (x1 ) , ..., Fp (xp )) .

2.3. Distribución de valores extremos multiva-


riante

2.3.1. Caso bivariante


Coles (2001) propone una caracterización asintótica de la distribución de valores
extremos, a partir de la función de distribución conjunta de M ∗ = (M1 /n, M2 /n):

Teorema 2.3.1 Sea M ∗ = (M1∗ , M2∗ ), con (Xi1 , Xi2 ), i = 1, ..., n, vectores
independientes con distribución marginal Fréchet. Entonces, si para n → ∞
d
P r {M1∗ ≥ x1 , M2∗ ≥ x2 } −
→ G (x1 , x2 ) , (2.1)

donde G es una función de distribución no degenerada, se tiene que

G (x1 , x2 ) = exp {−V (x1 , x2 )} , x1 > 0, x2 > 0,

donde
1  
ω 1−ω
Z
V (x1 , x2 ) = 2 máx , dH (ω) (2.2)
0 x1 x2
y H es una función de distribución en [0, 1] que satisface la condición
Z 1
ωdH (ω) = 1/2.
0

A diferencia del caso univariante, esta distribución engloba a infinitas familias


de distribuciones.
Ejemplo 1

Si se distribuye la masa de la fución H en los puntos ω = 0 y ω = 1, se obtiene

V (x1 , x2 ) = x−1 −1
1 + x2

y la correspondiente distribución de valores extremos generalizada será

G (x1 , x2 ) = exp − x−1 −1


 
1 + x2 , x1 > 0, x2 > 0.

Esta función se puede factorizar en función de x1 y x2 , por lo que ambas variables


serı́an independientes.

Ejemplo 2

Si se concentra toda la masa de la función H en ω = 0,5, se obtiene,

G (x1 , x2 ) = exp − máx x−1 −1


 
1 , x2 , x1 > 0, x2 > 0,

que corresponderı́a al caso de variables perfectamente dependientes, X1 = X2 .

La clase completa de funciones de distribución de valores extremos bivariante


se puede obtener mediante una generalización de las distribuciones marginales.
Tomando

  1/ξx1
x1 − µx1
x̃1 = 1 + ξx1 ,
σx1
  1/ξx2
x2 − µx2
x̃2 = 1 + ξx2 ,
σx2
se obtiene

G (x1 , x2 ) = exp {−V (x˜1 , x˜2 )} ,

donde la función V satisface (2.2) para cualquier H.



Las distribuciones marginales serán GEV con parámetros µxj , σxj , ξxj para
j = 1, 2.
La función V será homogénea de orden -1, en el sentido de que para cualquier
constante a > 0 se tiene

1  
1−ω
Z
ω
V a−1 x1 , a−1 x2 =2

máx ,
0 a−1 x1 a−1 x2
Z 1  
ω 1−ω
=2a máx , = aV (x1 , x2 ) .
0 x1 x2
2.3.2. Extensión para n variables
La siguiente caracterización del comportamiento lı́mite del vector de máximos
se va a desarrollar basándose en la pertenencia a un dominio de atracción.
Sea x = (x1 , ..., xp ) ∈ Rd . Si Xi = (Xi,1 , ..., Xi,p ) , i = 1, ..., n, son vectores alea-
torios independientes e idénticamente distribuidos de dimensión p con función
de distribución F , se puede considerar la existencia de secuencias (an )n y (bn )n
pertenecientes a Rp , con an,j > 0 y bn,j ∈ R, ∀ j = 1, ..., p, y una función de
distribución G con marginales no degeneradas tales que, cuando n → ∞,
  
Pr máx Xi − bn /an 6 x = F n (an x + bn ) → G (x) ,
i=1,...,n

donde G es la función de distribución multivariante de valores extremos y F


pertenece al dominio de atracción de G para el máximo.

Ejemplo 1

Se considera la función de distribución normal multivariante FN , con marginales


univariantes N (0, 1) y correlaciones EXi Xj < 1, ∀i, j = 1, ..., p. El dominio de
atracción será el producto de distribuciones marginales Gumbel:

p
Y
n
exp −e−xj ,

FN (an x + bn ) → G (x) =
j=1

con constantes

−1/2
an = (2 log n) ,

bn = bn 1,

donde
1/2 1/2
bn = (2 log n) − 1/2 (log (log n) + log 2π) / (2 log n)
1 = (1, ..., 1) .

2.3.3. Transformaciones
Algunas transformaciones que serán de utilidad son las siguientes:
Transformación de un modelo log-Gumbel en un modelo Fréchet. Una
variable aleatoria Gumbel con parametros de localización y escala log σ y 1/α
se corresponde con una distribución Fréchet con parámetros de forma y escala
α y σ, respectivamente.
Transformación de un modelo Weibull en un modelo Fréchet. Una
variable aleatoria Weibull con parametros de forma y escala −α y 1/σ se corres-
ponde con una distribución Fréchet con parámetros de forma y escala α y σ,
respectivamente.
Trasformación cuantil. Si U es una variable aleatoria distribuida uniforme-
mente en (0, 1), entonces F −1 (U ) tiene función de distribución F .
Transformación de probabilidad. De forma inversa a la transformación
cuantil, si X es una variable aleatoria con función de distribución continua F ,
entonces F (X) tiene función de distribución uniforme (0,1).

2.4. Convergencia del vector de máximos


La convergencia débil de una secuencia de vectores aleatorios implica la conver-
gencia débil de cada uno de los componentes. Por lo tanto, de la misma manera
que en el caso univariante, será razonable aplicar transformaciones a cada uno
de los vectores marginales y considerar la secuencia de variables aleatorias
 
Mn − bn,j
: j = 1, ..., p ,
an,j

para constantes de normalización an,j > 0 y bn,j . Cada uno de los componentes
(Mn − bn,j /an, j), j = 1, ..., p, converge débilmente cuando n → ∞ a una distri-
bución max-estable univariante, si se cumplen las condiciones de convergencia.
Sin embargo, la convergencia de cada una de las p componentes es estrictamen-
te más débil que la convergencia del vector conjunto de máximos. Se necesi-
tará además establecer la estructura de dependencia de la distribución conjunta
F de los vectores aleatorios Xi . Una manera de describir esta dependencia será a
través de cópulas:

P r [Xi 6 x] = F (x) = C1 (F1 (x1 ) , ..., Fp (xp )) .

Suponiendo continuidad de las funciones marginales, la cópula C1 de la función


de distribución F es única y puede obtenerse como la función de distribución
conjunta de los vectores aleatorios (F1 (Xi1 ) , ..., Fp (Xip )).
La cópula del vector de máximos y, por tanto, de cualquier vector creado me-
diante componentes transformadas, vendrá dada por

n  oi
1/i
Ci (u) = C1 u1 , ..., u1/i
p .

Puesto que los vectores de máximos convergen en distribución a una distribución


lı́mite no degenerada, la secuencia de cópulas Ci también debe converger.
La cópulas obtenidas a través de los lı́mites de Ci cuando i → ∞ se denominan
cópulas de valores extremos, esto es, una cópula C es una cópula de valores
extremos si existe una cópula C1 tal que, cuando i → ∞,

n  oi
1/i
lı́m C1 u1 , ..., u1/i
p = C (u1 , ..., up ) .
n→∞
La cópula C1 forma parte del dominio de atracción de C, al pertener esta última
a la clase de posibles cópulas lı́mite.
d
Una cópula C es max-estable si, para todo u ∈ [0, 1] y k = 1, 2, ...,

n  ok
1/k 1/k
C (u) = C u1 , ..., ud .

Por tanto, la clase de cópulas de valores extremos coincide con la clase de cópulas
max-estables.
En resumen, las distribuciones lı́mite no degeneradas de vectores de máximos
transformados tienen marginales de valores extremos y cópulas de valores ex-
tremos o max-estables. En concreto, si
 
p  
\ M n,j − b n,j ω
Pr  6 xj  → G (x1 , ..., xd ) , n → ∞,
j=1
a n,j

entonces,

G (x1 , ..., xd ) = C (G1 (x1 ) , ..., Gp (xp )) .

con marginales de valores extremos G1 , ..., Gp y cópula de valores extremos C.

2.5. Funciones de dependencia


Se partirá del lı́mite que caracteriza a la cópula de valores extremos para, a
continuación, tomar logaritmos y aplicar una expansión lineal. De esta forma,
se llegará a la expresión equivalente

lı́m n 1 − C1 1 − n−1 x1 , ..,1 − i−1 xp


 
n→∞
d
= − log C e−x1 , ..., e−xp = ` (x) ,

x ∈ [0, ∞) .

El lı́mite ` se conoce como función de dependencia de cola estable de C.


Si se considera el vector aleatorio X = (X1 , ..., Xp ),

1 − C1 (1 − x1 /n, ..., 1 − xp /n)


= P r [F1 (X1 ) > 1 − x1 /n ∪ ... ∪ Fp (Xp ) > 1 − xp /n] = ` (x) .

Esta probabilidad hace referencia al evento en el que al menos una de las p


componentes excede un alto percentil de su distribución.
La cópula de cola, R, introducida por Schmidt and Stadtmüler (2006), plan-
tea el caso opuesto donde las p componentes exceden simultáneamente un alto
percentil,
lı́m n P r [F1 (X1 ) > 1 − x1 /n ∩ ... ∩ Fp (Xp ) > 1 − xp /m]
n→∞
= R (x) ,

p
para x ∈ [0, ∞) .
Ambas funciones, ` (x) y R (x), son homogéneas.
Se aplicará la restricción de estas funciones al 1-sı́mplex
p
∆p−1 = {(ω1 , ..., ωp ) ∈ [0, 1] : ω1 + · · · + ωp = 1} .

La restricción ∆p−1 sobre ` se denomina función de dependencia de Pickands.


Por la propiedad de homogeneidad de `,

xj
` (x) = (x1 + · · · xp ) D (ω1 , ..., ωp ) , ωj = .
x1 + · · · + xp

La función ` (x) hace referencia a la probabilidad de la unión de los eventos


{Fj (Xj ) > 1 − xj /n}, siendo la probabilidad de cada uno de estos eventos por
separado xj /n, con xj ≤ n; como consecuencia, las cotas elementales de ` (x)
serán

máx (x1 /n, ..., xp /n) ≤


≤ Pr [F1 (X1 ) > 1 − x1 /n ∪ · · · ∪ Fp (Xp ) > 1 − xp /n] ≤ x1 /n + · · · + xp /n.

Multiplicando por n y haciéndola tender a infinito, se obtiene

p
máx (x1 , ..., xp ) ≤ ` (x1 , ..., xp ) ≤ x1 + · · · + xp , x ∈ [0, ∞) .

Según esta relación anterior, una cópula de valores extremos deberá satisfacer
la relación

u1 · · · up ≤ C (u1 , ..., up ) ≤ máx (u1 , ..., ui ) .

Las cotas inferior y superior de los dos casos anteriores se corresponden, respec-
tivamente, con los casos extremos de independencia y dependecia perfecta.
La cópula C puede ser dada en términos de la función de dependencia de cola
a través de

p
C (u1 , ..., up ) = exp {−` (− log u1 , ..., − log up )} , u ∈ (0, 1] .

En teorı́a de valores extremos, a menudo es conveniente estandarizar a otras


distribuciones distintas a la uniforme (0, 1). La tres formas más comunes son las
distribuciones de Fréchet, de Gumbel y exponencial inversa, que se representan
respectivamente como
 
p
C e−1/x1 , ..., e−1/xp = exp {−` (1/x1 , ..., 1/xp )} , x ∈ (0, ∞) ,
 −x −xp

C e−e , ..., e−e = exp −` e−x1 , ..., e−xp ,
1
x ∈ Rp ,
 

p
(ex1 , ..., exp ) = exp {−` (−x1 , ..., −xp )} , x ∈ (−∞, 0) ,

2.6. Modelización y estimación de parámetros

2.6.1. Caso bivariante

A partir de las series de realizaciones independientes (x1,1 , x2,1 ) , ..., (x1,n , x2,n )
se obtendrá la secuencia de bloques de máximos (z1,1 , z2,1 ) , ..., (z1,m , z2,m ).
Las series de bloques de máximos se considerarán de forma independiente y
se modelizarán como en el caso univariante usando la distribución GEV, es
decir, cada zi,j será tratada como una realización independiente de una variable
aleatoria Zi , para cada i = 1, 2, donde

Zi ∼ GEV (µi , σi , ξi ) .

A continuación se obtendrán las estimaciones de máxima verosimilitud de los


parámetros según la metodologı́a univariante, de forma que la variable transfor-
mada

  ˆi
1/xi
ˆ i Zi − µ̂i
Z̃i = 1 + xi
σ̂i

se distribuye aproximadamente según una distribución Fréchet estándar. Los pa-


res (z̃1,j , z̃2,j ) se obtienen sustituyendo en la expresión anterior las observaciones
(z1,j , z2,j ).
La función de densidad de este modelo es

g (x, y) = {Vx (x, y) Vy (x, y) − Vxy (x, y)} exp {−V (x, y)} , x > 0, y > 0,

donde Vx , Vy y Vx,y son las derivadas parciales y mixtas de V , respectivamente.


Seleccionando un modelo con parámetro θ para V , y a partir de su función de
verosimilitud

m
Y
L (θ) = g (z̃1,i , z̃2,i ) ,
i=1

se obtendrán las estimaciones de máxima verosimilitud de θ.


2.6.2. Extensión para n variables. Estimación mediante
cópulas
Sean xi = (xi,1 , .., xi,d ), i = 1, ...k, realizaciones de una variable que sigue una
distribución de valores extremos G con vectores de parámetros de localización,
escala y forma µ, σ, ξ, respectivamente.
A través de la cópula G se tiene la representación
  
G (x) = C Gµj ,σj ,ξj (xj ) j6d ,

donde Gµj ,σj ,ξj es la j-ésima distribución marginal de G.


Para estimar G se deberán construir estimadores de µ, σ, ξ y C. Las estimaciones
de los parámetros se calcularán de forma independiente para cada función de
distribución marginal Gµj ,σj ,ξj mediante máxima verosimilitud.
La estimación de C se basará en los vectores transformados:
 
zi = Gµ̂j ,σ̂j ,ξ̂j (xi,j ) .
j6d0

2.7. Niveles de retorno

2.7.1. Caso bivariante


Se denominará variable estructural a la operación

Z = φ (Mx1 , Mx2 ) ,

con función de distribución


Z
P r {Z 6 z} = g (x1 , x2 ) dx1 x2 ,
Az

donde Az = {(x1 , x2 ) : φ (x1 , x2 ) 6 z}.


El nivel de retorno en N unidades de tiempo de la variable estructural Z es la
solución de

GZ (z) = 1 − 1/N,

donde GZ es la función de distribución de Z.

2.7.2. Extensión para n variables. Medida de intensidad


Aplicando la teorı́a univariante se tiene que el tiempo de retorno de la ob-
servación Xj es Yj = 1/ {1 − Fj (Xj )}, con distribución de Pareto, ya que
P r [Yj > y] = P r [Fj (Xj ) > 1 − 1/y] = 1/y para y 6 1.
Suponiendo que la cópula C1 se encuentra en el dominio de atracción de una
cópula de valores extremos con función de dependencia `, el vector aleatorio
Y = (Y1 , ..., yp ) satisface

 
p
t {1 − C1 (1 − x1 /t, ..., 1 − xp /t)} = t P r ∪ {Yj > t/xj }
j=1
p
= t P r [Y /t ∈ ([0, ∞] \ [0, 1/x])] → ` (x) , t → ∞.

p
En el espacio Ed = [0, ∞] \ {0}, existe una medida µ, que se conocerá como
medida de intensidad, tal que
" n #
v
X
E I (Yi /n ∈ ·) = n P r [Y /n ∈ ·] → µ (·) , n → ∞,
i=1
R
siendo lı́m n E [f (Y /n)] = Ep f (x) dµ (x) para toda función f continua y aco-
n→∞
tada en Ep y que se hace cero en su origen, de forma que, cuando n crece, el
vector Y /n es desplazado cerca del origen, en las cercanı́as donde la función f
es cero. Por tanto, la medida de intensidad solo afecta a la cola superior de la
distribución de Y .
La medida µ representa el número esperado de observaciones en el suconjunto
sobre el que se aplique esta medida.
La función ` de dependencia actúa como función de distribución de la medida
de intensidad, de forma que

p
` (x) = µ ([0, ∞] \ [0, 1/x]) .

2.8. Modelos de excedencias de umbrales


La aproximación mediante métodos de máximos por bloques presenta varios in-
covenientes; en primer lugar, todos los datos excepto el máximo de cada bloque
son descartados; además, el vector compuesto por estos máximos no se corres-
ponde con observaciones reales.
Una aproximación más flexible y eficiente es considerar las excedencias sobre un
umbral.

2.8.1. Caso bivariante


Sean (x1,1 , x1,2 ) , ..., (xn,1 , xn,2 ) realizaciones independientes de un par de varia-
bles aleatorias (X1 , X2 ) con función de distribución conjunta F .
Para dos umbrales ux1 y ux2 se obtendrán dos funciones de probabilidad de
Pareto con sus correspondientes parámetros para cada una de las distribuciones
marginales de F .
Mediante la siguiente transformación,
(  −1/ξx−1 )!−1
ξx1 (X1 − ux1 )
X̃1 = − log 1 − ζx1 1 + ,
σx1

(  −1/ξx2 )!−1
xix2 (X2 − x2 )
X̃2 = − log 1 − ζx2 1 + ,
σx2
 
se obtiene el par X̃1 , X̃2 con función de distribución conjunta F̃ y marginales
Fréchet para X1 > ux1 y X2 > ux2 . Por tanto,

n o1/n
−1/n
F̃ (x˜1 , x˜2 ) = F̃ n (x˜1 , x˜2 ) ≈ [exp {−V (x˜1 /n, x˜2 /n)}]
= exp {−V (x˜1 , x˜2 )} ,

por la propiedad de homogeneidad de V . Finalmente, dado que F (x1 , x2 ) =


F̃ (x˜1 , x˜2 ), se tiene que

F (x1 , x2 ) ≈ G (x1 , x2 ) = exp {−V (x˜1 , x˜2 )} , x1 > ux1 , x2 > ux2 .

La inferencia en este caso es complicada, ya que dado un par bivariante, puede


ocurrir que solo uno de sus componentes supere el umbral. Para solucionar este
problema, se definen los siguientes espacios:

R0,0 = (−∞, ux1 ) × (−∞, ux2 ) ,

R1,0 = [ux1 , ∞) × (−∞, ux2 ) ,

R0,1 = (−∞, ux1 ) × [ux2 , ∞) ,

R1,1 = [ux1 , ∞) × [ux2 ∞) ,

de forma que, por ejemplo, un punto (x1 , x2 ) ∈ R0,1 si la componente x1 no


supera el umbral ux1 y la componente x2 supera el umbral ux2 .
La inferencia del modelo aproximado podrá llevarse a cabo por máxima verosimi-
litud para puntos pertenecientes a la región R1,1 ; para puntos en otras regiones,
al no poder aplicarse F̃ , será necesario censurar la componente correspondiente
en la función de verosimilitud.
Por ejemplo, si se supone que (x1 , x2 ) ∈ R1,0 , entonces, x1 > ux1 pero x2 < ux2
y los datos solo proporcionarán información de la componente x1 ; de esta forma,
la contribución a la verosimilitud de este punto será

∂F
P r {X1 = x1 , X2 ≤ ux2 } = .
∂x (x1 ,ux )
2
Ası́, para cada una de las regiones se usará la función de verosimilitud:

n
Y
L (θ; (x1,1 , x1,2 ) , ..., (xn,1 , xn,2 )) = ψ (θ; (xi,1 , xi,2 )) ,
i=1

donde

∂ 2 F


 (x1 , x2 ) ∈ R1,1



 ∂x1 ∂x2 (x1 ,x2 )

∂F


(x1 , x2 ) ∈ R1,0



ψ (θ; (x1 , x2 )) = ∂x1 (x1 ,ux ) .
2

∂F


(x1 , x2 ) ∈ R0,1





 ∂x2 (ux ,x2 )
 1


F (ux1 , ux2 ) (x1 , x2 ) ∈ R0,0

Mediante la maximización de la función de log-verosimilitud se obtendrán las


estimaciones y los errores estándar de los parámetros de F . De la misma forma
que en el caso de máximos por bloques, la inferencia puede simplificarse llevando
a cabo estimaciones marginales para después aplicar las transformaciones a las
variables. De este modo, al realizar la maximización, la función de verosimilitud
únicamente dependerá de los parámetros de V .
Se puede utilizar otro método alternativo si existe una variable estructural Z =
φ (X1 , X2 ); en este caso se aplicarán técnicas univariantes de excedencias de
umbrales a cada una de las series zi = φ (xi,1 , xi,2 ).

2.8.2. Extensión a n variables


Sea Xi = (Xi,1 , ..., Xi,p ), con i = 1, ..., n, una muestra de observaciones de
dimensión p, independientes e idénticamente distribuidas con función de distri-
bución conjunta F , perteneciente al dominio de atracción de una distribución
de valores extremos multivariante G.
Se han propuesto varios métodos de análisis de excedencias de umbrales, todos
ellos basados en la aproximación de F en una región extrema definida usando
las condiciones del dominio de atracción especificado en el siguiente teorema.

Teorema 2.8.1 Sea X un vector aleatorio de dimensión p con función de disti-


bución F y marginales Fréchet, y sea {Xi } , i = 1, ..., n, una secuencia extraı́da
de F. Sea también 1 (·) una función indicadora tal que 1 (z) es 1 si z ∈ C y
0 en otro caso. La función de distribución F pertenece al dominio de atrac-
ción de una distribución GEV multivariada con función de distribución G con
marginales Fréchet si y solo si

log F (sx) log G(x)


(a) lı́m =
s→∞ log F (s1) log G(1)

1−F (sx) log G(x)


(b) lı́m =
s→∞ 1−F (s1) log G(1)
 
1 G(x)
(c) lı́m P (X ≤ sx|X > s1) = − log G(1) log G(mı́n(x,1))
s→∞

(d) El proceso puntual Nn , definido como


n
1 (Xi /n ∈ ·) ,
X
Nn (·) =
i=1

d
converge débilmente a un proceso de Poisson no homogéneo en E = [0, ∞] \
{0} con medida de intensidad µ definida como
 
y H (A)
µ y ∈ E : ||y|| > r, = .
||y|| ∈ A r

Smith et al (1997) utilizan la condición (b) para obtener la aproximación


F (x1 , ..., xp ) ≈ 1 + log G (x1 , ..., xp ), válida en una región donde fj (xj ) es cer-
cana a 1 para cada xj , j = 1, ..., p. Sin embargo, este método presenta el in-
coveniente de que cuando F es asintóticamente independiente el sesgo de la
aproximación es significativo.
Ledford y Tawn (1996) corrigen el sesgo de la aproximación anterior usando la
condición (a) para obtener F (x1 , ..., xp ) = G (x1 , ..., xp ), de forma que fj (xj )
es cercana a 1 para cada xj , j = 1, ..., p
En ambos métodos se utiliza un modelo paramétrico para G, obteniéndose a
través de las aproximaciones una parametrización de la cola de F . Los paráme-
tros de los modelos se estimarán por máxima verosimilitud mediante un meca-
nismo de censura. En concreto, si uj es un umbral tal que Fj (uj ) es cercana a
1, todas las observaciones Xi,j ≤ uj estarán censuradas a la derecha por uj ; de
esta forma, la contribución a la verosimilitud de una variable Xi es

∂q
li = F (x1 , ..., xp )
∂xj1 · · · ∂xjq

evaluada en (máx (Xi,1 , u1 ) , ..., máx (Xi,p , up )), donde jl es el ı́ndice de las va-
riables marginales para las que Xi,jl > ujl y q es el número de variables que
cumplen la condición. Normalmente, se consideran las Fj marginales Fréchet,
siendo los uj umbrales localizados en altos cuantiles de esta distribución.
Coles y Tawn (1991) proponen otro método de análisis de excendencias de um-
brales partiendo de la condición (d). Para n grande, el proceso puntual Nn de
esta condición se comporta como un proceso de Poisson no homogéneo en E.
Si existen un número nA ≤ n de observaciones, X(1) , ..., X(nA ) , pertenecien-
c
tes a A = (0, u] , donde u es un vector de umbrales, la verosimilitud para la
aproximación Poisson es

n nA
 o Y  
L = exp −µ Ã µ dX̃(k) ,
k=1

donde à = A/n y X̃(k) = X(k) /n. Los parámetros se estimarán por máxima
verosimilitud.
Se puede desarrollar otro método alternativo de excedencias de umbrales a partir
de la condición (c) aproximando la distribución condicional de Xi dado que
Xi ∈ A mediante
 
1 G (x, α)
W (x) = log , x > 0,
− log G (1, α) G (mı́n (x, 1) , α)
c
donde A = (0, u1] y u es un vector de umbrales.
Mediante una parametrización de G se obtiene un modelo parámetrico para W
cuya estimación se realizara por máxima verosimilitud.
Hasta ahora se ha considerado que la distribución F tiene marginales Fréchet;
sin embargo, en la práctica las distribuciones marginales Fj son desconocidas y
deberán ser estimadas. Una solución habitual a este problema es la realización
de una aproximación semiparamétrica, donde Fj (xj ) se estima a través de la
distribución empı́rica para xj ≤ uj , y una modelización mediante la distribución
de Pareto generalizada para xj > uj :
   −1/ξj
 xj − uj
 1 − λ j 1 + ξj , si xj > uj


σj

+
Fj (xj ) = n
−1
1 (Xi,j ≤ xj ) ,
 X
(n + 1) si xj ≤ uj ,




i=1

donde λj > 0, σj > 0 y j ∈ R son parámetros desconocidos que deberán esti-


marse. Mediante las transformaciones −1/ log Fj (Xi,j ) , i = 1, ..., n, j = 1, ..., p
se obtienen las observaciones estandarizadas a escala Fréchet y se pueden aplicar
los procedimientos vistos en esta sección.
Capı́tulo 3

Extremos espaciales

3.1. Introducción

En algunos casos, principalmente en el estudio de fenómenos medioambientales,


se tendrá una dimensión espacial y el objetivo será modelizar la dependencia
espacial entre extremos a partir de observaciones localizadas en una malla. Para
ello, se usarán procesos max-estables, es decir, procesos estocásticos en los que
las distribuciones finito-dimensionales serán distribuciones de valores extremos.
Los extremos espaciales son de interés en el estudio de riesgos naturales como
olas de calor, altas precipitaciones y nevadas, inundaciones, etc.
Como ejemplo de aplicaciones del estudio de extremos espaciales, en el contexto
de las precipitaciones se tiene a Davidson, Padoan y Ribatet (2012), que reali-
zaron una comparación entre modelos de variables latentes y modelos basados
en procesos max-estables para un conjunto de precipitaciones estivales máximas
recogidas entre 1962-2008 en 52 localizaciones en la región de Plateau (Suiza).
En cuanto a inundaciones costeras, Coles y Tawn (1990), consideraron los niveles
máximos de la marea en la costa Británica, construyendo un modelo logı́stico
bivariante para localizaciones cercanas.

3.2. Procesos max-estables espaciales

Para el caso de extremos espaciales son fundamentales los procesos max-estables,


análogos infinito-dimensionales a los vectores aleatorios max-estables, definidos
en la sección anterior.
En este caso, los extremos estarán definidos en el dominio espacial R ⊂ Rd ,
d ≥ 1, y se supondrá en el caso teórico que R es un subconjunto compacto
de Rd y que todos los procesos estocásticos considerados siguen patrones de
muestreo continuo, aunque estas últimas condiciones se podrán relajar en la
práctica.
Sea Z (r) un proceso estocastico con marginales no degeneradas definido en

37
un conjunto R. Este proceso sera max-estable si todas sus distribuciones n-
dimensionales satisfacen la propiedad de max-estabilidad, es decir:

Definición 3.2.1 Un proceso estocástico {Z (r) : r ∈ R} se dice max-estable si


existen sucesiones de funciones continuas {an (r) > 0 : r ∈ R, n ≥ 1} y
{bn (r) ∈ R : r ∈ R, n ≥ 1} tales que, para todo n ≥ 1,
 
 máx Zi (r) − bn (r) 
i=1,...,n d
:r∈R = {Z (r) : r ∈ R} ,
 an (r) 

donde {Zi (r) : r ∈ R, i ≥ 1} es una sucesión de réplicas independientes de {Z (r) : r ∈ R}.

La relevancia de los procesos max-estables en la modelización de extremos es-


paciales estará basada en argumentos asintóticos, según según de Haan (1984):

Teorema 3.2.1 Sea X1 (r) , X2 (r) , ... una sucesión de réplicas independientes
de un proceso estocástico {Xn (r) : r ∈ R}. Si existen funciones continuas cn > 0
y dn ∈ R tales que el proceso lı́mite {Zn (r) : r ∈ R}, definido como

máx Xi (r) − dn (r)


i=1,...,n
→ Zn (r) , r ∈ R, n → ∞,
cn (r)

es no degenerado, entonces, el proceso {Zn (r) : r ∈ R} es un proceso max-


estable.

Una descripción más precisa de los procesos max-estables se puede obtener a


partir de su representación espectral.

3.2.1. Representación espectral de los procesos max-estables


El proceso aleatorio {Z (r) : r ∈ R} será un proceso max-estable simple si tiene
marginales Fréchet.

Teorema 3.2.2 (de Haan, 1984; Penrose, 1992) Cualquier proceso max-estable
simple no degenerado {Z (r) : r ∈ R} definido en un conjunto compacto R ⊂ Rd ,
d ≥ 1, con patrones de muestreo continuos satisface

d
Z (r) = máxζi fi (r) , r ∈ R,
i≥1

donde {(ζi , fi ) : i ≥ 1} son las coordenadas de un proceso de Poisson en (0, ∞)×


C con medida de intensidad ζ −2 dζν (df ) para alguna medida local finita ν defi-
nida en el espacio de las funciones continuas no negativas C en R, tal que
Z
f (r) ν (df ) = 1, r ∈ R.
3.2.2. Modelos max-estables
En esta sección se van a introducir algunos modelos max-estables que serán
relevantes en el estudio de extremos espaciales.

Modelo de Smith (1990)

Smith (1990) introdujo la siguiente representación paramétrica de un proceso


max-estable continuo en el tiempo:

Z (r) = máxXi f (r, Ui ) ,


i≥1

donde {Xi , Ui }i≥1 es una realización de un proceso de Poisson definido en


(0, ∞) × RU con medida de intensidad x−2 dxν (du) y f (·) es una función de
densidad no negativa en R × RU , con integral finita respecto a una medida
positiva ν en RU .
El
R proceso Z (·) será un proceso max-estable simple si imponemos la restricción
RU
f (r, u) ν (du) = 1.
En este proceso, que Smith (1990) interpretó en términos de precipitaciones, Xi
y Ui son, respectivamente, el tamaño y el tipo de tormenta i.
Smith (1990) definió el proceso de valores extremos gaussianos como

  
1 a (h) 1 z2
P (Z (r) ≤ z1 , Z (r + h) ≤ z2 ) = exp − Φ + log
z1 2 a (h) z1
 
1 a (h) 1 z1
− Φ + log , h ∈ Rd ,
z2 2 a (h) z2

donde Φ (·) es la función de distribución normal estándar y a (h) = h0 Σ−1 h.
Esta representación se usará en modelos medioambientales.

Modelo de Schlather (2002)

Schlather (2002) propuso una extensión del modelo de Smith (1990), sustitu-
yendo la función determinı́stica f (·) por una forma aleatoria Y (·):

Z (r) = máxXi Yi (r − Ui ) ,
i≥1

donde {Xi , Ui }i≥1 es una realización de un proceso de Poisson en (0, ∞) × Rd


con medida de intensidad x−2 dx × µ−1 du e {Yi (·)}i≥1 son réplicas indepen-
dientes
R de una función aleatoria no negativa Y (·) definida en Rd , tal que µ =
E R Y (r) dr ∈ (0, ∞).
Schlather (2002) también introdujo una clase de procesos max-estables basada
en procesos estocásticos estacionarios con media finita.
Sea {Xi }i≥1 una realización de un proceso de Poisson definido en (0, ∞) con
medida de intensidad µ−1 x−2 dx y sean {Wi (·)}i≥1 replicas independientes de un
proceso estocástico estacionario W (·) definido en R2 , con µ = E (máx {0, W (0)}) ∈
(0, ∞).
El proceso estocástico Z (·) definido por

Z (r) = máxXi Wi (r)


i≥1

es un proceso max-estable simple estacionario, con distribución multivariante,


para cualquier subconjunto {r1 , ..., rm } de R:
  
W (ri )
P (Z (r1 ) ≤ z1 , ..., Z (rm ) ≤ zm ) = exp −E sup .
1≤i≤m zi

Schlather (2002) propuso tomar Yi (·) como un proceso estacionario gaussiano


con función de correlación ρ. El proceso max-estable resultante se conoce como
proceso de extremos gaussiano y su distribución bivariante viene dada por

  
1 1 1
P (Z (r) ≤ z1 , Z (r + h) ≤ z2 ) = exp − +
2 z1 z2
!)
z1 z2
r
× 1 + 1 − 2 (ρ (h) + 1) 2 .
(z1 + z2 )

3.3. Medidas de dependencia para extremos es-


paciales

3.3.1. Función del coeficiente extremal


El coeficiente extremal se utilizará para medir la dependencia entre los extre-
mos de dos vectores aleatorios (X, Y ) con distribución de extremos bivariante y
distribución marginal común F (r). Este coeficiente se define por la identidad

P (máx (X, Y ) ≤ x) = F θ (x) ,

donde el rango de θ será [1, 2], siendo θ = 1 dependencia perfecta y θ = 2


independencia entre los máximos.
Para un proceso max-estable Z (r), la función del coeficiente extremal, θ (·),
viene dada por (Schlater y Tawn, 2003)

P (máx (Z (r) , Z (r + h)) ≤ z) = exp {−θ (h) /z} .

Algunas propiedades de θ (·) son:


1. 2 − θ (·) es una función semi-definida positiva.

2. La función θ (·) no es diferenciable en 0, excepto en el caso θ (h) = 1, ∀h

3. Si Z (·) es un proceso estocástico isotropico, θ (·) solo tiene un punto de


discontinuidad en 0.

3.3.2. F-madograma

En geoestadı́stica, la herramienta básica para medir la dependencia espacial es


el semi-variograma, que se define como

1 h 2
i
γ (h) = E {W (r) − W (r + h)} , r, r + h ∈ R.
2

La herramienta análoga al semi-variograma para la teorı́a de extremos fue defi-


nida por Cooley et al. (2006) y se conoce como F -madograma:

1
νF (h) = E [|F {Z (r + h)} − F {Z (r)} |] , r, r + h ∈ R1
2

donde F es la función de distribución acumulativa de Z (r).


El rango de valores del F -madrograma varı́a de 0 a 1/6, lo que se corresponde
con dependencia completa e idependencia respectivamente.
En la práctica se estimará de la siguiente forma:

n n
1 X X
ν̂F (h) = |Ri (r) − Ri (r + h) |, Ri (r) = 1{Z` (r)≤Zi (r)} ,
2n (n + 1) i=1
`=1

donde Z1 , ..., Zn son réplicas independientes de un proceso max-estable Z.

3.3.3. λ-madograma

El coeficiente extremal, ası́ como el F -madograma, no representan de forma


completa la dependencia espacial de un campo aleatorio.
Para resolver este problema, Naveau et al. (2009) introdujo el λ-madograma,
que se define como

1  λ
E |F {Z (r1 )} − F 1−λ {Z (r2 )} | ,

νλ (r1 , r2 ) = ∀ λ ∈ [0, 1] .
2

De esta forma, mediante la variación de λ se tendrá P [Z (r1 ) ≤ z1 , Z (r2 ) ≤ z2 ],


donde z1 = λz y z2 = (1 − λ) z, y de esta forma se podrá explorar todo el
espacio.
3.4. Simulación de procesos espaciales max-estables
Las simulaciones espaciales permiten la recuperación de información en cual-
quier localización.
Se llamarán simulaciones condicionadas aquellas que se ajustan a los datos ob-
servados; en otro caso, se denominarán no condicionadas. La simulación de las
primeras será especialmente útil ya que se podrán combinar diferentes simula-
ciones condicionadas para obtener estimadores de cualquier medida de interés,
como cuatiles, probabilidades de excedencia de umbrales o niveles de retorno.

3.4.1. Simulaciones no condicionadas


Los procesos max-estables se han definido basándose en infinitas réplicas de un
proceso estocástico; sin embargo, en la práctica solo se podrá obtener un número
finito de realizaciones de un proceso. A pesar de esto existen resultados teóricos
y procedimientos ad-hoc mediante los que será posible obtener simulaciones
exactas o aproximadas de algunos modelos max-estables.
En Schlather (2002) se dan algunas condiciones para la obtención de simulacio-
nes exactas para la clase de procesos max-estables

Z (r) = máxXi Wi (r) .


i≥1

Suponemos que el proceso estocástico W (·) está uniformemente acotado por


una constante positiva C < ∞. Para T0 = 0 y k = 1, 2, ..., el algoritmo de
simulación será el siguiente:

1. Se genera Ek ∼ E y se hace Tk = Tk−1 + Ek , Xk = Tk−1 ,

2. se genera Wk (·) ∼ W (·),

3. si CXk > máx Xi Wi (r), volver a 1; en otro caso Z (r) = máx Xi Wi (r),
1≤i≤k 1≤i≤k

donde E es una distribución exponencial con media unitaria, (Ti )i ≥ 1 es un


proceso de Poisson definido en (0, ∞) con medida de intensidad dt y (Xi )i ≥ 1
es un proceso de Poisson definido en (0, ∞) con medida de intensidad x−2 dx.
Se obtendrá una buena aproximación al elegir un C tal que P (W (r) > C) sea
lo suficientemente pequeña.
Para el proceso Z (r) = máxXi Yi (r − Ui ) se puede aplicar un algoritmo análogo,
i≥1
ya que Y (·) es un proceso uniformemente acotado.

3.4.2. Simulaciones condicionadas


La simulación condicionada, a diferencia de la no condicionada, no puede obte-
nerse de forma exacta.
Wang y Stoev (2011) dan una solución aproximada a este problema. Suponga-
mos que se observa un proceso max-estable Z (·) en las ubicaciones r1 , ..., rn .
0
La distribución de un vector max-estable (Z (r1 ) , ..., Z (rn )) se podrá apro-
 0
ximar mediante el vector multivariante Z̃ (r1 ) , ..., Z̃ (rn ) , donde Z̃ (ri ) =
máx φk (ri ) Yk , para i = 1, ..., n, siendo φk (·) funciones determinı́sticas no ne-
k=1,...,p
gativas e Yk variables aleatorias independientes con distribución α-Fréchet, es
decir, P (Yk ≤ y) = exp (−σkα y −α ), para α , σk , y > 0.
El algoritmo de Wang y Stoev (2011) genera muestras a partir de la probabilidad
 0
0 0
condicionada de (Y1 , ..., Yp ) dadas Ẑ (r1 ) , ..., Ẑ (rn ) = (Z (r1 ) , ..., Z (rn )) y
predice Z (r∗) en cualquier ubicación arbitraria r∗ a partir de las muestras
 0
(sim) (sim)
Yi , ..., Yp , usando la siguiente aproximación:

(sim)
Ẑ (sim) (r∗) = máx φk (r∗) Yk .
k=1,...,p

3.5. Ajuste de un proceso max-estable Fréchet


a los datos

En esta sección se van a presentar dos aproximaciones diferentes para realizar


el ajuste de procesos max-estables a los datos. La primera estará basada en
mı́nimos cuadrados y para la segunda se usará el estimador compuesto máximo-
verosı́mil propuesto por Lindsay (1988).

3.5.1. Mı́nimos cuadrados

El procedimiento consiste en minimizar la función

 2
X θ −
 i,j  θ̃ i,j 
C (ψ) = ,
i<j s θ̃i,j

donde ψ es el vector de parámetros del proceso max-estable, θi,j es el coeficien-


te extremal predicho para el modelo y θ̃i,j es un estimador semi-paramétrico
definido como

n
θ̃i,j = P n o
n −1 −1
k=1 mı́n Z k (ri ) , Zk (rj )

 
con desviación estándar s θ̃i,j .
3.5.2. Máxima verosimilitud
Este método presenta el incoveniente de que la verosimilitud del modelo no
se conoce de forma analı́tica para dimensiones mayores o iguales que tres. Sin
embargo, dado que la densidad bivariante sı́ es conocida, Ribatet (2009) sugiere
el uso funciones de verosimilitud dos a dos. La expresión de la log-verosimilitud
dos a dos viene dada por

ni,j  
(i) (j)
XX
`p (z; ψ) = log f zk , zk ; ψ ,
i<j k=1

donde z es el conjunto total de observaciones de la región, ni,j es el número


(i)
de observaciones entre las localizaciones i y j, yk es la k-ésima observación en
la localización i y f (·, ·) es la densidad bivariante de un proceso max-estable
Fréchet.
Esta estimación por máxima verosimilitud estará especificada parcialmente;
es
 decir, la variable
de la que se han extraı́do las observaciones del modelo,
f (y; ψ) , ψ ∈ Rd , tiene una desidad g desconocida. El modelo quedará espe-
cificado si existe ψ∗ ∈ Rd tal que f (y; ψ∗ ) = g (y) para todo y.
La estimación del vector de parámetros se distribuirá como
 
−1 −1
ψ̂ = N ψ, H (ψ) J (ψ) H (ψ) ,

donde

∂ 2 log f (y; ψ)
Z
H (ψ) = n g (y) dy
∂ψ∂ψ T
Z
∂ log f (y; ψ) ∂ log f (y; ψ)
J (ψ) = n .
∂ψ ∂ψ T

3.6. Selección del modelo


Dados varios modelos a los que se ajustan bien a los datos, será útil disponer
de herramientas que nos permitan elegir el más adecuado. Si los modelos a
comparar tienen la misma log-verosimilud, se seleccionará el modelo más simple;
sin embargo, si difieren en una mı́nima cantidad se deberá disponer de algún
criterio para discriminar un modelo frente a otro.
Se van a presentar dos criterios de selección de modelos; el Criterio de Infor-
mación de Takeuchi (Akaike, 1974) y el Estadı́stico de la Tasa de Verosimilud
(Davison, 2003).

3.6.1. Criterio de información de Takeuchi


Este criterio se utiliza para comparar modelos dos a dos.
Se considera una muestra aleatoria Y1 , ..., Yn extraı́da de una variable con fun-
ción de densidad g desconocida.
Se ajustará un modelo estadı́stico f (y; ψ) maximizando la log-verosimilitud. La
medida de divergencia de Kullback-Leibler (Varin, Vidoni, 2005) determina la
discrepancia entre el modelo ajustado f y el modelo real g,
Z  
g (y)
D (fψ , g) = log g (y) dy
f (y; ψ)

Esta medida siempre será positiva; y se tomará como modelo aquel que minimice
D (fψ , g).
Dado que algunos modelos satisfacen D (fψ , g) = 0, este método no es lo sufi-
cientemente discriminante; para solucionar este problema se utiliza el Criterio
de Informacion de Takeuchi (TIC),
   
T IC = −2` ψ̂ − 2tr JˆĤ −1 ,

donde ψ̂ es la estimación de ψ y Ĥ y Jˆ son estimaciones consistentes de las


matrices H (ψ) y J (ψ), respectivamente.
El mejor modelo se corresponderá con aquel que minimice el TIC.

3.6.2. Estadı́stico de la tasa de verosimilitud


Este criterio es útil para comparar modelos anidados.

Para un modelo estadı́stico {f (x, ψ)}, con ψ T = κT , φT , se quiere comprobar
que los datos son consistentes con la hipótesis de que κ = κ0 ; para ello, se
usará el estadı́stico de la tasa dde verosimilitud W (κ0 ),

n    o p
X
W (κ0 ) = 2 ` κ̂, φ̂ − ` κ0 , φ̂κ0 → λi Xi , n → ∞,
i=1

donde φ̂κ0 es la estimación por máxima verosimilitud de φ bajo la restricción


κ = κ0 , p es la dimensión del parámetro κ0 , Xi son variables independientes
con distribución χ21 y λi son los autovalores de

−1
H −1 JH −1 H −1
 
κ κ
.
Capı́tulo 4

Aplicación con R. Extremos


espaciales

Se va a trabajar con el paquete de R SpatialExtremes y con el conjunto de da-


tos wind que recoje los máximos anuales de la velocidad del viento, registrados
en 35 estaciones meteorológicas de los Paises Bajos en el periodo 1971-2012.

Figura 4.1: Localizaciones de las estaciones meteorológicas

47
maps::map(xlim=c(0,9),ylim=c(50,54),fill=TRUE,col="green3")
points(coord,pch=15)

Para el tratamiento de los datos se seguirá el siguiente proceso propuesto por


Smith (1990):

1. Ajuste de la distribución de valores extremos generalizada a los datos.

2. Tranformación de los datos a la distribución Fréchet.

3. Ajuste de un proceso max-estable a los datos transformados.

4. Diagnosis y elección del modelo.

5. Predicciones.

4.1. Ajuste de la distribución de valores extre-


mos generalizada a los datos

Para el ajuste de la distribución GEV será necesario comprobar si los datos


presentan dependencia espacial y localizar tendencias espaciales si las hubiera.
La forma de comprobar la existencia de dependencia espacial es a través de la
representación del F -madograma y el coeficiente extremal.
● ●

● ●
1.7

1.7
● ●
● ●
● ●
● ●
● ●



● ●● ●
● ●

● ● ●
●● ●

● ● ● ●
1.6

● ● ●
● ● ●

1.6
● ●● ● ●● ●
● ●

Ext.Coeff
Ext.Coeff

● ● ●
●● ● ● ● ●
● ● ●● ● ●● ● ● ●

●● ● ●
● ● ● ●
● ● ● ●
● ● ● ● ●●● ● ●
●● ● ● ● ● ● ●
● ● ●● ● ●
●●● ● ● ●● ●

●●
●● ●●
● ● ● ● ●●
● ●

● ● ● ●
● ● ●● ● ● ●● ● ●
● ● ● ●● ● ● ● ●●
● ●
1.5


● ●● ● ●● ●
●●
●● ● ●●● ●● ●
● ● ●● ● ● ● ●● ●● ● ●● ●
●●

1.5
●● ● ● ●

●●●● ●●

● ●
● ● ●●● ●


●●

● ●● ● ● ●
●● ● ●
●● ●
● ● ● ● ●● ● ●● ●
● ●● ● ●● ●
● ●● ●●
● ●
● ● ● ●
● ●
●●● ●●●● ● ●

● ● ● ●
●● ●
● ●●● ● ●
● ● ●
● ● ●
1.4

● ●
● ●●●●

1.4

● ●

0 20 40 60 80 100 120 0 20 40 60 80 100 120

Distancia Distancia

(a) Longitud vs. Altitud (b) Latitud vs. Altitud


1.7

● ● ●
● ● ●
● ●●

● ● ● ●
●● ● ●●
1.6

● ● ●●

● ● ● ●
●●
Ext.Coeff

● ● ●
● ●
● ● ●
● ● ●
● ● ● ●● ● ● ● ● ● ●

● ●
●● ●●
● ● ●● ● ● ●● ● ●
● ● ●
●● ● ● ●●
● ● ● ●

● ● ● ●
● ●● ●
1.5

●● ● ● ● ●

● ●● ● ●●
● ●
● ● ● ●
● ● ●●
● ●
●● ●
● ● ● ● ●
● ● ● ●
● ● ●
●● ●●
● ● ● ●
● ●
● ●
1.4


● ●

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5

Distancia

(c) Longitud vs. Latitud

Figura 4.2: Representación del coeficiente extremal frente a la distancia

En estos gráficos no se observa dependencia espacial en función de la altitud,


únicamente para las dimensiones longitud y latitud se observa una variación del
coeficiente extremal en función de la distancia.

Los siguientes gráficos ayudarán a detectar tendencias espaciales:


● ● ● ● ● ● ● ● ● ● ● ● ● ●
± 10 ± 20 ± 30 ± 40 ± 50 ± 10 ± 20 ± 30 ± 40 ±1 ±2 ±3 ±4 ±5

Figura 4.3: Gráficos de tendencias para µ (x), σ (x) y xi (x), respectivamente.

symbolplot(wind,coord,which="mean")
symbolplot(wind,coord)
border<-function(add=FALSE)maps::map(xlim=c(0,9),ylim=c(47.5,57.5),add=add)
symbolplot(wind,coord,plot.border=border,scale=2)
En el caso de los parámetros σ y ξ no se observa ningún patrón. En el gráfico de
la izquierda, correspondiente al parámetro de escala µ, se observa una separación
de los datos, por lo que este parámetro se definirá en función de las coordenadas.
La definición de los parámetros propuesta por Ribatet (2015) es la siguiente:

µ (x) = β0,µ + β1,µ lon (x) + β2,µ lat (x) + β3,µ lon (x) lat (x)
σ (x) = β0,σ
ξ (x) = β0,ξ

loc.M1<-y lon*lat; scale.M1<-shape.M1<-y 1

El ajuste de la distribución GEV a los datos da como resultado el valor de los


parámetros estimados junto con sus errores estándar y la matriz de varianzas
covarianzas.
# Ajuste de la distribución GEV a los datos
M1<-fitspatgev(wind,scale(coord,scale=FALSE),loc.M1,scale.M1,shape.M1);
M1
Model: Spatial GEV model
Deviance: 10735.05
TIC: 10768.3

Location Parameters:
locCoeff1 locCoeff2 locCoeff3 locCoeff4
267.56 -15.96 22.85 2.23
Scale Parameters:
scaleCoeff1
33.86
Shape Parameters:
shapeCoeff1
-0.09237
Standard Errors
locCoeff1 locCoeff2 locCoeff3 locCoeff4 scaleCoeff1 shapeCoeff1
3.12232 1.04371 1.71802 1.58605 1.03977 0.01854

Asymptotic Variance Covariance


locCoeff1 locCoeff2 locCoeff3 locCoeff4 scaleCoeff1
locCoeff1 9.7488927 -0.5203683 1.5094327 -0.4169546 1.6639010
locCoeff2 -0.5203683 1.0893215 -0.7382272 0.2126821 -0.0542326
locCoeff3 1.5094327 -0.7382272 2.9515839 0.9878967 0.3578833
locCoeff4 -0.4169546 0.2126821 0.9878967 2.5155629 -0.2183085
scaleCoeff1 1.6639010 -0.0542326 0.3578833 -0.2183085 1.0811170
shapeCoeff1 -0.0319854 0.0102704 -0.0072285 0.0031509 0.0018621
shapeCoeff1
locCoeff1 -0.0319854
locCoeff2 0.0102704
locCoeff3 -0.0072285
locCoeff4 0.0031509
scaleCoeff1 0.0018621
shapeCoeff1 0.0003439

Optimization Information
Convergence: successful
Function Evaluations: 529

El ajuste de este modelo nos permitirá predecir los niveles de retorno para un
determinado número de años.

x<-seq(min(coord[,1]),max(coord[,1]),length=100)
y<-seq(min(coord[,2]),max(coord[,2]),length=100)
grid<-expand.grid(x,y);
colnames(grid)<-c("lon","lat")
grid[,1]<-grid[,1]-mean(coord[,1]);
grid[,2]<-grid[,2]-mean(coord[,2]);
ans<-predict(M1,newdata=grid,ret.per=100)$Q100
maps::map(xlim=range(x),ylim=range(y))
image(x,y,matrix(ans,100),add=TRUE,col=cm.colors(64))
contour(x,y,matrix(ans,100),add=TRUE);
maps::map(add=TRUE)

0 0 0
42 41 40
0
44

0
43
0
39

0
38

0
37

0
36
0
35
0
34

Figura 4.4: Predición de los niveles de retorno a 100 años

4.2. Transformación de los datos a la distribu-


ción Fréchet.
En primer lugar se transformarán los datos partiendo de la distribución GEV
calculada en la sección anterior; para ello, se tendrán en cuenta los valores de
los parámetros ajustados:

loc.coeff.M1<-M1$fitted.values[1:4];
sigma.M1<-M1$fitted.values[5];
xi.M1<-M1$fitted.values[6]
mu.M1<-numeric(nrow(coord))
for(i in 1:nrow(coord))
mu.M1[i]<-loc.coeff.M1[1]+loc.coeff.M1[2]*coord[i,1]+loc.coeff.M1[3]*coord[i,2]+
loc.coeff.M1[4]*coord[i,1]*coord[i,2]
>mu.M1 [1] 1903.066 1928.571 1965.237 1943.246 1992.469 1975.527 2038.791
1976.932 [9] 1987.656 2025.385 2024.461 2074.960 2069.942 2046.741
2127.005 2098.129 [17] 2141.369 2156.563 2123.485 2217.123 2155.867
1798.158 1817.627 1829.274 [25] 1867.950 1872.704 1898.695 1947.514
1934.438 1965.384 1977.807 2014.147 [33] 2003.479 1993.622 2056.935
>sigma.M1
scaleCoeff1
33.85732
>xi.M1
shapeCoeff1
-0.09236842

Con el siguiente código y el uso de la función gev2frec() se transforma el con-


junto de datos original en otro conjunto de datos (wind.frechet.M1) del mismo
tamaño, en el que cada uno de los vectores tiene marginales Fréchet:

wind.frechet.M1<-numeric(length(wind));
wind.frechet.M1<-matrix(wind.frechet.M1,nrow=nrow(wind),ncol=ncol(wind));
for(j in 1:ncol(wind))
for(i in 1:nrow(wind))
wind.frechet.M1[i,j]<-gev2frech(wind[i,j],mu.M1[j],sigma.M1,xi.M1)

Además se realizara una segunda transformación en la que los datos serán trans-
formados usando la distribución acumulada empı́rica, obteniendose el conjunto
de datos wind.frechet.M2.

wind.frechet.M2<-apply(wind,2,gev2frech,emp=TRUE)

4.3. Ajuste del proceso max-estable a los datos


transformados.

4.3.1. Modelo 1

Se ajustará sobre el conjunto de datos wind.frechet.M1.


De la misma forma que en el ajuste de la distribución GEV, habrá que definir
la forma de los parámetros del modelo:
symbolplot(log(wind.frechet.M1),coord[,-3],which="mean")
symbolplot(log(wind.frechet.M1),coord[,-3])
border<-function(add=FALSE)maps::map(xlim=c(0,9),ylim=c(47.5,57.5),add=add)
symbolplot(log(wind.frechet.M1),coord[,-3],plot.border=border,scale=10)

± 0.5

±1

± 1.5
● ● ● ●
± 0.05 ± 0.1 ± 0.15 ± 0.2
● ● ●●●●
±1 ±2 ±3 ±4 ±5 ±6

Figura 4.5: Gráficos de tendencias para µ (x), σ (x) y ξ (x), respectivamente.

En este caso, el parámetro de localización va a depender de las coordenadas,


mientras que los parámetros de escala y forma serán constantes ya que no se
aprecia una tendencia clara en sus respectivos gráficos.

µ (x) = β0,µ + β1,µ lon (x) + β2,µ lat (x) + β3,µ lon (x) lat (x) ,
σ (x) = β0,σ
ξ (x) = β0,ξ

Equivalentemente en R.

loc.M1<-y ∼ lon*lat; scale.M1<-shape.M1<-y ∼ 1;

El ajuste del modelo se hará con la función fitmaxstab().

M1.maxstab<-fitmaxstab(log(wind.frechet.M1),coord[,-3],"whitmat",nugget=0,
loc.M1,scale.M1,shape.M1)

4.3.2. Modelo 2
Se repite el mismo proceso para el conjunto de datos wind.frechet.M2.
symbolplot(wind.frechet.M2,coord[,-3],which="mean")
symbolplot(wind.frechet.M2,coord[,-3])
border<-function(add=FALSE)maps::map(xlim=c(0,9),ylim=c(47.5,57.5),add=add)
symbolplot(wind.frechet.M2,coord[,-3],plot.border=border,scale=3)

● ● ● ● ● ● ● ● ● ● ● ● ● ●
± 0.05
± 0.1
± 0.15
± 0.2
± 0.25
± 0.3 ± 0.1 ± 0.2 ± 0.3 ± 0.4 ± 0.2 ± 0.4 ± 0.6 ± 0.8

Figura 4.6: Gráficos de tendencias para µ (x), σ (x) y ξ (x), respectivamente.

No se aprecia ningún patrón en la distribución espacial de los puntos, todos los


parámetros serán constantes:

loc.M2<-scale.M2<-shape.M2<-y 1;
M2.maxstab<-fitmaxstab(wind.frechet.M2,coord[,-3],"whitmat",nugget=0,loc.M2,scale.M2,shape.M2)

4.3.3. Diagnosis y elección del modelo


Para verificar la validez de cada uno de los modelos se va a analizar que las
observaciones estén bien modelizadas en cada una de las localizaciones y que la
estructura de dependencia esté bien definida; esto se hará mediante dos funcio-
nes: qqgev() en el caso del Modelo 1 y fmadogram() para el Modelo 2.

Modelo 1

qqgev(M1.maxstab)
σModel ξModel

−17.5

15

−18.0

15

●●


−18.5



●●

10
Frequency

Frequency

10

µModel

●●



−19.0

●●

●● ●




5

5

−19.5

●●


−20.0


0

−20.0 −19.0 −18.0 0.00 0.10 0.20 0.30 −1 0 1 2 3 4 5 6

µMLE σMLE ξMLE

Figura 4.7: Diagnosis de los parámetros del modelo.

Con los gráficos qq se podrá comprobar si los parámetros predichos a través


de las tendencias de superficie son relevantes; para ello, se compararán con los
estimados por máxima verosimilitud en cada una de las localizaciones:
fmadogram(wind.frechet.M1,coord[,-3],M1.maxstab,which=c(.ext"),col=c("grey"))
fmadogram(wind.frechet.M1,coord[,-3],M1.maxstab,n.bins=200,which=c(.ext"),col=c(1,2),m

2.0

● ●


● ●

● ●


● ●●
1.8

● ● ●●
● ● ● ●
● ●● ● ●●●● ● ●
●● ● ●● ● ●● ● ●
● ● ●
● ●
● ●● ● ●● ● ●
● ●●●●● ● ●● ● ● ● ● ●
● ● ● ●●● ●
● ●● ● ●● ● ● ● ● ●●●●● ● ● ● ● ●
● ● ●● ● ●
●●● ● ●● ● ●
●● ● ● ● ● ●●
● ● ● ●● ● ●● ● ●●

● ●● ●
●●
● ●● ●
●●●● ●● ●●
● ●● ● ● ●●● ● ● ● ● ●● ●● ● ● ●
1.6

● ● ● ● ● ● ●● ● ●
● ●●● ●
● ●● ● ● ●● ●● ●● ●● ●● ● ●●
● ●●●●●● ● ● ● ● ●●●● ●
●● ● ●●●● ● ●

●●●●●● ●
●●
●● ●
●●● ●● ●●
●● ● ● ●●
● ●● ● ● ● ●●
● ● ● ●● ● ●● ●●● ●●●● ● ●● ● ● ●● ●●●●● ●● ●

● ●●● ● ●●
●● ●
●●● ● ●
●●●● ●●●
● ●
●● ● ● ●● ● ●● ● ●
● ●● ●●● ● ● ●● ● ●● ●● ● ●
θ(h)

●● ●● ● ●
● ● ●
●●●
● ● ●
● ●
●● ●● ●
●● ● ● ● ●● ● ●
● ● ● ● ● ●●
●●●
● ●●●
● ● ● ●●● ● ●● ● ●● ● ●
● ●● ● ●
● ●
● ●● ● ● ●●●
● ● ●●● ● ● ● ● ●● ●● ● ●
●● ● ● ●●●● ●●● ●● ● ●
●●● ● ●●●●● ● ● ● ●● ● ● ● ●

● ●●

● ● ●●
●●
●● ●● ●●●
● ● ● ●●●● ● ● ●● ● ● ● ●● ●


●● ● ● ●● ●● ●● ●
●●
●●● ●●●● ● ●● ● ● ●● ● ● ● ●●
●●● ● ●●●
●● ●●●


●●●●
● ●● ●●

●●● ● ● ●●
● ●●●●●●● ● ●● ●●●
● ●
●●● ● ●● ● ●
●●
● ●
● ● ● ●
● ● ●
● ● ● ● ● ●
● ● ● ● ●
●● ●● ● ● ● ●
● ●● ●● ● ●●●●●● ● ●●●● ● ● ● ●●● ●
● ●● ● ● ●
1.4

● ● ● ●●●● ● ●● ● ● ●● ● ●●●
● ●●
●●
●●●● ● ●● ●●● ● ● ●● ● ●

● ●● ● ● ● ● ● ● ●
●● ● ●● ●●
● ●● ● ●● ●● ● ● ● ● ●
●● ● ● ●
●● ● ● ●
●● ● ● ● ● ●●
● ●●
● ● ●
● ●● ● ●
● ● ● ● ●
● ●


1.2
1.0

0 1 2 3 4

Figura 4.8: F -madograma Modelo 1.

Mediante el F -madograma se comprobará si el modelo es capaz de reprodu-


cir la estructura de dependencia espacial; con este gráfico se compararán las
estimaciones del coeficiente extremal con las predichas por el modelo ajustado.

Modelo 2

fmadogram(wind.frechet.M2,coord[,-3],M2.maxstab,which=c(.ext"),col=c("grey"))
fmadogram(wind.frechet.M2,coord[,-3],M2.maxstab,n.bins=200,which=c(.ext"),col=c(1,2),marge=.emp"

2.0
● ●


● ●

● ●


● ●●
1.8

● ● ●●
● ● ● ●
● ●● ● ●●●● ● ●
●● ● ●● ● ●● ● ●
● ● ●
● ●
● ●● ● ●● ● ●
● ●●●
●● ● ●● ● ● ●● ●
● ● ● ●●● ●
● ●● ● ●● ● ● ● ● ●●●●● ● ● ● ● ●
● ● ●● ● ●
●●● ● ●● ● ●
●● ● ● ● ● ●●
● ● ● ●● ● ●
●● ● ●● ● ●● ●
●●
● ●● ●
●●●● ●● ●●
● ●● ● ● ●●● ● ● ● ● ●● ●● ● ● ●
1.6

● ● ● ● ● ● ●● ● ●● ●●● ●
● ●● ● ● ●● ●● ●●●● ●● ● ●●
● ●●●●●● ● ● ● ● ●
●●● ●●
● ● ●●●● ● ●

●●●●●● ●
●●
●● ●
●●● ●● ●●
●● ● ● ●●
● ●● ● ● ● ●●
● ● ● ●● ● ●● ●●● ●●●● ● ●● ●● ●● ●●●●● ●● ●

● ●●● ● ●●
●● ●
●●● ● ● ●● ●●●

●● ●
●● ● ● ●● ● ●● ● ●
● ●● ●●● ● ● ●● ● ●● ●● ● ●
θ(h)

●● ●● ● ●
● ● ●●
●●
● ● ●
● ●
●● ●● ●
●● ● ● ● ●● ● ●
● ●● ● ● ●●●●
●● ●
●●●
● ● ●●● ● ●● ● ●● ● ●
● ●● ● ●
● ●
● ●● ● ● ●●●
● ● ●●● ●● ● ● ●● ●● ● ●
●● ● ● ●●●● ●●● ●● ● ●
●●● ● ●●●●● ● ● ● ●● ● ● ● ●

● ●●

● ● ●●
●●
●● ●● ●●●
● ● ● ●●●● ● ● ●● ● ● ● ●● ●


●● ● ● ●● ●● ●● ●
●●
●●● ●●●● ● ●● ● ● ●● ● ● ● ●●
●●● ● ●●●
●● ●●●


●●●●
● ●● ●●

●●● ● ● ●●
● ●●●●●●● ● ●● ●●●
● ●
●●● ● ●● ● ●
●●
● ●
● ● ● ●
● ● ●
● ● ● ●● ●● ● ● ● ●
●● ●● ● ● ● ●
● ●● ●● ● ●●●●●● ● ●●●● ● ● ● ●●● ●
● ●● ● ● ●
1.4

● ● ● ●●●● ● ●● ● ● ●● ● ●●●
● ●●
●●
●●●● ● ●● ●●● ● ● ●● ● ●

● ●● ● ● ● ● ● ● ●

● ● ●● ●●
● ●● ● ●● ●● ● ● ● ● ●
●● ● ● ●
●● ● ● ●
●● ● ● ● ● ●●
● ●●
● ● ●
● ●● ● ●
● ● ● ● ●
● ●


1.2
1.0

0 1 2 3 4

Figura 4.9: F -madograma.

Para realizar la selección del modelo se podrá usar el estadı́stico TIC (Criterio
de Información de Takeuchi); para los modelos calculados anteriormente se tiene:

TIC(M1.maxstab,M2.maxstab)
M1.maxstab M2.maxstab
-14857.13 75049.33

Según este criterio, se elegirı́a el Modelo 1, ya que es el que minimiza el resultado


del TIC; este resultado es previsible, ya que el Modelo 1 es más restrictivo que
el Modelo 2 al haber predefinido los parámetros.
Existen otros métodos de selección de modelos que penalizan el número de
parámetros, como el Estadı́stico de la Tasa de Verosimilitud. En R este test se
llevará a cabo mediante la función anova() usando como método la aproxima-
ción de Rotnitzky y Jewell RJ”.
Eigenvalue(s): 28.69 14.72 9.37
Analysis of Variance Table
MDf Deviance Df Chisq Pr(>sum lambda Chisq)
M2.maxstab 5 74417
M1.maxstab 8 -15503 3 89920 <2.2e-16 ***
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

4.3.4. Predicciones.
Las predicciones para los modelos espaciales se llevan a cabo mediante la función
predict, que ya se usó en la sección del ajuste del modelo GEV; sin embar-
go, habrá que tener en cuenta que estas predicciones se realizan localización a
localización, por lo que no se tendrá en cuenta la dependencia espacial.
Capı́tulo 5

Conclusión

El principal objetivo de este trabajo ha sido realizar una revisión de los princi-
pales enfoques en el análisis de valores extremos, centrándose principales en los
procesos max-estables y en la dimensión espacial.

En la última década han surgido muchos avances en materia de extremos, fun-


damentalmente en el campo de la geoestadı́stica a través del uso de extremos
espaciales, siendo los procesos max-estables una solución elegante para modeli-
zar estos valores.

Algunos de los últimos avances en el campo del ánalisis de los valores extremos
son los modelos bayesianos jerárquicos (Ribatet et al., 2012; Reich and Shaby,
2012; Thibaud et al., 2015), que proporcionan algunas mejoras en la definición
de las tendencias de superficie y modelos eficientes para las predicciones pun-
tuales; extremos espacio-temporales (Huser, 2013; Turkman and Pereira, 2010)
y análisis de valores extremos no estacionarios (Cheng y AghaKuochak, 2014).

Sin embargo, el desarrollo de esta teorı́a es relativamente reciente y aún siguen


surgiendo nuevos avances teóricos y se esperan progresos en el futuro.

61
Bibliografı́a

[1] S. A. Padoan A. C. Davison and M. Ribatet. Statistical modeling of spatial


extremes. Institute of Mathematical Statistics, Vol. 27, No. 2, 161-186,
2012.
[2] H. Akaike. Selected Papers of Hirotugu Akaike, chapter A new look at the
statistical model identification, pages 215–222. Springer, 1974.
[3] J. N. Bacro and C. Gaetan. Advances and Challenges in Space-time Mo-
delling of Natural., chapter 5 - A Review on Spatial Extreme Modelling,
pages 103–124. Springer-Verlag, 2012.
[4] L. Cheng and A. AghaKouchak. Nonstationary precipitation intensity-
duration-frequency curves for infrastucture design in a changing climate.
Scientific Reports, No. 4, 7093, 2014.
[5] S. Coles. An Introduction to Statistical Modeling of Extreme Values. Sprin-
ger, 2001.
[6] S. Coles and J. A. Tawn. Statistics of coastal flood prevention. Philosophical
Transactions of the Royal Society of London, No. 332, 457-476, 1990.
[7] S. Coles and J. A. Tawn. Modelling extreme multivariate events. Journal
of the Royal Statistical Society, No. 53, 377-392, 1991.
[8] A. C. Davison. Statistical Models. Cambridge Series in Statistical and
Probabilistic Mathematics. Cambridge University Press, 2003.
[9] L. de Haan. On Regular Variation and its Applications to the Weak Con-
vergence of Sample Extremes. PhD thesis, Mathematical Centre Tract, 32,
Amsterdam, 1970.
[10] L. de Haan. A spectral representations for max-stable processes. Annals
of Probability, No. 12, 1194-1204, 1984.
[11] L. de Haan and A. Ferreira. Extreme Value Theory. An Introduction. Sprin-
ger, 2006.
[12] R. Fisher and L. H. C. Tippet. Limiting forms of the frecuency distribution
of the largest or smallest member of a sample. Proceedings of the Cambridge
Philosophical Society. No. 24, 180-190, 1928.
[13] M. Fréchet. Sur la loi de probabilité de l’écart maximum. Ann.Soc. Math.
Polon. No. 6, 93-116, 1927.

63
[14] R. Huser and A. C. Davison. Space-time modelling of extreme events. Royal
Statistical Society, No. 76, 439-461, 2013.

[15] J. Pickand III. Statistical inference using extreme order statistics. Annals
of Statistic, No. 3, 119-131, 1975.

[16] M. A. Turkman K. F. Turkman and J. M. Pereira. Asymptotic models


and inference for extremes of spatio-temporal data. Extremes, No. 13, 375,
2010.

[17] T. W. Yee L. Thibaut and D. Arnaud. Expectation particle belief propa-


gation. Technical report, University of Oxford. Department of Statistics.,
2015.

[18] A. W. Ledford and J. A. Tawn. Statistics for near independence in multi-


variate extreme values. Biometrika, No. 83, 169-187, 1996.

[19] B. Lindsay. Composite likelihood methods. Contemporary Mathematics,


No. 80, 220-239, 1988.

[20] A. Guillou P. Naveau and D. Cooley. Modelling pairwise dependence of


maxima in space. Biometrika, No. 1, 1-17, 2009.

[21] M. D. Penrose. Semi-min-stable processes. Annals of Probability, No. 3,


1450-1463, 1992.

[22] B. J. Reich and B. A. Shaby. A hierarchical max-stable spatial model for


extreme precipitation. The Annals of Applied Statistics, No. 4, 1430-1451,
2012.

[23] R. D. Reiss and M. Thomas. Statistical Analysis of Extreme Values with


Applications to Insurance, Finance, Hydrology and Other Fields. Birkhau-
ser, 1997.

[24] M. Ribatet. Modelling spatial extremes with the spatialextreme packa-


ge. In The 9th International Conference on Extreme Value Analysis and
Application. University of Montpellier.

[25] M. Ribatet. A User’s Guide to the SpatialExtremes Package. École Poly-


technique Fédérale de Lausanne, 2009.

[26] M. Ribatet. Spatial extremes: Max-stable processes at work. Journal of


Société Française de Statistique, Vol .154,No. 2,156-177, 2013.

[27] M. Ribatet and M. Sedki. Extreme value copulas and max-stable processes.
Journal of Société Française de Statistique, Vol. 152, No. 3,138-150, 2012.

[28] M. Hofmann S. Aulbach, M. Falk and M. Zott. On max-stable processes


and the functional d-norm revisited. Extremes, Vol. 16, No. 3, 255-283,
2013.

[29] M. Schlather. Models for stationary max-stable random fields. Extremes,


Vol. 1, No.5, 33-44, 2002.
[30] M. Schlather and J. Tawn. A dependence measure for multivariate and
spatial extreme values: properties and inference. Biometrika, No. 90, 139-
156, 2003.

[31] J. Segers. Max-stable models for multivariate extremes. In Copula Models


and Dependence. Centre de Recherches Mathématiques, 2011.
[32] R. L. Smith. Max-stable processes and spatial extremes. Technical report,
1990.

[33] S. A. Stoev. Max-stable processes: Representations, ergodic properties and


statistical applications. Technical report, Department of Statistics, The
University of Michigan, 2010.
[34] C. Varin and P. Vidoni. A note on composite likelihood inference and model
selction. Biometrika. No. 92, 519-528, 2005.

[35] R. von Mises. la distribution de la plus grande de n valeurs. Revue


Mathématique de l’Union Interbalcanique, No. 1, 141-160, 1936.
[36] Y. Wang and S. A. Stoev. Conditional sampling for spspectral discrete
max-stable random fields. Advances in Applied Probability, No. 2, 461-
483, 2011.
[37] Z. Zhang. On approximating max-stable processes and constructing extre-
mal copula functions. Statistical Inference for Stochastic Processes, No. 12,
89-114, 2009.

También podría gustarte