MedialdeaVillanueva TFM
MedialdeaVillanueva TFM
MedialdeaVillanueva TFM
Extremos
Modelización Espacial
12 de Septiembre de 2016
Índice general
Resumen 5
3
4 ÍNDICE GENERAL
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
5. Conclusión 61
6 ÍNDICE GENERAL
Resumen
Mn = máx{X1 , ..., Xn },
n−1
fMn (z) = n (FMn (z)) f (z) . (1.2)
9
Mn − bn
Mn∗ = (1.3)
an
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
Teorema 1.1.2 Si existen sucesiones de constantes {an > 0} y {bn } tales que
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)
σ
k k
X zi − µ X zi − µ
log L (µ, σ) = −k log σ − − exp − . (1.10)
i=1
σ i=1
σ
µ − σ 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.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)
σ̃ = σ + ξ (u − µ) . (1.17)
k
1X
log L (σ, ξ) = −k log σ − yi , si ξ = 0. (1.19)
σ i=1
σ
(σ, ξ) → (τ, ξ) con τ = ,
ξ
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
r s
Mp,r,s = E {Y p [F (Y )] [1 − F (Y )] } ,
k s
1X j
M̂1,0,s = 1− yj,k .
k j=1 k+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
− 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
ξˆ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 ,
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.
σ
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.
σu = σu0 + ξ (u − u0 ) si ξ 6= 0. (1.25)
σ ∗ = σu − ξu.
ψ̂ ∼ Nr ψ; 5Tθ ψΣθ̂ 5θ ψ ,
V ar (σ ∗ ) ≈ 5σ ∗T V σ ∗ ,
donde
∂σ ∗ ∂σ ∗
∗T
5σ = , = [1, −u] .
∂σu ∂ξ
−1/ξ
x−u
P r {X > x|X > u} = 1 + ξ . (1.26)
σ
−1/ξ
xk − u 1
ζu 1 + ξ = . (1.27)
σ k
u + σ (kζu )ξ − 1
h i
para ξ 6= 0,
xk = ξ (1.28)
u + σ log (kζu ) para ξ = 0.
τ1 = mı́n {k : Xk > u} ,
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.
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
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
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
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:
M = {M1 , ..., Mp } .
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,
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
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
V (x1 , x2 ) = x−1 −1
1 + x2
Ejemplo 2
1/ξx1
x1 − µx1
x̃1 = 1 + ξx1 ,
σx1
1/ξx2
x2 − µx2
x̃2 = 1 + ξx2 ,
σx2
se obtiene
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
Ejemplo 1
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).
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:
n oi
1/i
Ci (u) = C1 u1 , ..., u1/i
p .
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,
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} .
xj
` (x) = (x1 + · · · xp ) D (ω1 , ..., ωp ) , ωj = .
x1 + · · · + xp
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
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] .
p
(ex1 , ..., exp ) = exp {−` (−x1 , ..., −xp )} , x ∈ (−∞, 0) ,
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 ) .
ˆi
1/xi
ˆ i Zi − µ̂i
Z̃i = 1 + xi
σ̂i
g (x, y) = {Vx (x, y) Vy (x, y) − Vxy (x, y)} exp {−V (x, y)} , x > 0, y > 0,
m
Y
L (θ) = g (z̃1,i , z̃2,i ) ,
i=1
Z = φ (Mx1 , Mx2 ) ,
GZ (z) = 1 − 1/N,
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]) .
( −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 )} ,
F (x1 , x2 ) ≈ G (x1 , x2 ) = exp {−V (x˜1 , x˜2 )} , x1 > ux1 , x2 > ux2 .
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
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
∂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
Extremos espaciales
3.1. Introducción
37
un conjunto R. Este proceso sera max-estable si todas sus distribuciones n-
dimensionales satisfacen la propiedad de max-estabilidad, es decir:
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
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
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.
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
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.2. F-madograma
1 h 2
i
γ (h) = E {W (r) − W (r + h)} , r, r + h ∈ R.
2
1
νF (h) = E [|F {Z (r + h)} − F {Z (r)} |] , r, r + h ∈ R1
2
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
3.3.3. λ-madograma
1 λ
E |F {Z (r1 )} − F 1−λ {Z (r2 )} | ,
νλ (r1 , r2 ) = ∀ λ ∈ [0, 1] .
2
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
(sim)
Ẑ (sim) (r∗) = máx φk (r∗) Yk .
k=1,...,p
2
X θ −
i,j θ̃ i,j
C (ψ) = ,
i<j s θ̃i,j
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
∂ 2 log f (y; ψ)
Z
H (ψ) = n g (y) dy
∂ψ∂ψ T
Z
∂ log f (y; ψ) ∂ log f (y; ψ)
J (ψ) = n .
∂ψ ∂ψ T
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 ,
n o p
X
W (κ0 ) = 2 ` κ̂, φ̂ − ` κ0 , φ̂κ0 → λi Xi , n → ∞,
i=1
−1
H −1 JH −1 H −1
κ κ
.
Capı́tulo 4
47
maps::map(xlim=c(0,9),ylim=c(50,54),fill=TRUE,col="green3")
points(coord,pch=15)
5. Predicciones.
● ●
1.7
1.7
● ●
● ●
● ●
● ●
● ●
●
●
●
● ●● ●
● ●
●
● ● ●
●● ●
●
● ● ● ●
1.6
● ● ●
● ● ●
1.6
● ●● ● ●● ●
● ●
Ext.Coeff
Ext.Coeff
● ● ●
●● ● ● ● ●
● ● ●● ● ●● ● ● ●
●
●● ● ●
● ● ● ●
● ● ● ●
● ● ● ● ●●● ● ●
●● ● ● ● ● ● ●
● ● ●● ● ●
●●● ● ● ●● ●
●
●●
●● ●●
● ● ● ● ●●
● ●
●
● ● ● ●
● ● ●● ● ● ●● ● ●
● ● ● ●● ● ● ● ●●
● ●
1.5
●
● ●● ● ●● ●
●●
●● ● ●●● ●● ●
● ● ●● ● ● ● ●● ●● ● ●● ●
●●
1.5
●● ● ● ●
●
●●●● ●●
●
● ●
● ● ●●● ●
●
●
●●
●
● ●● ● ● ●
●● ● ●
●● ●
● ● ● ● ●● ● ●● ●
● ●● ● ●● ●
● ●● ●●
● ●
● ● ● ●
● ●
●●● ●●●● ● ●
●
● ● ● ●
●● ●
● ●●● ● ●
● ● ●
● ● ●
1.4
● ●
● ●●●●
1.4
●
● ●
●
Distancia Distancia
●
1.7
● ● ●
● ● ●
● ●●
●
● ● ● ●
●● ● ●●
1.6
● ● ●●
●
● ● ● ●
●●
Ext.Coeff
● ● ●
● ●
● ● ●
● ● ●
● ● ● ●● ● ● ● ● ● ●
●
● ●
●● ●●
● ● ●● ● ● ●● ● ●
● ● ●
●● ● ● ●●
● ● ● ●
●
● ● ● ●
● ●● ●
1.5
●● ● ● ● ●
●
● ●● ● ●●
● ●
● ● ● ●
● ● ●●
● ●
●● ●
● ● ● ● ●
● ● ● ●
● ● ●
●● ●●
● ● ● ●
● ●
● ●
1.4
●
● ●
●
Distancia
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,ξ
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
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
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
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.1. Modelo 1
± 0.5
●
±1
●
± 1.5
● ● ● ●
± 0.05 ± 0.1 ± 0.15 ± 0.2
● ● ●●●●
±1 ±2 ±3 ±4 ±5 ±6
µ (x) = β0,µ + β1,µ lon (x) + β2,µ lat (x) + β3,µ lon (x) lat (x) ,
σ (x) = β0,σ
ξ (x) = β0,ξ
Equivalentemente en R.
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
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)
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
● ●
●
●
●
● ●
●
● ●
●
●
● ●●
1.8
● ● ●●
● ● ● ●
● ●● ● ●●●● ● ●
●● ● ●● ● ●● ● ●
● ● ●
● ●
● ●● ● ●● ● ●
● ●●●●● ● ●● ● ● ● ● ●
● ● ● ●●● ●
● ●● ● ●● ● ● ● ● ●●●●● ● ● ● ● ●
● ● ●● ● ●
●●● ● ●● ● ●
●● ● ● ● ● ●●
● ● ● ●● ● ●● ● ●●
●
● ●● ●
●●
● ●● ●
●●●● ●● ●●
● ●● ● ● ●●● ● ● ● ● ●● ●● ● ● ●
1.6
● ● ● ● ● ● ●● ● ●
● ●●● ●
● ●● ● ● ●● ●● ●● ●● ●● ● ●●
● ●●●●●● ● ● ● ● ●●●● ●
●● ● ●●●● ● ●
●
●●●●●● ●
●●
●● ●
●●● ●● ●●
●● ● ● ●●
● ●● ● ● ● ●●
● ● ● ●● ● ●● ●●● ●●●● ● ●● ● ● ●● ●●●●● ●● ●
●
● ●●● ● ●●
●● ●
●●● ● ●
●●●● ●●●
● ●
●● ● ● ●● ● ●● ● ●
● ●● ●●● ● ● ●● ● ●● ●● ● ●
θ(h)
●● ●● ● ●
● ● ●
●●●
● ● ●
● ●
●● ●● ●
●● ● ● ● ●● ● ●
● ● ● ● ● ●●
●●●
● ●●●
● ● ● ●●● ● ●● ● ●● ● ●
● ●● ● ●
● ●
● ●● ● ● ●●●
● ● ●●● ● ● ● ● ●● ●● ● ●
●● ● ● ●●●● ●●● ●● ● ●
●●● ● ●●●●● ● ● ● ●● ● ● ● ●
●
● ●●
●
● ● ●●
●●
●● ●● ●●●
● ● ● ●●●● ● ● ●● ● ● ● ●● ●
●
●
●● ● ● ●● ●● ●● ●
●●
●●● ●●●● ● ●● ● ● ●● ● ● ● ●●
●●● ● ●●●
●● ●●●
●
●
●●●●
● ●● ●●
●
●●● ● ● ●●
● ●●●●●●● ● ●● ●●●
● ●
●●● ● ●● ● ●
●●
● ●
● ● ● ●
● ● ●
● ● ● ● ● ●
● ● ● ● ●
●● ●● ● ● ● ●
● ●● ●● ● ●●●●●● ● ●●●● ● ● ● ●●● ●
● ●● ● ● ●
1.4
● ● ● ●●●● ● ●● ● ● ●● ● ●●●
● ●●
●●
●●●● ● ●● ●●● ● ● ●● ● ●
●
● ●● ● ● ● ● ● ● ●
●● ● ●● ●●
● ●● ● ●● ●● ● ● ● ● ●
●● ● ● ●
●● ● ● ●
●● ● ● ● ● ●●
● ●●
● ● ●
● ●● ● ●
● ● ● ● ●
● ●
●
●
1.2
1.0
0 1 2 3 4
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
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
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.
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).
61
Bibliografı́a
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.
[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.