xparetotvtcov<-function(lon,lat,array3d,span,p=0.9,index,covariates=NULL,sigl=NULL,shl=NULL,siglink=identity,shlink=identity,nonmissing=0.5) {

# Fit Generalized Pareto distribution with time-varying threshold at each grid 
# point for a given three-dimensional array of montly data with first two space 
# dimensions (e.g. longitude and latitude) and third time dimension. Allows
# linear modelling of the paramters.
#
# Description:
#
#      Returns fitted parameters. 
#
# Usage:
#
#      xparetotvtcov(lon,lat,array3d,span,p,index,covariates,sigl,shl,siglink,shlink,nonmissing)
#
# Input:
#
#       lon: vector with p longitude values
#
#       lat: vector with q latitude values 
# 
#   array3d: a three-dimensional array of monthly data with p longitude points
#            and q latitude points as the first two dimensions and n as the
#            third time dimension 
#
#      span: fraction of the total number of points 'n' to be used to compute
#            the long term mean
#
#         p: a value between 0 and 1 that informs the percentage of point
#            that will be left below the time-varying threshold 
# 
#     index: index vector of length n composed of logicals (TRUE and FALSE) 
#            defining the months to be considered when computing the
#            time-varying threshold
#
#covariates: Matrix of covariates for generalized linear modelling of
#            the parameters. The number of rows should be 'n' (i.e. the same as
#	     the time dimension of 'array3d'.
#
# sigl, shl: Numeric vectors of integers, giving the columns of 'ydat'
#            that contain covariates for generalized linear modelling of
#            the scale and shape parameters repectively .
#
# siglink, shlink: Inverse link functions for generalized linear
#          modelling of the scale and shape parameters repectively.
#
#nonmissing: Only grid points with fraction given by 'nonmissing' 
#            (between 0 and 1) of non-missing values are used to estimate the
#            Generalized Pareto distribution parameters. Default is 0.5.
#
# Output:
#
#     $output: n' x p x q array with estimates parameters, where n' is
#              the total number of estimated parameters
#
# Authors:
#
#      Caio Coelho <c.a.d.s.coelho@reading.ac.uk> 28 Feb 2006
#      Chris Ferro <c.a.t.ferro@reading.ac.uk>


# reshape three-dimensional array into a data matrix with
# time as first dimension and space as sencond dimension
data <- reshapefield(lon,lat,array3d)$out

# compute percentage of non-missing values at each grid point
aux <- apply(data[index,],2,function(x)sum(!is.na(x))/(length(x)))

# identify grid points with less than 50% missing values
indexgrid <- (1:length(aux))[aux >= nonmissing]

outaux<-matrix(ncol=dim(data)[2],nrow=(length(sigl)+length(shl)+2))

# Estimate Generalized Pareto distribution parameters (scale and shape) for
# exceedances above the threshold p
estim.par <- function(y,p,index,covariates,sigl,shl,siglink,shlink) {

threshold<-tvt(y,span,index,1-p)$threshold

yy<-y[index]
covs<-as.matrix(covariates[index,])

index1 <- !(is.na(yy) | apply(t(apply(covs,1,is.na)),1,any))

params<-mygpd.fit(yy[index1],threshold[index1],ydat=as.matrix(covs[index1,]),sigl=sigl,shl=shl,siglink=siglink,shlink=shlink,show=FALSE)$mle

list(params=params)

}

out<-apply(data[,indexgrid],2,estim.par,p,index,covariates,sigl,shl,siglink,shlink)

for(j in 1:(length(sigl)+length(shl)+2)){
for (i in 1:length(indexgrid)){
outaux[j,indexgrid[i]]<-out[[i]]$params[j]
}
}

output<-aperm(reshapefield(lon,lat,outaux)$out,c(3,1,2))

invisible(list(output=output))

}