R/loglikelihood.R

In ExtremalDep: Extremal Dependence Models

#######################################################
### Authors: Giulia Marcon and Simone Padoan        ###
### Emails: giulia.marcon@phd.unibocconi.it,        ###
### simone.padoan@unibocconi.it                     ###
### Institution: Department of Decision Sciences,   ###
### University Bocconi of Milan                     ###
### File name: loglikelihood.r                      ###
### Description:                                    ###
### This file provides the log-likelihood function  ###
### corresponding to the bivariate max-stable       ###
### distribution, using the angular distribution    ###
### or the Pickands dependence function,            ###
### in Bernstein form.                              ###
### See eq. 3.16 in Marcon et al. (2016)            ###
### Last change: 15/08/2016                         ###
#######################################################

###################### START LOG-LIKELIHOOD ####################################

################################################################################
# INPUT:                                                                     ###
# coef is a vector of the eta coefficients                                   ###
# k is the polynomial order                                                  ###
# data is a list of:                                                         ###
#      z = 1/x + 1/y (Frechet scale)                                         ###
#      w = angular of data                                                   ###
#      r = radius of the data                                                ###
#      w2 = w^2                                                              ###
#      r2 = r^2                                                              ###
#      den = x^2, y^2                                                        ###
# pm is a vector of point masses at zero and one                             ###
# bpb, bpb1, bpb2 are matrices of the Bernstein polynomial basis             ###
# nsim is the number of the iterations of the chain                          ###
# if approx = TRUE, only coef, k and bp1 are needed.                         ###
################################################################################ 

llik <- function(coef, k, data = NULL, pm = NULL, bpb = NULL, bpb1, bpb2 = NULL, approx = FALSE){
  # z: 1/x + 1/y (Frechet scale)
  # w angular of data
  # r: radius of the data
  # w2 <- w^2
  # r2 <- r^2
  # den <- x^2, y^2
  
  if(approx){
    # compute the angular density:
    h <- k*c(bpb1 %*% diff(coef))
    llik <- log(h)
  }
  
  else{
    p0 <- pm$p0
    p1 <- pm$p1
    # derive the betas coefficients
    beta <- net(coef, from='H')$beta
  
    # compute the difference of the coefficients:
    beta1 <- diff(beta); beta2 <- diff(beta1)
  
    # compute the Pickands and its derivatives:
    A <- c(bpb %*% beta)
    A1 <- c((k+1) * (bpb1 %*% beta1))
    A2 <- c(k * (k+1) * (bpb2 %*% beta2))
  
    # GEV distribution
    lG <- -data$z*A
  
    C <- data$w*A1
    A1.s <- (A - C[,1])*(A + C[,2]) / data$den[,2]
    A2.s <- A2*data$w2[,1]/data$r
    llik <- lG + log(A1.s+A2.s) - log(data$den[,1])
  }
    return(sum(llik))
}

###################### END LOG-LIKELIHOOD ####################################


