# R/EstimationBR.r In spatialTailDep: Estimation of spatial tail dependence models

#### Defines functions MestimatorBRAsymVarBRpsidotBRintBRcase5intBRcase4intBRcase3intBRcase2intBRcase1I3I1v4I1v3I1v2I1I3BRI2BRIBRellBR4dv4ellBR4dv3ellBR4dv2tailBRcorrfuncgfuncellBR4dellBRtominimizeBRtailBRIntparsBRtoTauparsTauToBR

```#' @include gridUtilities.r
#' @include tailIntEmp.r
NULL

parsTauToBR<-function(tau){
wortel<-sqrt(1 + (4*(tau[3]^2))/((tau[2]-tau[1])^2))
rat<-(tau[2]+tau[1])/(tau[1]-tau[2])
beta<-0.5*atan((2*tau[3])/(tau[2]-tau[1]))
if(beta >=0){
cc<-sqrt((rat-wortel)/(rat+wortel))
rho<-sqrt((cc^2 - 1)/(wortel*(tau[2]-tau[1])))
} else{
beta<-0.5*(atan((2*tau[3])/(tau[2]-tau[1]))+pi)
cc<-1/sqrt((rat-wortel)/(rat+wortel))
rho<-sqrt((cc^2 - 1)/(wortel*(tau[1]-tau[2])))
}
return(c(rho,beta,cc))
}

parsBRtoTau<-function(pars){
tau11<-(1/(2*(pars[1]^2)))*(1 + pars[3]^2 - (pars[3]^2 - 1)*cos(2*pars[2]))
tau22<-(1/(2*(pars[1]^2)))*(1 + pars[3]^2 + (pars[3]^2 - 1)*cos(2*pars[2]))
tau12<-(((pars[3]^2 - 1)/(2*(pars[1]^2)))*sin(2*pars[2]))
return(c(tau11,tau22,tau12))
}

tailBRInt <- function(loc, tau, alpha){
locdiff<-loc[3:4]-loc[1:2]
a <- sqrt(2)*((t(locdiff) %*% tau %*% locdiff)^(alpha/4))
logres<-(a^2 + pnorm((-3*a)/2,log.p=TRUE) - log(3))
return(as.vector(pnorm(a/2) + exp(logres)))
}

tominimizeBR<-function(pars,totlist,pairs,w){
tau <- rbind(c(pars[2], pars[4]), c(pars[4], pars[3]))
res<-apply(pairs,1,tailBRInt,tau=tau,alpha=pars[1]) - totlist
return(t(res) %*% w %*% res)
}

ellBR<-function(t,loc,tau,alpha){
if((loc[1]==0) && (loc[2]==0)){
return(pmax(t[1],t[2]))
} else{
locdiff<-loc[3:4]-loc[1:2]
a <- sqrt(2)*((t(loc) %*% tau %*% loc)^(alpha/4))
res<-.C("ellsmith2",as.double(t),as.double(a),result=double(1),PACKAGE="spatialTailDep")\$result
return(res)
}
}

ellBR4d<-function(t,rlist,gamlist,d){
result<-vector(length=d)
for(i in 1:d){
const<-(gamlist[[i]]/2 + log(t[i]/t[-i])/gamlist[[i]])
result[i]<-t[i]*pmvnorm(lower=rep(-Inf,d-1),upper=const,corr=rlist[[i]],algorithm=TVPACK)[1]
}
return(sum(result))
}
gfunc<-function(tau,alpha,loc){
a2 <- ((t(loc) %*% tau %*% loc)^(alpha/2))
return(a2)
}

corrfunc<-function(tau,alpha,loc,d,i){
result<-matrix(1,nrow=(d-1),ncol=(d-1))
for(j in 1:(d-2)){
for(k in (j+1):(d-1)){
lo<-loc[,-i]
result[j,k]<-result[k,j]<-((gfunc(tau,alpha,loc[,i]-lo[,j]) + gfunc(tau,alpha,loc[,i]-lo[,k]) -
gfunc(tau,alpha,lo[,j]-lo[,k]))/(2*sqrt(gfunc(tau,alpha,loc[,i]-lo[,j])*gfunc(tau,alpha,loc[,i]-lo[,k]))))
}
}
return(result)
}

tailBR <- function(t, loc, tau, alpha, d){
if(d==2){
return(ellBR(t,loc[,2]-loc[,1],tau,alpha))
} else{
rlist<-gamlist<-vector('list',length=d)
for(i in 1:d){
rlist[[i]]<-corrfunc(tau,alpha,loc,d,i)
gamlist[[i]]<-sapply(c(1:d)[-i], function(j) sqrt(2*gfunc(tau,alpha,loc[,i]-loc[,j])))
}
return(ellBR4d(t,rlist,gamlist,d))
}
}

ellBR4dv2<-function(t,rlist,gamlist,d,ind){
return(2*t[ind]*ellBR4d(t,rlist,gamlist,d))
}
ellBR4dv3<-function(t,coord,tau,alpha,ind){
return(ellBR(c(pmax(t[ind[1]],t[ind[2]]),t[ind[3]]),coord[,2]-coord[,1],tau,alpha))
}
ellBR4dv4<-function(t,coord,tau,alpha){
return(4*t[1]*t[2]*ellBR(t,coord[,2]-coord[,1],tau,alpha))
}

IBR<-function(a){ #  This function can be implemented in C as well.
return(pnorm(a/2) + exp(a*a)*(1/3)*pnorm(-(3/2)*a))
}
I2BR<-function(x,a){  #  This function can be implemented in C as well.
return(0.5*pnorm(a/2 - log(x)/a) + x*pnorm(a/2 + log(x)/a) + 0.5*x*x*exp(a*a)*pnorm(-(3/2)*a - log(x)/a))
}
I3BR<-function(x,a){
res<-.C("I3smith2",as.double(x),as.double(a),result=double(1),PACKAGE="spatialTailDep")\$result
return(res)
}

I1<-function(u,a,rlist,gamlist,d,ind){ return(I3BR(u[ind],a)*ellBR4d(u,rlist,gamlist,d))}
I1v2<-function(u,a,tau,alpha,coordt,ind,ind2){ return(I3BR(u[ind],a)*ellBR4dv3(u,coordt,tau,alpha,ind2))}
I1v3<-function(u,a,rlist,gamlist,d){return(I3BR(u[1],a)*ellBR4d(c(u[1],u[3],u[2]),rlist,gamlist,d))}
I1v4<-function(u,a,rlist,gamlist,d){return(I3BR(u[1],a)*ellBR4d(c(u[2],u[1],u[3]),rlist,gamlist,d))}
I3<-function(u,coord,tau,alpha,coordt){
loc1<-coord[,2]-coord[,1]
loc2<-coord[,4]-coord[,3]
anr1 <- sqrt(2)*((t(loc1) %*% tau %*% loc1)^(alpha/4))
anr2 <- sqrt(2)*((t(loc2) %*% tau %*% loc2)^(alpha/4))
return(I3BR(u[1],anr1)*I3BR(u[2],anr2)*ellBR(u,coordt[,2]-coordt[,1],tau,alpha))}

####### case 1: locations s,t,u,v all different
intBRcase1<-function(coord,tau,alpha,Tol){ #coord=(s,t,u,v) where s,t,u,v, are 2-dim columnvectors
loc1<-coord[,2]-coord[,1]
loc2<-coord[,4]-coord[,3]
anr1 <- sqrt(2)*((t(loc1) %*% tau %*% loc1)^(alpha/4))
anr2 <- sqrt(2)*((t(loc2) %*% tau %*% loc2)^(alpha/4))
########### T1
d<-4
rlist<-gamlist<-vector('list',length=d)
for(i in 1:d){
rlist[[i]]<-corrfunc(tau,alpha,coord,d,i)
gamlist[[i]]<-sapply(c(1:d)[-i], function(j) sqrt(2*gfunc(tau,alpha,coord[,i]-coord[,j])))
}
############# T2
d<-3
rlistpt1<-rlistpt2<-gamlistpt1<-gamlistpt2<-vector('list',length=d)
coordt1<-coord[,-4]
coordt2<-coord[,-3]
for(i in 1:d){
rlistpt1[[i]]<-corrfunc(tau,alpha,coordt1,d,i)
gamlistpt1[[i]]<-sapply(c(1:d)[-i], function(j) sqrt(2*gfunc(tau,alpha,coordt1[,i]-coordt1[,j])))
rlistpt2[[i]]<-corrfunc(tau,alpha,coordt2,d,i)
gamlistpt2[[i]]<-sapply(c(1:d)[-i], function(j) sqrt(2*gfunc(tau,alpha,coordt2[,i]-coordt2[,j])))
}
T2<-(2*IBR(anr1)*(I2BR(1,anr2) -0.5) + 2*I2BR(1,anr2) - 2*IBR(anr2) -
############# T3
d<-3
rlistpt1<-rlistpt2<-gamlistpt1<-gamlistpt2<-vector('list',length=3)
coordt1<-coord[,-2]
coordt2<-coord[,-1]
for(i in 1:d){
rlistpt1[[i]]<-corrfunc(tau,alpha,coordt1,d,i)
gamlistpt1[[i]]<-sapply(c(1:d)[-i], function(j) sqrt(2*gfunc(tau,alpha,coordt1[,i]-coordt1[,j])))
rlistpt2[[i]]<-corrfunc(tau,alpha,coordt2,d,i)
gamlistpt2[[i]]<-sapply(c(1:d)[-i], function(j) sqrt(2*gfunc(tau,alpha,coordt2[,i]-coordt2[,j])))
}
T3<-(2*IBR(anr2)*(I2BR(1,anr1) -0.5) + 2*I2BR(1,anr1) - 2*IBR(anr1) -
################ T4
T4<-((IBR(anr1) - I2BR(1,anr1))*(1 - 2*I2BR(1,anr2)) + (IBR(anr2) -  I2BR(1,anr2))*(1 - 2*I2BR(1,anr1)) -
################# T5
T5<-((IBR(anr1) - I2BR(1,anr1))*(1 - 2*I2BR(1,anr2)) + (IBR(anr2) -  I2BR(1,anr2))*(1 - 2*I2BR(1,anr1)) -
return(T1 - T2 - T3 + T4 + T5)
}

## case 2: t = u
intBRcase2<-function(coord,tau,alpha,Tol){ #coord=(s,t,u,v) where s,t,u,v, are 2-dim columnvectors
loc1<-coord[,2]-coord[,1]
loc2<-coord[,4]-coord[,3]
anr1 <- sqrt(2)*((t(loc1) %*% tau %*% loc1)^(alpha/4))
anr2 <- sqrt(2)*((t(loc2) %*% tau %*% loc2)^(alpha/4))
########### T1
d<-3
rlist<-gamlist<-vector('list',length=d)
coordt<-coord[,-3]
for(i in 1:d){
rlist[[i]]<-corrfunc(tau,alpha,coordt,d,i)
gamlist[[i]]<-sapply(c(1:d)[-i], function(j) sqrt(2*gfunc(tau,alpha,coordt[,i]-coordt[,j])))
}
T1<-(IBR(anr1) + IBR(anr2) -
############# T2
T2<-(2*IBR(anr1)*(I2BR(1,anr2) -0.5) + 2*I2BR(1,anr2) - 2*IBR(anr2) -
############# T3
T3<-(2*IBR(anr2)*(I2BR(1,anr1) -0.5) + 2*I2BR(1,anr1) - 2*IBR(anr1) -
################ T4
T4<-((IBR(anr1) - I2BR(1,anr1))*(1 - 2*I2BR(1,anr2)) + (IBR(anr2) -  I2BR(1,anr2))*(1 - 2*I2BR(1,anr1)) -
################# T5
T5<-((IBR(anr1) - I2BR(1,anr1))*(1 - 2*I2BR(1,anr2)) + (IBR(anr2) -  I2BR(1,anr2))*(1 - 2*I2BR(1,anr1)) -
return(T1 - T2 - T3 + T4 + T5)
}

## case 3: s = u, t neq v
intBRcase3<-function(coord,tau,alpha,Tol){ #coord=(s,t,u,v) where s,t,u,v, are 2-dim columnvectors
loc1<-coord[,2]-coord[,1]
loc2<-coord[,4]-coord[,3]
anr1 <- sqrt(2)*((t(loc1) %*% tau %*% loc1)^(alpha/4))
anr2 <- sqrt(2)*((t(loc2) %*% tau %*% loc2)^(alpha/4))
########### T1
d<-3
rlist<-gamlist<-vector('list',length=d)
coordt<-coord[,-3]
for(i in 1:d){
rlist[[i]]<-corrfunc(tau,alpha,coordt,d,i)
gamlist[[i]]<-sapply(c(1:d)[-i], function(j) sqrt(2*gfunc(tau,alpha,coordt[,i]-coordt[,j])))
}
############# T2
T2<-(2*IBR(anr1)*(I2BR(1,anr2) - 0.5) + 2*I2BR(1,anr2) - 2*IBR(anr2) -
############# T3
T3<-(2*IBR(anr2)*(I2BR(1,anr1) -0.5) + 2*I2BR(1,anr1) - 2*IBR(anr1) -
################ T4
T4<-((IBR(anr1) - I2BR(1,anr1))*(1 - 2*I2BR(1,anr2)) + (IBR(anr2) -  I2BR(1,anr2))*(1 - 2*I2BR(1,anr1)) -
################# T5
T5<-((IBR(anr1) - I2BR(1,anr1))*(1 - 2*I2BR(1,anr2)) + (IBR(anr2) -  I2BR(1,anr2))*(1 - 2*I2BR(1,anr1)) -
return(T1 - T2 - T3 + T4 + T5)
}

## case 4: s neq u, t = v
intBRcase4<-function(coord,tau,alpha,Tol){ #coord=(s,t,u,v) where s,t,u,v, are 2-dim columnvectors
loc1<-coord[,2]-coord[,1]
loc2<-coord[,4]-coord[,3]
anr1 <- sqrt(2)*((t(loc1) %*% tau %*% loc1)^(alpha/4))
anr2 <- sqrt(2)*((t(loc2) %*% tau %*% loc2)^(alpha/4))
########### T1
d<-3
rlist<-gamlist<-vector('list',length=d)
coordt<-coord[,-4]
for(i in 1:d){
rlist[[i]]<-corrfunc(tau,alpha,coordt,d,i)
gamlist[[i]]<-sapply(c(1:d)[-i], function(j) sqrt(2*gfunc(tau,alpha,coordt[,i]-coordt[,j])))
}
############# T2
T2<-(2*IBR(anr1)*(I2BR(1,anr2) - 0.5) + 2*I2BR(1,anr2) - 2*IBR(anr2) -
############# T3
T3<-(2*IBR(anr2)*(I2BR(1,anr1) - 0.5) + 2*I2BR(1,anr1) - 2*IBR(anr1) -
################ T4
T4<-((IBR(anr1) - I2BR(1,anr1))*(1 - 2*I2BR(1,anr2)) + (IBR(anr2) -  I2BR(1,anr2))*(1 - 2*I2BR(1,anr1)) -
################# T5
T5<-((IBR(anr1) - I2BR(1,anr1))*(1 - 2*I2BR(1,anr2)) + (IBR(anr2) -  I2BR(1,anr2))*(1 - 2*I2BR(1,anr1)) -
return(T1 - T2 - T3 + T4 + T5)
}

## case 5: s = u, t = v
intBRcase5<-function(coord,tau,alpha,Tol){ #coord=(s,t,u,v) where s,t,u,v, are 2-dim columnvectors
loc<-coord[,2]-coord[,1]
a <- sqrt(2)*((t(loc) %*% tau %*% loc)^(alpha/4))
########### T1
coordt<-coord[,1:2]
############# T2
T2<-(2*IBR(a)*(I2BR(1,a) - 0.5) + 2*I2BR(1,a) - 2*IBR(a) -
coordt=coordt,ind=3,ind2=c(1,3,2),tol=Tol)\$integral -
coordt=coordt,ind=3,ind2=c(2,3,1),tol=Tol)\$integral)
############# T3
T3<-(2*IBR(a)*(I2BR(1,a) - 0.5) + 2*I2BR(1,a) - 2*IBR(a) -
coordt=coordt,ind=1,ind2=c(1,2,3),tol=Tol)\$integral -
coordt=coordt,ind=1,ind2=c(1,3,2),tol=Tol)\$integral)
################ T4
T4<-((IBR(a) - I2BR(1,a))*(1 - 2*I2BR(1,a)) + (IBR(a) -  I2BR(1,a))*(1 - 2*I2BR(1,a)) -
################# T5
T5<-((IBR(a) - I2BR(1,a))*(1 - 2*I2BR(1,a)) + (IBR(a) -  I2BR(1,a))*(1 - 2*I2BR(1,a)) -
return(T1 - T2 - T3 + T4 + T5)
}

psidotBR<-function(tau,alpha,pairs){ #nrow(pairs) = number of pairs, ncol(pairs)=4
npairs<-nrow(pairs)
loc<-pairs[,3:4]-pairs[,1:2]
psi<-matrix(,nrow=npairs,ncol=4)
deriv<-vector(length=4)
for(i in 1:npairs){
a <- sqrt(2)*((t(loc[i,]) %*% tau %*% loc[i,])^(alpha/4))
temp<-sqrt(2)*(alpha/4)*((t(loc[i,]) %*% tau %*% loc[i,])^(alpha/4 - 1))
deriv[1]<-(a*log((t(loc[i,]) %*% tau %*% loc[i,])^(1/4)))
deriv[2]<-temp*(loc[i,1]^2)
deriv[3]<-temp*(loc[i,2]^2)
deriv[4]<-temp*(2*loc[i,1]*loc[i,2])
psi[i,]<-deriv*((1/2)*dnorm(a/2) - (1/2)*exp(a^2)*dnorm(-(3/2)*a) + (2/3)*a*exp(a^2)*pnorm(-(3/2)*a))
}
return(psi)
}

AsymVarBR<-function(pairs,tau,alpha,Tol){
npairs<-nrow(pairs)
BRMat<-matrix(,nrow=npairs,ncol=npairs)
for(i in 1:npairs){
for(j in i:npairs){
temp<-matrix(c(pairs[i,],pairs[j,]),nrow=2)
if(all(temp[,1]==temp[,4]) && (! all(temp[,2]==temp[,3]))){
BRMat[i,j]<-intBRcase2(cbind(temp[,3:4],temp[,1:2]),tau,alpha,Tol=Tol)
}else if(all(temp[,1]==temp[,3]) && all(temp[,2]==temp[,4])){
BRMat[i,j]<-intBRcase5(temp,tau,alpha,Tol=Tol)
}else if(all(temp[,2]==temp[,4]) && (! all(temp[,1]==temp[,3]))){
BRMat[i,j]<-intBRcase4(temp,tau,alpha,Tol=Tol)
}else if(all(temp[,1]==temp[,3]) && (! all(temp[,2]==temp[,4]))){
BRMat[i,j]<-intBRcase3(temp,tau,alpha,Tol=Tol)
}else if(all(temp[,2]==temp[,3]) && (! all(temp[,1]==temp[,4]))){
BRMat[i,j]<-intBRcase2(temp,tau,alpha,Tol=Tol)
}else if((! all(temp[,1]==temp[,3])) && (! all(temp[,2]==temp[,4]))){
BRMat[i,j]<-intBRcase1(temp,tau,alpha,Tol=Tol)
}
BRMat[j,i]<-BRMat[i,j]
}
}
return(BRMat)
}

MestimatorBR <- function(x, locations, pairIndices, k,Tol, startingValue, Omega,iterate,covMat) {
n <- dim(x)[1]
if (k < 1 || k > (n-1)) {
warning("k must be between 1 and n")
}
npairs <- nrow(pairIndices)
pairs <- pairCoordinates(locations, pairIndices)
ranks <- apply(x, 2, rank)
totL <- numeric(npairs)
for (l in 1:npairs) {
pairsOfRanks <- cbind(ranks[, pairIndices[l, 1]],ranks[, pairIndices[l, 2]])
totL[l] <- tailIntEmp(pairsOfRanks, n, k)
}
tau<-parsBRtoTau(startingValue[2:4])
alpha<-startingValue[1]
theta <- thetaPilot <- optim(par = c(alpha,tau), fn = tominimizeBR, method="BFGS", totlist = totL,
pairs = pairs,w = Omega,control=list(maxit=10000,reltol=1e-12))\$par
if (iterate) {
Omega <- solve(AsymVarBR(pairs, rbind(c(theta[2],theta[4]),c(theta[4],theta[3])), theta[1], Tol))
theta <- optim(par = thetaPilot,fn = tominimizeBR,method="BFGS",totlist = totL,
pairs = pairs,w = Omega,control=list(maxit=10000,reltol=1e-12))\$par
}
covMatrix<-NULL
if(covMat){
taum<-rbind(c(theta[2],theta[4]),c(theta[4],theta[3]))
psidot <- psidotBR(taum,theta[1], pairs)
temp <- solve(t(psidot) %*% Omega %*% psidot) %*% t(psidot)
GammaMat <- AsymVarBR(pairs, taum, theta[1], Tol)
covMatrix <- (temp %*% Omega %*% GammaMat %*% Omega %*% t(temp)) / k
}
return(list(theta = c(theta[1],parsTauToBR(theta[2:4])),
theta_pilot = c(thetaPilot[1],parsTauToBR(thetaPilot[2:4])),
covMatrix = covMatrix,
Omega = Omega))
}
```

## Try the spatialTailDep package in your browser

Any scripts or data that you put into this service are public.

spatialTailDep documentation built on May 29, 2017, 5:53 p.m.