Inferencia Bayesiana para Valores Extremos
Mario M. Pizarro
- Distribuciones de Valores Extremos
- Algoritmo de Metropolis-Hastings
- Estrategia Baseline Distribution para estimar los parámetros de las distribuciones de Valores Extremos
- Estrategia Improved Baseline Distribution para estimar los parámetros de la distribución Gumbel con distribuciones base: Gumbel, Exponencial y Normal
- Estrategia Informative Priors Baseline Distribution para estimar los parámetros de la distribución de Pareto Generalizada
- Algoritmo para estimar las Medidas de Riesgo Financieras (VaR y CVaR) con distribuciones base: Normal, Cauchy, Exponencial y Gamma
- Algoritmo de un Modelo Jerárquico Bayesiano
- Agradecimientos
Distribuciones de Valores Extremos
Distribución de Valores Extremos Generalizada
pGEV <- function(q,xi=0,sigma=1,mu=0,lower.tail=TRUE){
s <- (q-mu)/sigma
if(xi==0) out <- exp(-exp(-s)) else
out <- exp(-(1+xi*s)^(-1/xi))*as.numeric((1+xi*s)>0)
if(!lower.tail) out <- 1-out
out
}
dGEV <- function(x,xi=0,sigma=1,mu=0){
s <- (x-mu)/sigma
if(xi==0) out <- exp(-exp(-s))*exp(-s)/sigma else
out <- exp(-(1+xi*s)^(-1/xi))*(1+xi*s)^(-1/xi-1)/sigma*as.numeric((1+xi*s)>0)
out
}
rGEV <- function(n=1,xi=0,sigma=1,mu=0){
if(xi==0) out <- mu-sigma*log(-log(runif(n))) else
out <- mu-sigma/xi*(1-(-log(runif(n)))^(-xi))
out
}
qGEV <- function(p,xi=0,sigma=1,mu=0,lower.tail=TRUE){
if(!lower.tail) p <- 1-p
if(xi==0) out <- mu-sigma*log(-log(p)) else
out <- mu-sigma/xi*(1-(-log(p))^(-xi))
out
}Distribución Gumbel
pGumbel <- function(q,sigma=1,mu=0,lower.tail=TRUE){
out <- exp(-exp(-(q-mu)/sigma))
if(!lower.tail) out <- 1-out
out
}
dGumbel <- function(x,sigma=1,mu=0)
exp(-(x-mu)/sigma)*exp(-exp(-(x-mu)/sigma))/sigma
rGumbel <- function(n=1,sigma=1,mu=0)
mu - sigma*log(-log(runif(n)))
qGumbel <- function(p,sigma=1,mu=0,lower.tail=TRUE){
if(!lower.tail) p <- 1-p
mu - sigma*log(-log(p))
}Distribución de Pareto Generalizada
pGPD <- function(q,xi=0,sigma=1,lower.tail=TRUE){
s <- q/sigma
if(xi==0) out <- 1-exp(-s) else
out <- 1-(1+xi*s)^(-1/xi)
if(xi>=0) out <- out*as.numeric(s>=0) else
out <- out*as.numeric((s>=0)&&(s<=-1/xi))
if(!lower.tail) out <- 1-out
out
}
dGPD <- function(x,xi=0,sigma=1){
s <- x/sigma
if(xi==0) out <- exp(-s)/sigma else
out <- (1+xi*s)^(-1/xi-1)/sigma
if(xi>=0) out <- out*as.numeric(s>=0) else
out <- out*as.numeric((s>=0)&&(s<=-1/xi))
out
}
rGPD <- function(n=1,xi=0,sigma=1){
ifelse(xi==0,-sigma*log(runif(n)),
-sigma/xi*(1+runif(n)^(-xi)))
}
qGPD <- function(p,xi=0,sigma=1,lower.tail=TRUE){
if(!lower.tail) p <- 1-p
if(xi==0) out <- -sigma*log(1-p) else
out <- -sigma/xi*(1+(1-p)^(-xi))
}Algoritmo de Metropolis-Hastings
Para la distribución Gumbel
MHGumbel <- function(data,k=10,nthin=25,ndraw=1000,nburn=1/3*ndraw){
m <- length(data)
n <- m%/%k
Xmax <- apply(matrix(data,ncol=k,nrow=n),1,max)
sigmaI <- abs((mean(data)-median(data)))/(-digamma(1)+log(log(2)))+10^(-7)
muI <- mean(data)+digamma(1)*sigmaI
sdmu <- 1; sdsigma <- sigmaI/4
mu0 <- muI; sigma0 <- 4
lambda0 <- 4
mui <- muI+sigmaI*log(k); sigmai <- sigmaI
for(i in 1:nburn){
muc <- rnorm(1,mui,sdmu)
sigmac <- rnorm(1,sigmai,sdsigma)
while(sigmac<=0) sigmac <- rnorm(1,sigmai,sdsigma)
rmu <- n*(muc-mui)/sigmai + (mui-muc)/sigma0 + exp((mu0-mui)/sigma0) -
exp((mu0-muc)/sigma0) + sum(exp((mui-Xmax)/sigmai)-exp((muc-Xmax)/sigmai))
rsigma <- (n-1)*log((sigmai/sigmac)) + (sigmai^2-sigmac^2)/(2*lambda0^2) +
sum((Xmax-mui)*(1/sigmai-1/sigmac)) + sum(exp((mui-Xmax)/sigmai)-exp((mui-Xmax)/sigmac))
u <- log(runif(2))
if(u[1]<rmu) mui <- muc
if(u[2]<rsigma) sigmai <- sigmac
}
draw <- matrix(c(mui,sigmai),ncol=2,nrow=ndraw,byrow=TRUE)
for(i in 2:ndraw){
for(j in 1:nthin){
muc <- rnorm(1,mui,sdmu)
sigmac <- rnorm(1,sigmai,sdsigma)
while(sigmac<=0) sigmac <- rnorm(1,sigmai,sdsigma)
rmu <- n*(muc-mui)/sigmai + (mui-muc)/sigma0 + exp((mu0-mui)/sigma0) -
exp((mu0-muc)/sigma0) + sum(exp((mui-Xmax)/sigmai)-exp((muc-Xmax)/sigmai))
rsigma <- (n-1)*log((sigmai/sigmac)) + (sigmai^2-sigmac^2)/(2*lambda0^2) +
sum((Xmax-mui)*(1/sigmai-1/sigmac)) + sum(exp((mui-Xmax)/sigmai)-exp((mui-Xmax)/sigmac))
u <- log(runif(2))
if(u[1]<rmu) mui <- muc
if(u[2]<rsigma) sigmai <- sigmac
}
draw[i,] <- c(mui,sigmai)
}
print(paste("Location = ",mean(draw[,1])," Scale = ",mean(draw[,2]),sep=" "))
list(blockmax=Xmax,
MCMC=list(location=draw[,1],scale=draw[,2]))
}Para la distribución Exponencial
MHExp <- function(data,p=0.9,nthin=25,ndraw=1000,nburn=1/3*ndraw){
u <- as.numeric(quantile(data,p))
Xu <- data[data>u]-u
m <- length(Xu)
sigma <- 1/rgamma(ndraw,1+m,1+sum(Xu))
print(paste("Scale = ",mean(sigma),sep=" "))
list(threshold=u,
excess=Xu+u,
MCMC=list(scale=sigma))
}Para la distribución de Pareto Generalizada
MHGPDTail <- function(data,shape,p=0.9,nthin=25,ndraw=1000,nburn=1/3*ndraw){
u <- as.numeric(quantile(data,p))
Xu <- data[data>u]-u
m <- length(Xu)
if(shape=="light"){
MXu <- max(Xu)
k0 <- 0.15; delta0 <- sd(data)/k0
a1 <- 0.1; b1 <- 0.1; sdk <- 0.1
a2 <- 0.1; b2 <- 0.1; sddelta <- 0.1
ki <- k0; deltai <- delta0
for(i in 1:nburn){
kc <- rnorm(1,ki,sdk)
while(kc<0) kc <- rnorm(1,ki,sdk)
deltac <- rnorm(1,deltai,sddelta)
while(deltac<MXu) deltac <- rnorm(1,deltai,sddelta)
rk <- (a1-m-1)*log(kc/ki) + b1*(ki-kc) + (1/kc-1/ki)*sum(log(1-Xu/deltai))
rdelta <- (a2-m-1)*log(deltac/deltai) + b2*(deltai-deltac) +
(1/ki-1)*sum(log(1-Xu/deltac)) - (1/ki-1)*sum(log(1-Xu/deltai))
v <- log(runif(2))
if(v[1]<rk) ki <- kc
if(v[2]<rdelta) deltai <- deltac
}
draw <- matrix(c(-ki,ki*deltai),ncol=2,nrow=ndraw,byrow=TRUE)
for(i in 2:ndraw){
for(j in 1:nthin){
kc <- rnorm(1,ki,sdk)
while(kc<0) kc <- rnorm(1,ki,sdk)
deltac <- rnorm(1,deltai,sddelta)
while(deltac<MXu) deltac <- rnorm(1,deltai,sddelta)
rk <- (a1-m-1)*log(kc/ki) + b1*(ki-kc) + (1/kc-1/ki)*sum(log(1-Xu/deltai))
rdelta <- (a2-m-1)*log(deltac/deltai) + b2*(deltai-deltac) +
(1/ki-1)*sum(log(1-Xu/deltac)) - (1/ki-1)*sum(log(1-Xu/deltai))
v <- log(runif(2))
if(v[1]<rk) ki <- kc
if(v[2]<rdelta) deltai <- deltac
}
draw[i,] <- c(-ki,ki*deltai)
}
} else {
xi0 <- 1; sigma0 <- 100
alpha <- 10^(-6); cc <- 10^(-3); sdxi <- 0.1
a1 <- 0.1; b1 <- 0.1; sdsigma <- 5
xii <- xi0; sigmai <- sigma0
for(i in 1:nburn){
xic <- rnorm(1,xii,sdxi)
while(xic<cc) xic <- rnorm(1,xii,sdxi)
sigmac <- rnorm(1,sigmai,sdsigma)
while(sigmac<0) sigmac <- rnorm(1,sigmai,sdsigma)
rxi <- (alpha+1)*log(xii/xic) + (1+1/xii)*sum(log(1+xii*Xu/sigmai)) -
(1+1/xic)*sum(log(1+xic*Xu/sigmai))
rsigma <- (m+a1+1)*log(sigmai/sigmac) + b1*(1/sigmai-1/sigmac) +
(1+1/xii)*sum(log(1+xii*Xu/sigmai)) - (1+1/xii)*sum(log(1+xii*Xu/sigmac))
v <- log(runif(2))
if(v[1]<rxi) xii <- xic
if(v[2]<rsigma) sigmai <- sigmac
}
draw <- matrix(c(xii,sigmai),ncol=2,nrow=ndraw,byrow=TRUE)
for(i in 2:ndraw){
for(j in 1:nthin){
xic <- rnorm(1,xii,sdxi)
while(xic<cc) xic <- rnorm(1,xii,sdxi)
sigmac <- rnorm(1,sigmai,sdsigma)
while(sigmac<0) sigmac <- rnorm(1,sigmai,sdsigma)
rxi <- (alpha+1)*log(xii/xic) + (1+1/xii)*sum(log(1+xii*Xu/sigmai)) -
(1+1/xic)*sum(log(1+xic*Xu/sigmai))
rsigma <- (m+a1+1)*log(sigmai/sigmac) + b1*(1/sigmai-1/sigmac) +
(1+1/xii)*sum(log(1+xii*Xu/sigmai)) - (1+1/xii)*sum(log(1+xii*Xu/sigmac))
v <- log(runif(2))
if(v[1]<rxi) xii <- xic
if(v[2]<rsigma) sigmai <- sigmac
}
draw[i,] <- c(xii,sigmai)
}
}
print(paste("Shape = ",mean(draw[,1])," Scale = ",mean(draw[,2]),sep=" "))
list(threshold=u,
excess=Xu+u,
MCMC=list(shape=draw[,1],scale=draw[,2]))
}MHGPD <- function(data,p=0.9,nthin=25,ndraw=1000,nburn=1/3*ndraw){
shape <- readline(prompt = "Indique si desea estimar cola Exp (exp), cola ligera (light) o cola pesada (heavy): ")
if(shape=="exp"){
MHExp(data,p,nthin,ndraw,nburn)
} else{
MHGPDTail(data,shape,p,nthin,ndraw,nburn)
}
}Estrategia Baseline Distribution para estimar los parámetros de las distribuciones de Valores Extremos
Para las distribuciones Gumbel, Exponencial y Normal con el objetivo de estimar los parámetros de la distribución Gumbel
BDM <- function(data,k=10,ndraw=1000){
m <- length(data)
dist <- readline(prompt = "Indique si desea emplear la distribución Gumbel (gumbel), Exponencial (exp) o Normal (norm): ")
switch(dist,
"gumbel" = {
out <- MHGumbel(data,k=1,ndraw=ndraw)
draw <- matrix(c(out$MCMC$location+out$MCMC$scale*log(k),out$MCMC$scale),ncol=2)
print(paste("Location Block Maxima = ",mean(draw[,1])," Scale Block Maxima = ",mean(draw[,2]),sep=" "))
return(list(MCMC=list(location=out$MCMC$location,scale=out$MCMC$scale,locationBlockMax=draw[,1],scaleBlockMax=draw[,2])))
},
"exp" = {
lambda <- rgamma(ndraw,1+m,1+sum(data))
draw <- matrix(c(log(k)/lambda,1/lambda),ncol=2)
print(paste("Rate =",mean(lambda),sep=""))
print(paste("Location = ",mean(draw[,1])," Scale = ",mean(draw[,2]),sep=" "))
return(list(MCMC=list(rate=lambda,location=draw[,1],scale=draw[,2])))
},
"norm" = {
mu <- rnorm(ndraw,m*mean(data)/(m+sd(data)),sqrt(sd(data)^2/(m+sd(data))))
sigma <- sqrt(1/rgamma(ndraw,m/2+1,sum((data-mean(data))^2/2)+1))
draw <- matrix(numeric(),ncol=2,nrow=ndraw)
draw[,1] <- mu + sigma*(sqrt(2*log(k))-(log(log(k))+log(4*pi))/(2*sqrt(2*log(k))))
draw[,2] <- sigma/sqrt(2*log(k))
print(paste("Mean =",mean(mu)," Sd = ",mean(sigma),sep=" "))
print(paste("Location = ",mean(draw[,1])," Scale = ",mean(draw[,2]),sep=" "))
return(list(MCMC=list(mean=mu,sd=sigma,location=draw[,1],scale=draw[,2])))
})
}Para las distribuciones Estables: Normal, Cauchy, y Lévy con el fin de estimar los parámetros de la distribución de Pareto Generalizada
BDMStable <- function(data,p=0.9,nthin=25,ndraw=1000,nburn=1/3*ndraw){
n <- length(data)
dist <- readline(prompt = "Indique si desea estimar la distribución Normal (norm), Cauchy (cauchy) o Lévy (levy): ")
draw <- matrix(numeric(),ncol=2,nrow=ndraw)
switch(dist,
"norm" = {
scale <- sqrt(1/rgamma(ndraw,1+n/2,1+sum(data^2)/2))
draw[,1] <- rep(-0.7+0.61*p,ndraw)
draw[,2] <- (0.34+3.18*(1-p)-12.4*(1-p)^2)*scale
},
"levy" = {
scale <- rgamma(ndraw,1+n/2,1+sum(1/data)/2)
draw[,1] <- rep(2,ndraw)
draw[,2] <- 4/(pi*(1-p)^2)*scale
},
"cauchy" = {
a <- 1; b <- 1; sdscale <- 0.1
scalei <- as.numeric((quantile(data,3/4)-quantile(data,1/4))/2)
for(i in 1:nburn){
scalec <- rnorm(1,scalei,sdscale)
while(scalec<0) scalec <- rnorm(1,scalei,sdscale)
rscale <- (a-n-1)*log(scalec/scalei) - b*(scalec-scalei) -
sum(log(1+data^2/scalec^2)) + sum(log(1+data^2/scalei^2))
if(log(runif(1))<rscale) scalei <- scalec
}
scale <- scalei
draw[1,] <- c(1,1/(pi*(1-p))*scalei)
for(i in 2:ndraw){
for(j in 1:nthin){
scalec <- rnorm(1,scalei,sdscale)
while(scalec<0) scalec <- rnorm(1,scalei,sdscale)
rscale <- (a-n-1)*log(scalec/scalei) - b*(scalec-scalei) -
sum(log(1+data^2/scalec^2)) + sum(log(1+data^2/scalei^2))
if(log(runif(1))<rscale) scalei <- scalec
}
scale <- c(scale,scalei)
draw[i,] <- c(1,1/(pi*(1-p))*scalei)
}
}
)
print(paste("Scale of ",dist," = ",mean(scale)))
print(paste("Shape = ",mean(draw[,1])," Scale = ",mean(draw[,2]),sep=" "))
list(MCMC=list(scaledist=scale,shape=draw[,1],scale=draw[,2]))
}Estrategia Improved Baseline Distribution para estimar los parámetros de la distribución Gumbel con distribuciones base: Gumbel, Exponencial y Normal
IBDM <- function(data,k=10,nthin=25,ndraw=1000,nburn=1/3*ndraw){
m <- length(data)
n <- m%/%k
Xmax <- apply(matrix(data,ncol=k,nrow=n),1,max)
dist <- readline(prompt = "Indique si desea emplear la distribución Gumbel (gumbel), Exponencial (exp) o Normal (norm): ")
if(dist=="gumbel") sigmaI <- abs((mean(data)-median(data)))/(-digamma(1)+log(log(2)))+10^(-7)
switch(dist,
"gumbel" = {
muI <- mean(data)+digamma(1)*sigmaI
mui <- muI+sigmaI*log(k)
sigmai <- sigmaI
},
"exp" = {
lambdaI <- rgamma(1,1+m,1+sum(data))
mui <- log(n)/lambdaI
sigmai <- 1/lambdaI
},
"norm" = {
muI <- rnorm(1,m*mean(data)/(m+sd(data)),sqrt(sd(data)^2/(m+sd(data))))
sigmaI <- sqrt(1/rgamma(1,m/2+1,sum((data-mean(data))^2/2)+1))
mui <- muI+sigmaI*(sqrt(2*log(k))-(log(log(k))+log(4*pi))/(2*sqrt(2*log(k))))
sigmai <- sigmaI/sqrt(2*log(k))
}
)
for(i in 1:nburn){
switch(dist,
"gumbel" = {
muc <- rnorm(1,mui,1)
sigmac <- rnorm(1,sigmai,sigmaI/4)
while(sigmac<=0) sigmac <- rnorm(1,sigmai,sigmaI/4)
muc <- muc + sigmac*log(k)
},
"exp" = {
lambdac <- rgamma(1,1+m,1+sum(data))
muc <- log(k)/lambdac
sigmac <- 1/lambdac
},
"norm" = {
muc <- rnorm(1,m*mean(data)/(m+sd(data)),sqrt(sd(data)^2/(m+sd(data))))
sigmac <- sqrt(1/rgamma(1,m/2+1,sum((data-mean(data))^2/2)+1))
muc <- muc+sigmac*(sqrt(2*log(k))-(log(log(k))+log(4*pi))/(2*sqrt(2*log(k))))
sigmac <- sigmac/sqrt(2*log(k))
}
)
r <- n*(log(sigmai)-log(sigmac)) +
sum(-exp(-(Xmax-muc)/sigmac) + exp(-(Xmax-mui)/sigmai) -
(Xmax-muc)/sigmac + (Xmax-mui)/sigmai)
if(log(runif(1))<r){
mui <- muc
sigmai <- sigmac
}
}
draw <- matrix(c(mui,sigmai),ncol=2,nrow=ndraw,byrow=TRUE)
for(i in 2:ndraw){
for(j in 1:nthin){
switch(dist,
"gumbel" = {
muc <- rnorm(1,mui,1)
sigmac <- rnorm(1,sigmai,sigmaI/4)
while(sigmac<=0) sigmac <- rnorm(1,sigmai,sigmaI/4)
muc <- muc + sigmac*log(k)
},
"exp" = {
lambdac <- rgamma(1,1+m,1+sum(data))
muc <- log(k)/lambdac
sigmac <- 1/lambdac
},
"norm" = {
muc <- rnorm(1,m*mean(data)/(m+sd(data)),sqrt(sd(data)^2/(m+sd(data))))
sigmac <- sqrt(1/rgamma(1,m/2+1,sum((data-mean(data))^2/2)+1))
muc <- muc+sigmac*(sqrt(2*log(k))-(log(log(k))+log(4*pi))/(2*sqrt(2*log(k))))
sigmac <- sigmac/sqrt(2*log(k))
}
)
r <- n*(log(sigmai)-log(sigmac)) +
sum(-exp(-(Xmax-muc)/sigmac) + exp(-(Xmax-mui)/sigmai) -
(Xmax-muc)/sigmac + (Xmax-mui)/sigmai)
if(log(runif(1))<r){
mui <- muc
sigmai <- sigmac
}
}
draw[i,] <- c(mui,sigmai)
}
print(paste("Location = ",mean(draw[,1])," Scale = ",mean(draw[,2]),sep=" "))
list(blockmax=Xmax,
MCMC=list(location=draw[,1],scale=draw[,2]))
}Estrategia Informative Priors Baseline Distribution para estimar los parámetros de la distribución de Pareto Generalizada
Distribución base: distribuciones estables Normal, Cauchy y Lévy
IPBDMStable <- function(data,dist,p=0.9,nthin=25,ndraw=1000,nburn=1/3*ndraw){
n <- length(data)
u <- as.numeric(quantile(data,p))
Xu <- data[data>u] - u
m <- length(Xu)
switch(dist,
"norm" = {
MXu <- max(Xu)
muxi0 <- -0.7 + 0.61*p; sdxi0 <- 0.03
musigma0 <- (0.34+3.18*(1-p)-12.4*(1-p)^2)*sd(data)
c1 <- -46.24; c2 <- 83.55; c3 <- -41.58
},
"cauchy" = {
muxi0 <- 1; sdxi0 <- 0.065
c1 <- 323.57; c2 <- -588.51; c3 <- 266.13
a <- 1; b <- 1
scalei <- as.numeric((quantile(data,3/4)-quantile(data,1/4))/2)
musigma0 <- 1/(pi*(1-p))*scalei
sdscale <- 0.1
},
"levy" = {
muxi0 <- 2; sdxi0 <- 0.1
c1 <- 500.2; c2 <- -900.9; c3 <- 408.2
scalei <- rgamma(1,1+n/2,1+sum(1/data)/2)
musigma0 <- 4/(pi*(1-p)^2)*scalei
}
)
sdsigma0 <- exp(c1*p^2+c2*p+c3)
sdxi <- 0.1; sdsigma <- 0.1
xii <- rnorm(1,muxi0,sdxi0)
sigmai <- rnorm(1,musigma0,sdsigma0)
if(dist=="norm"){
while(MXu>=(-sigmai/xii)||sigmai<=0||xii>=0){
xii <- rnorm(1,muxi0,sdxi0)
sigmai <- rnorm(1,musigma0,sdsigma0)
}
} else {
while(xii<0) xii <- rnorm(1,muxi0,sdxi0)
while(sigmai<=0) sigmai <- rnorm(1,musigma0,sdsigma0)
}
for(i in 1:nburn){
xic <- rnorm(1,xii,sdxi)
sigmac <- rnorm(1,sigmai,sdsigma)
if(dist=="norm"){
while(MXu>=(-sigmac/xic)||sigmac<=0||xic>=0){
xic <- rnorm(1,xii,sdxi)
sigmac <- rnorm(1,sigmai,sdsigma)
}
} else {
while(xic<0) xic <- rnorm(1,xii,sdxi)
while(sigmac<=0) sigmac <- rnorm(1,sigmai,sdsigma)
}
switch (dist,
"norm" = {
scale <- sqrt(1/rgamma(1,1+n/2,1+sum(data^2)/2))
musigma0 <- (0.34+3.18*(1-p)-12.4*(1-p)^2)*scale
},
"cauchy" = {
scale <- rnorm(1,scalei,sdscale)
while(scale<=0) scale <- rnorm(1,scalei,sdscale)
rscale <- (a-n-1)*log(scale/scalei) - b*(scale-scalei) -
sum(log(1+data^2/scale^2))+sum(log(1+data^2/scalei^2))
if(log(runif(1))<rscale) scalei <- scale
musigma0 <- 1/(pi*(1-p))*scalei
},
"levy" = {
scale <- rgamma(1,1+n/2,1+sum(1/data)/2)
musigma0 <- 4/(pi*(1-p)^2)*scale
}
)
if(dist=="norm"){
r <- m*log(sigmai/sigmac) + 0.5*((xii-muxi0)^2 - (xic-muxi0)^2)/sdxi0^2 +
0.5*((sigmai-musigma0)^2 - (sigmac-musigma0)^2)/sdsigma0^2 -
(1+1/xic)*sum(log(1+xic*Xu/sigmac)) + (1+1/xii)*sum(log(1+xii*Xu/sigmai))
if(log(runif(1))<r){
xii <- xic
sigmai <- sigmac
}
} else {
rxi <- 0.5*((xii-muxi0)^2 - (xic-muxi0)^2)/sdxi0^2 -
(1+1/xic)*sum(log(1+xic*Xu/sigmai)) +
(1+1/xii)*sum(log(1+xii*Xu/sigmai))
rsigma <- m*log(sigmai/sigmac) + 0.5*(sigmai - musigma0)^2/sdsigma0^2 -
0.5*(sigmac-musigma0)^2/sdsigma0^2 -
(1+1/xii)*sum(log(1+xii*Xu/sigmac)) +
(1+1/xii)*sum(log(1+xii*Xu/sigmai))
v <- log(runif(2))
if(v[1]<rxi) xii <- xic
if(v[2]<rsigma) sigmai <- sigmac
}
}
draw <- matrix(c(xii,sigmai),ncol=2,nrow=ndraw,byrow=TRUE)
for(i in 2:ndraw){
for(j in 1:nthin){
xic <- rnorm(1,xii,sdxi)
sigmac <- rnorm(1,sigmai,sdsigma)
if(dist=="norm"){
while(MXu>=(-sigmac/xic)||sigmac<=0||xic>=0){
xic <- rnorm(1,xii,sdxi)
sigmac <- rnorm(1,sigmai,sdsigma)
}
} else {
while(xic<0) xic <- rnorm(1,xii,sdxi)
while(sigmac<=0) sigmac <- rnorm(1,sigmai,sdsigma)
}
switch (dist,
"norm" = {
scale <- sqrt(1/rgamma(1,1+n/2,1+sum(data^2)/2))
musigma0 <- (0.34+3.18*(1-p)-12.4*(1-p)^2)*scale
},
"cauchy" = {
scale <- rnorm(1,scalei,sdscale)
while(scale<=0) scale <- rnorm(1,scalei,sdscale)
rscale <- (a-n-1)*log(scale/scalei) - b*(scale-scalei) -
sum(log(1+data^2/scale^2))+sum(log(1+data^2/scalei^2))
if(log(runif(1))<rscale) scalei <- scale
musigma0 <- 1/(pi*(1-p))*scalei
},
"levy" = {
scale <- rgamma(1,1+n/2,1+sum(1/data)/2)
musigma0 <- 4/(pi*(1-p)^2)*scale
}
)
if(dist=="norm"){
r <- m*log(sigmai/sigmac) + 0.5*((xii-muxi0)^2 - (xic-muxi0)^2)/sdxi0^2 +
0.5*((sigmai-musigma0)^2 - (sigmac-musigma0)^2)/sdsigma0^2 -
(1+1/xic)*sum(log(1+xic*Xu/sigmac)) + (1+1/xii)*sum(log(1+xii*Xu/sigmai))
if(log(runif(1))<r){
xii <- xic
sigmai <- sigmac
}
} else {
rxi <- 0.5*((xii-muxi0)^2 - (xic-muxi0)^2)/sdxi0^2 -
(1+1/xic)*sum(log(1+xic*Xu/sigmai)) +
(1+1/xii)*sum(log(1+xii*Xu/sigmai))
rsigma <- m*log(sigmai/sigmac) + 0.5*(sigmai - musigma0)^2/sdsigma0^2 -
0.5*(sigmac-musigma0)^2/sdsigma0^2 -
(1+1/xii)*sum(log(1+xii*Xu/sigmac)) +
(1+1/xii)*sum(log(1+xii*Xu/sigmai))
v <- log(runif(2))
if(v[1]<rxi) xii <- xic
if(v[2]<rsigma) sigmai <- sigmac
}
}
draw[i,] <- c(xii,sigmai)
}
print(paste("Shape = ",mean(draw[,1])," Scale = ",mean(draw[,2]),sep=" "))
list(threshold=u,
excess=Xu+u,
MCMC=list(shape=draw[,1],scale=draw[,2]))
}Distribución base: distribución Exponencial
IPBDMExp <- function(data,p=0.9,nthin=25,ndraw=1000,nburn=1/3*ndraw){
n <- length(data)
u <- as.numeric(quantile(data,p))
Xu <- data[data>u] - u
m <- length(Xu)
sigmai <- mean(data)
for(i in 1:nburn){
lambdai <- rgamma(1,1+n,1+sum(data))
sigmac <- 1/lambdai
rsigma <- m*log(sigmai/sigmac) + sum(Xu)*(1/sigmai-1/sigmac)
if(log(runif(1))<rsigma) sigmai <- sigmac
}
draw <- sigmai
for(i in 2:ndraw){
for(j in 1:nthin){
lambdai <- rgamma(1,1+n,1+sum(data))
sigmac <- 1/lambdai
rsigma <- m*log(sigmai/sigmac) + sum(Xu)*(1/sigmai-1/sigmac)
if(log(runif(1))<rsigma) sigmai <- sigmac
}
draw <- c(draw,sigmai)
}
print(paste("Scale = ",mean(draw),sep=" "))
list(threshold=u,
excess=Xu+u,
MCMC=list(scale=draw))
}Distribución base: distribución Gamma
IPBDMGamma <- function(data,p=0.9,nthin=25,ndraw=1000,nburn=1/3*ndraw){
n <- length(data)
u <- as.numeric(quantile(data,p))
Xu <- data[data>u] - u
m <- length(Xu)
alphai <- mean(data)^2/var(data); sdalpha <- 0.5
betai <- mean(data)/var(data); sdbeta <- 0.5
xii <- -0.032+0.014/alphai
sigmai <- 1/betai*(0.5+0.5*sqrt(alphai))
for(i in 1:nburn){
alphac <- rnorm(1,alphai,sdalpha)
betac <- rnorm(1,betai,sdbeta)
while(alphac<0||betac<0){
alphac <- rnorm(1,alphai,sdalpha)
betac <- rnorm(1,betai,sdbeta)
}
rB <- -n*log(gamma(alphac)) + (n*alphac-1)*log(betac) +
(alphac-1)*sum(log(data)) - betac*sum(data) +
0.5*log(alphac*trigamma(alphac)-1) +
n*log(gamma(alphai)) - (n*alphai-1)*log(betai) -
(alphai-1)*sum(log(data)) + betai*sum(data) -
0.5*log(alphai*trigamma(alphai)-1)
if(log(runif(1))<rB){
alphai <- alphac
betai <- betac
}
xic <- -0.032+0.014/alphai
sigmac <- 1/betai*(0.5+0.5*sqrt(alphai))
rIP <- -m*log(sigmac) - (1+1/xic)*sum(log(1+xic*Xu/sigmac)) +
m*log(sigmai) + (1+1/xii)*sum(log(1+xii*Xu/sigmai))
if(log(runif(1))<rIP){
xii <- xic
sigmai <- sigmac
}
}
draw <- matrix(c(xii,sigmai),ncol=2,nrow=ndraw,byrow=TRUE)
for(i in 2:ndraw){
for(j in 1:nthin){
alphac <- rnorm(1,alphai,sdalpha)
betac <- rnorm(1,betai,sdbeta)
while(alphac<0||betac<0){
alphac <- rnorm(1,alphai,sdalpha)
betac <- rnorm(1,betai,sdbeta)
}
rB <- -n*log(gamma(alphac)) + (n*alphac-1)*log(betac) +
(alphac-1)*sum(log(data)) - betac*sum(data) +
0.5*log(alphac*trigamma(alphac)-1) +
n*log(gamma(alphai)) - (n*alphai-1)*log(betai) -
(alphai-1)*sum(log(data)) + betai*sum(data) -
0.5*log(alphai*trigamma(alphai)-1)
if(log(runif(1))<rB){
alphai <- alphac
betai <- betac
}
xic <- -0.032+0.014/alphai
sigmac <- 1/betai*(0.5+0.5*sqrt(alphai))
rIP <- -m*log(sigmac) - (1+1/xic)*sum(log(1+xic*Xu/sigmac)) +
m*log(sigmai) + (1+1/xii)*sum(log(1+xii*Xu/sigmai))
if(log(runif(1))<rIP){
xii <- xic
sigmai <- sigmac
}
}
draw[i,] <- c(xii,sigmai)
}
print(paste("Shape = ",mean(draw[1,])," Scale = ",mean(draw[,2]),sep=" "))
list(threshold=u,
excess=Xu+u,
MCMC=list(shape=draw[,1],scale=draw[,2]))
}Algoritmo para cualquier distribución base de las mencionadas anteriormente
IPBDM <- function(data,p=0.9,nthin=25,ndraw=1000,nburn=1/3*ndraw){
dist <- readline(prompt = "Indique si desea estimar la cola de la distribución Normal (norm), Cauchy (cauchy), Lévy (levy), Exponencial (Exp) o Gamma (gamma): ")
if(dist=="exp"){
IPBDMExp(data,p,nthin,ndraw,nburn)
} else if(dist=="gamma"){
IPBDMGamma(data,p,nthin,ndraw,nburn)
} else {
IPBDMStable(data,dist,p,nthin,ndraw,nburn)
}
}Algoritmo para estimar las Medidas de Riesgo Financieras (VaR y CVaR) con distribuciones base: Normal, Cauchy, Exponencial y Gamma
RiskMeasures <- function(data,p=0.9,prisk=0.95,nthin=25,ndraw=1000,nburn=1/3*ndraw){
dist <- readline(prompt = "Indique si desea estimar la cola de la distribución Normal (norm), Cauchy (cauchy), Exponencial (Exp) o Gamma (gamma): ")
u <- as.numeric(quantile(data,p))
ptail <- (1-prisk)/(1-p)
switch(dist,
"exp" = {
print("MH")
outMH <- MHExp(data,p,nthin,ndraw,nburn)
drawMH <- matrix(c(-log(1-ptail)*outMH$MCMC$scale + u, (1-log(1-ptail))*outMH$MCMC$scale + u), ncol=2,nrow=ndraw)
print(paste("VaR = ",mean(drawMH[,1])," CVaR = ",mean(drawMH[,2]),sep=" "))
print("IPBDM")
outIP <- IPBDMExp(data,p,nthin,ndraw,nburn)
drawIP <- matrix(c(-log(1-ptail)*outIP$MCMC$scale + u, (1-log(1-ptail))*outIP$MCMC$scale + u),ncol=2,nrow=ndraw)
print(paste("VaR = ",mean(drawIP[,1])," CVaR = ",mean(drawIP[,2]),sep=" "))
list(threshold=u,
MCMC=list(VaRMH=drawMH[,1],CVaRMH=drawMH[,2],
VaRIP=drawIP[,1],CVaRIP=drawIP[,2]))
},
"gamma" = {
print("MH")
outMH <- MHGPDTail(data,"light",p,nthin,ndraw,nburn)
drawMH <- matrix(c(outMH$MCMC$scale/outMH$MCMC$shape*((1-ptail)^(-outMH$MCMC$shape)-1) + u,outMH$MCMC$scale/outMH$MCMC$shape*((1-ptail)^(-outMH$MCMC$shape)/(1-outMH$MCMC$shape)-1) + u),ncol=2,nrow=ndraw)
print(paste("VaR = ",mean(drawMH[,1])," CVaR = ",mean(drawMH[,2]),sep=" "))
print("IPBDM")
outIP <- IPBDMGamma(data,p,nthin,ndraw,nburn)
drawIP <- matrix(c(outIP$MCMC$scale/outIP$MCMC$shape*((1-ptail)^(-outIP$MCMC$shape)-1) + u,outIP$MCMC$scale/outIP$MCMC$shape*((1-ptail)^(-outIP$MCMC$shape)/(1-outIP$MCMC$shape)-1) + u),ncol=2,nrow=ndraw)
print(paste("VaR = ",mean(drawIP[,1])," CVaR = ",mean(drawIP[,2]),sep=" "))
list(threshold=u,
MCMC=list(VaRMH=drawMH[,1],CVaRMH=drawMH[,2],
VaRIP=drawIP[,1],CVaRIP=drawIP[,2]))
},
"norm" = {
print("MH")
outMH <- MHGPDTail(data,"light",p,nthin,ndraw,nburn)
drawMH <- matrix(c(outMH$MCMC$scale/outMH$MCMC$shape*((1-ptail)^(-outMH$MCMC$shape)-1) + u,outMH$MCMC$scale/outMH$MCMC$shape*((1-ptail)^(-outMH$MCMC$shape)/(1-outMH$MCMC$shape)-1) + u),ncol=2,nrow=ndraw)
print(paste("VaR = ",mean(drawMH[,1])," CVaR = ",mean(drawMH[,2]),sep=" "))
print("IPBDM")
outIP <- IPBDMStable(data,dist,p,nthin,ndraw,nburn)
drawIP <- matrix(c(outIP$MCMC$scale/outIP$MCMC$shape*((1-ptail)^(-outIP$MCMC$shape)-1) + u,outIP$MCMC$scale/outIP$MCMC$shape*((1-ptail)^(-outIP$MCMC$shape)/(1-outIP$MCMC$shape)-1) + u),ncol=2,nrow=ndraw)
print(paste("VaR = ",mean(drawIP[,1])," CVaR = ",mean(drawIP[,2]),sep=" "))
list(threshold=u,
MCMC=list(VaRMH=drawMH[,1],CVaRMH=drawMH[,2],
VaRIP=drawIP[,1],CVaRIP=drawIP[,2]))
},
"cauchy" = {
print("MH")
outMH <- MHGPDTail(data,"heavy",p,nthin,ndraw,nburn)
drawMH <- outMH$MCMC$scale/outMH$MCMC$shape*((1-ptail)^(-outMH$MCMC$shape)-1) + u
print(paste("VaR = ",mean(drawMH),sep=" "))
print("IPBDM")
outIP <- IPBDMStable(data,dist,p,nthin,ndraw,nburn)
drawIP <- outIP$MCMC$scale/outIP$MCMC$shape*((1-ptail)^(-outIP$MCMC$shape)-1) + u
print(paste("VaR = ",mean(drawIP),sep=" "))
list(threshold=u,
MCMC=list(VaRMH=drawMH,
VaRIP=drawIP))
}
)
}Algoritmo de un Modelo Jerárquico Bayesiano
CopEstimates <- function(u,distance,cop.param,nT){
nS <- nrow(distance)
Ci <- cop.param[1]*exp(-distance/cop.param[2])
diag(Ci) <- rep(1,nS)
c0c <- runif(1)
c1c <- runif(1,0,1000)
Cc <- c0c*exp(-distance/c1c)
diag(Cc) <- rep(1,nS)
r <- -0.5*nT*log(det(Cc)) - 0.5*sum(sapply(1:nT,function(k) u[k,]%*%solve(Cc)%*%u[k,])) +
0.5*nT*log(det(Ci)) + 0.5*sum(sapply(1:nT,function(k) u[k,]%*%solve(Ci)%*%u[k,]))
if(log(runif(1))<r){
cop.param[1] <- c0c
cop.param[2] <- c1c
}
cop.param
}AlphaEstimates <- function(param,W,SigmaAlpha,tau,X){
library(mvtnorm)
B <- solve(t(X)%*%X/tau^2 + solve(SigmaAlpha))
b <- 1/tau^2*t(X)%*%(param-W)
if(length(param)==1)
rnorm(1,B*b,sqrt(B))
else
t(rmvnorm(1,B%*%b,B))
}WEstimates <- function(param,alpha,Sigma,tau,X){
library(mvtnorm)
nS <- length(param)
Id <- diag(1,nrow=nS,ncol=nS)
B <- solve(Id/tau^2 + solve(Sigma))
b <- 1/tau^2 * (param - X%*%alpha)
t(rmvnorm(1,B%*%b,B))
}BetaEstimates <- function(distance,param,W,Sigma,beta){
nS <- nrow(distance)
a0 <- 2; b0 <- 1
a1 <- 4; b1 <- 1/0.01; sdbeta <- 0.5
Di <- exp(-distance/beta[2])
a <- a0 + 0.5*nS
b <- b0 + 0.5*t(W)%*%solve(Di)%*%W
beta[1] <- 1/rgamma(1,a,b)
beta1c <- rnorm(1,beta[2],sdbeta)
Dc <- exp(-distance/beta1c)
r <- (a1-1)*log(beta1c/beta[2]) - 1/b1*(beta1c-beta[2]) -
0.5*log(det(Dc)) + 0.5*log(det(Di)) -
0.5/beta[1]*t(W)%*%solve(Dc)%*%W + 0.5/beta[1]*t(W)%*%solve(Di)%*%W
if(log(runif(1))<r) beta[2] <- beta1c
beta
}TauEstimates <- function(param,W,alpha,X){
nS <- length(param)
a0 <- 1; b0 <- 1
a <- a0 + 0.5*nS
b <- b0 + 0.5*t(param-X%*%alpha-W)%*%(param-X%*%alpha-W)
sqrt(1/rgamma(1,a,b))
}ParamEstimates <- function(name.param,sd,location,scale,shape,X,alpha,W,tau,data,u,C){
library(mvtnorm)
library(evd)
nS <- ncol(data)
nT <- nrow(data)
Sigma <- diag(sd,ncol=nS,nrow=nS)
switch(name.param,
"location" = {
param <- location
paramc <- t(rmvnorm(1,param,Sigma))
uc <- sapply(1:nS,function(k) qnorm(pgev(data[,k],paramc[k],scale[k],shape[k])))
rc <- -0.5/tau^2 * t(paramc-X%*%alpha-W)%*%(paramc-X%*%alpha-W) -
sum(sapply(1:nS,function(k) (1+1/shape[k])*log(1+shape[k]*(data[,k]-paramc[k])/scale[k]))) -
sum(sapply(1:nS,function(k) (1+shape[k]*(data[,k]-paramc[k])/scale[k])^(-1/shape[k]))) +
0.5*sum(uc*uc) - 0.5*sum(sapply(1:nT,function(k) uc[k,]%*%solve(C)%*%uc[k,]))
},
"scale" = {
param <- scale
paramc <- t(rmvnorm(1,param,Sigma))
uc <- sapply(1:nS,function(k) qnorm(pgev(data[,k],location[k],paramc[k],shape[k])))
rc <- -0.5/tau^2 * t(paramc-X%*%alpha-W)%*%(paramc-X%*%alpha-W) -
sum(sapply(1:nS,function(k) (1+1/shape[k])*log(1+shape[k]*(data[,k]-location[k])/paramc[k]))) -
sum(sapply(1:nS,function(k) (1+shape[k]*(data[,k]-location[k])/paramc[k])^(-1/shape[k]))) +
0.5*sum(uc*uc) - 0.5*sum(sapply(1:nT,function(k) uc[k,]%*%solve(C)%*%uc[k,]))
},
"shape" = {
param <- shape
paramc <- t(rmvnorm(1,param,Sigma))
uc <- sapply(1:nS,function(k) qnorm(pgev(data[,k],location[k],scale[k],paramc[k])))
rc <- -0.5/tau^2 * t(paramc-X%*%alpha-W)%*%(paramc-X%*%alpha-W) -
sum(sapply(1:nS,function(k) (1+1/paramc[k])*log(1+paramc[k]*(data[,k]-location[k])/scale[k]))) -
sum(sapply(1:nS,function(k) (1+paramc[k]*(data[,k]-location[k])/scale[k])^(-1/paramc[k]))) +
0.5*sum(uc*uc) - 0.5*sum(sapply(1:nT,function(k) uc[k,]%*%solve(C)%*%uc[k,]))
})
ri <- -0.5/tau^2 * t(param-X%*%alpha-W)%*%(param-X%*%alpha-W) -
sum(sapply(1:nS,function(k) (1+1/shape[k])*log(1+shape[k]*(data[,k]-location[k])/scale[k]))) -
sum(sapply(1:nS,function(k) (1+shape[k]*(data[,k]-location[k])/scale[k])^(-1/shape[k]))) +
0.5*sum(u*u) - 0.5*sum(sapply(1:nT,function(k) u[k,]%*%solve(C)%*%u[k,]))
r <- rc - ri
if(log(runif(1))<r) param <- paramc
param
}SpatialModel <- function(coord){
##Calcula las distancias entre observatorios
nS <- nrow(coord)
distance <- diag(0,nrow=nS,ncol=nS)
for(i in 1:(nS-1)){
for(j in (i+1):nS){
distance[i,j] <- sqrt((coord[i,1]-coord[j,1])^2 + (coord[i,2]-coord[j,2])^2)
distance[j,i] <- distance[i,j]
}
}
alt <- (coord[,3]-min(coord[,3]))/(max(coord[,3])-min(coord[,3]))
list(distance=distance,covariable=cbind(rep(1,nS),alt))
}HierBayesian <- function(data,coord,ndraw=1000,nthin=25,nburn=1/3*ndraw){
nS <- ncol(data)
nT <- nrow(data)
library(mvtnorm)
library(evd)
#variables espaciales
spatial <- SpatialModel(coord)
distance <- spatial$distance
X <- spatial$covariable
#copula
cop.params <- c(0.8,400)
#location
mu <- apply(data,2,mean)
alphamu <- c(50,-10)
betamu <- c(100,1000)
taumu <- 10
Sigmamu <- betamu[1]*exp(-distance/betamu[2])
Wmu <- t(rmvnorm(1,sigma=Sigmamu))
sdmu <- 2; SigmaAlphamu <- diag(10000,ncol=length(alphamu),nrow=length(alphamu))
#scale
sigma <- exp(rep(3,nS))
alphasigma <- matrix(10,nrow=1,ncol=1)
betasigma <- c(100,1000)
tausigma <- 10
Sigmasigma <- betasigma[1]*exp(-distance/betasigma[2])
Wsigma <- t(rmvnorm(1,sigma=Sigmasigma))
sdsigma <- 0.001; SigmaAlphasigma <- diag(100,ncol=length(alphasigma),nrow=length(alphasigma))
#shape
xi <- -0.25
sdxi <- 0.0001; a <- 5
nparam <- length(cop.params)+2*nS+length(alphamu)+length(alphasigma)+7
for(i in 1:nburn){
ui <- sapply(1:nS,function(k) qnorm(pgev(data[,k],mu[k],sigma[k],xi)))
##Cambio los parámetros de la copula (c0,c1)
cop.params <- CopEstimates(ui,distance,cop.params,nT)
C <- cop.params[1]*exp(-distance/cop.params[2])
diag(C) <- rep(1,nS)
##Cambio la localización con el modelo espacial
mu <- ParamEstimates("location",sdmu,mu,sigma,rep(xi,nS),X,alphamu,Wmu,taumu,data,ui,C)
alphamu <- AlphaEstimates(mu,Wmu,SigmaAlphamu,taumu,X)
Wmu <- WEstimates(mu,alphamu,Sigmamu,taumu,X)
betamu <- BetaEstimates(distance,mu,Wmu,Sigmamu,betamu)
Sigmamu <- betamu[1]*exp(-distance/betamu[2])
taumu <- TauEstimates(mu,Wmu,alphamu,X)
##Cambio la escala con el modelo espacial
sigma <- ParamEstimates("scale",sdsigma,mu,sigma,rep(xi,nS),X[,1],alphasigma,Wsigma,tausigma,data,ui,C)
alphasigma <- AlphaEstimates(sigma,Wsigma,SigmaAlphasigma,tausigma,X[,1])
Wsigma <- WEstimates(sigma,alphasigma,Sigmasigma,tausigma,X[,1])
betasigma <- BetaEstimates(distance,sigma,Wsigma,Sigmasigma,betasigma)
Sigmasigma <- betasigma[1]*exp(-distance/betasigma[2])
tausigma <- TauEstimates(sigma,Wsigma,alphasigma,X[,1])
##Cambio la forma sin modelo espacial (constante para todo el espacio)
ri <- -0.5/a^2*xi^2 -
sum(sapply(1:nS,function(k) (1+1/xi)*log(1+xi*(data[,k]-mu[k])/sigma[k]))) -
sum(sapply(1:nS,function(k) (1+xi*(data[,k]-mu[k])/sigma[k])^(-1/xi))) +
0.5*sum(ui*ui) - 0.5*sum(sapply(1:nT,function(k) ui[k,]%*%solve(C)%*%ui[k,]))
xic <- rnorm(1,xi,sdxi)
uc <- sapply(1:nS,function(k) qnorm(pgev(data[,k],mu[k],sigma[k],xic)))
rc <- -0.5/a^2*xic^2 -
sum(sapply(1:nS,function(k) (1+1/xic)*log(1+xic*(data[,k]-mu[k])/sigma[k]))) -
sum(sapply(1:nS,function(k) (1+xic*(data[,k]-mu[k])/sigma[k])^(-1/xic))) +
0.5*sum(uc*uc) - 0.5*sum(sapply(1:nT,function(k) uc[k,]%*%solve(C)%*%uc[k,]))
r <- rc - ri
if(log(runif(1))<r) xi <- xic
}
draw <- matrix(c(cop.params,
mu,alphamu,betamu,taumu,
log(sigma),log(alphasigma),betasigma,tausigma,
xi),
nrow=ndraw,ncol=nparam,byrow=TRUE)
for(i in 2:ndraw){
for(j in 1:nthin){
ui <- sapply(1:nS,function(k) qnorm(pgev(data[,k],mu[k],sigma[k],xi)))
##Cambio los parámetros de la copula (c0,c1)
cop.params <- CopEstimates(ui,distance,cop.params,nT)
C <- cop.params[1]*exp(-distance/cop.params[2])
diag(C) <- rep(1,nS)
##Cambio la localización con el modelo espacial
mu <- ParamEstimates("location",sdmu,mu,sigma,rep(xi,nS),X,alphamu,Wmu,taumu,data,ui,C)
alphamu <- AlphaEstimates(mu,Wmu,SigmaAlphamu,taumu,X)
Wmu <- WEstimates(mu,alphamu,Sigmamu,taumu,X)
betamu <- BetaEstimates(distance,mu,Wmu,Sigmamu,betamu)
Sigmamu <- betamu[1]*exp(-distance/betamu[2])
taumu <- TauEstimates(mu,Wmu,alphamu,X)
##Cambio la escala con el modelo espacial
sigma <- ParamEstimates("scale",sdsigma,mu,sigma,rep(xi,nS),X[,1],alphasigma,Wsigma,tausigma,data,ui,C)
alphasigma <- AlphaEstimates(sigma,Wsigma,SigmaAlphasigma,tausigma,X[,1])
Wsigma <- WEstimates(sigma,alphasigma,Sigmasigma,tausigma,X[,1])
betasigma <- BetaEstimates(distance,sigma,Wsigma,Sigmasigma,betasigma)
Sigmasigma <- betasigma[1]*exp(-distance/betasigma[2])
tausigma <- TauEstimates(sigma,Wsigma,alphasigma,X[,1])
##Cambio la forma sin modelo espacial (constante para todo el espacio)
ri <- -0.5/a^2*xi^2 -
sum(sapply(1:nS,function(k) (1+1/xi)*log(1+xi*(data[,k]-mu[k])/sigma[k]))) -
sum(sapply(1:nS,function(k) (1+xi*(data[,k]-mu[k])/sigma[k])^(-1/xi))) +
0.5*sum(ui*ui) - 0.5*sum(sapply(1:nT,function(k) ui[k,]%*%solve(C)%*%ui[k,]))
xic <- rnorm(1,xi,sdxi)
uc <- sapply(1:nS,function(k) qnorm(pgev(data[,k],mu[k],sigma[k],xic)))
rc <- -0.5/a^2*xic^2 -
sum(sapply(1:nS,function(k) (1+1/xic)*log(1+xic*(data[,k]-mu[k])/sigma[k]))) -
sum(sapply(1:nS,function(k) (1+xic*(data[,k]-mu[k])/sigma[k])^(-1/xic))) +
0.5*sum(uc*uc) - 0.5*sum(sapply(1:nT,function(k) uc[k,]%*%solve(C)%*%uc[k,]))
r <- rc - ri
if(log(runif(1))<r) xi <- xic
}
draw[i,] <- c(cop.params,
mu,alphamu,betamu,taumu,
log(sigma),log(alphasigma),betasigma,tausigma,
xi)
}
print(paste("Copula parameters: c0 = ",mean(draw[,1])," c1 = ",mean(draw[,2]),sep=" "))
print(paste("Spatial parameter of location: alpha0 = ",mean(draw[,nS+3])," alpha1 = ",mean(draw[,nS+4]),sep=" "))
print(paste("Spatial parameter of scale: alpha0 = ",mean(draw[,2*nS+8]),sep=" "))
for(k in 1:nS){
kk <- 2+k; kkk <- nS+k+7
print(paste("Observatory number ",k,sep=" "))
print(paste("location = ",mean(draw[,kk]),"scale = ",mean(draw[,kkk]),
"shape = ",mean(draw[,2*nS+12]),sep=" "))
}
list(MCMC=list(c0.copula=draw[,1],c1.copula=draw[,2],
location=draw[,3:(nS+2)],alpha0.location=draw[,nS+3],alpha1.location=draw[,nS+4],
beta0.location=draw[,nS+5],beta1.location=draw[,nS+6],tau.location=draw[,nS+7],
scale=draw[,(nS+8):(2*nS+7)],alpha0.scale=draw[,2*nS+8],beta0.scale=draw[,2*nS+9],
beta1.scale=draw[,2*nS+10],tau.scale=draw[,2*nS+11],shape=draw[,2*nS+12]))
}Agradecimientos
Ministerio de Ciencia e Innovación y Agencia Estatal de Investigación: proyecto PID2021-122209OB-C32.
European Union: proyecto PID2021-122209OB-C32.