###################### START CONSTRAINED MAXIMUM LOG-LIKELIHOOD FITTING ####################################
###################### FOR MAXIMA AND THRESHOLD EXCEEDANCES             ####################################
constmle <- function(data, k, w, start=NULL, type="maxima", r0=NULL, q=NULL){
  ################################
  # START auxilary functions
  
  llik_pickands_min <- function(coef){
    # compute the difference of the coefficients:
    beta1 <- diff(coef); beta2 <- diff(beta1)
    
    # compute the Pickands and its derivatives:
    A <- c(bpb %*% coef)
    A1 <- c(k * (bpb1 %*% beta1))
    A2 <- c(k * (k-1) * (bpb2 %*% beta2))
    
    # GEV distribution
    lG <- -pseudo$z*A
    
    C <- pseudo$w*A1
    A1.s <- (A - C[,1])*(A + C[,2]) / pseudo$den[,2]
    A2.s <- A2*pseudo$w2[,1]/pseudo$r
    if(any((A1.s+A2.s)<0)){
      llik <- 1e10
    }else{
      llik <- lG + log(A1.s+A2.s) - log(pseudo$den[,1])
      llik <- -sum(llik)  
    }
    
    if(is.na(llik)) llik <- 1e10
    return(llik)
  }
  
  llik_h_beta <- function(coef){
    # define the eta coefficients
    #eta <- net(beta=coef,from='A')$eta
    eta <- 1/2 + k*diff(coef)/2
    # compute the angular density:
    h <- (k-1)*c(bpb2 %*% diff(eta))
    # check if it is the same to dh in CODE_Simulation...
    # is it faster?
    if(any(h<0)){
      llik <- 1e100
    }else{
      llik <- log(h)
      llik <- -sum(llik)      
    }
  
    if(is.na(llik)) llik <- 1e100
    return(llik)
  }
  
  confun <- function(x) return(constraint$r - constraint$R%*%x)
  
  # END auxilary functions
  ################################
  kp <- k+1
  km <- k-1
  # define the linear constraints
  constraint <- constraints(k)
  # define the setup for the estimation
  pseudo <- list(z=rowSums(1/data), r=rowSums(data), w=data/rowSums(data), 
                 r2=(rowSums(data))^2, w2=(data/rowSums(data))^2, den=data^2)
  # define options
  #opts <- list(algorithm="NLOPT_LN_NELDERMEAD", xtol_rel=1e-25, maxeval=100000)
  opts <- list(algorithm="NLOPT_LN_COBYLA", xtol_rel=1e-25, maxeval=100000)
  #opts <- list(algorithm="NLOPT_LN_NEWUOA_BOUND", xtol_rel=1e-25, maxeval=100000)
  #opts <- list(algorithm="NLOPT_LN_SBPLX", xtol_rel=1e-08)
  # Compute starting values
  #if(is.null(start)){
  #  eta0 <- prior_eta_sampler(km)$eta
  #  beta0 <- net(eta=eta0, from='H')$beta}
  #else beta0 <- start
  if(is.null(start)) beta0 <- rcoef(km)$beta
  else beta0 <- start

  # maximize the constraint log-likelihood function
  if(type=="maxima"){
    bpb <- bpb1 <- bpb2 <- NULL
    bpb <- bp2d(pseudo$w, k)
    bpb1 <- bp2d(pseudo$w, km)
    bpb2 <- bp2d(pseudo$w, km-1)
    fit <- nloptr(beta0, eval_f=llik_pickands_min, eval_g_ineq=confun, 
                  lb=rep(0,kp), ub=rep(1,kp), opts=opts)
  }
  if(type=="rawdata"){
    bpb2 <- NULL 
    bpb2 <- bp2d(pseudo$w[pseudo$r>r0, ], km-1)
    fit <- list()
    fvalue <- c(Inf, numeric(14))
    fit[[1]] <- nloptr(beta0, eval_f=llik_h_beta, eval_g_ineq=confun, 
                       lb=rep(0,kp), ub=rep(1,kp), opts=opts)
    for(i in 2:15){
      fit[[i]]<- nloptr(fit[[i-1]]$solution[kp:1], eval_f=llik_h_beta, eval_g_ineq=confun, 
                        lb=rep(0,kp), ub=rep(1,kp), opts=opts)
      fvalue[i] <- fit[[i]]$objective}
    
    ind <- c(which(fvalue==min(fvalue)))[1]
    fit <- fit[[ind]]
  }
  beta <- fit$solution[kp:1]
  A <- beed(data, cbind(w, 1-w), 2, 'md', 'emp', k, beta=beta, plot=FALSE)
  
  return(list(beta=beta, A=A$A, status=fit$status, start=beta0, iterations=fit$iterations, message=fit$message))
}

###################### FOR MAXIMA AND THRESHOLD EXCEEDANCES             ####################################
###################### END CONSTRAINED MAXIMUM LOG-LIKELIHOOD FITTING ####################################