R/Densities.R

In ExtremalDep: Extremal Dependence Models

####################################################################################
### Authors: Boris Beranger and Simone Padoan        	 		                 ###
### 				   				                                             ###
###	Emails: borisberanger@gmail.com, simone.padoan@unibocconi.it				 ###
### 																			 ###
###	Institutions: Department of Decision Sciences, University Bocconi of Milan	 ###
### School of Mathematics and Statistics, University of New South Wales 		 ###
### 																			 ###
### File name: Densities.r	                 							         ###
### 																			 ###
### Description:                                  		              			 ###
### This file provides the Angular densities of extremal dependence models: 	 ###
### 1) Pairwise Beta (PB), Tilted Dirichlet (TD) and Husler-Reiss (HR)           ###
###     with only mass on the interior of the simplex)                           ###
### 2) Asymmetric logistic (AL), Extremal-t (ET) and extremal skew-t (EST)       ###
###     with mass on all subsets of the simplex)				                 ###
### It also provides the likelihood and log-likelihood functions				 ###
### 																			 ###
### Last change: 06/08/2019                         		  					 ###
### 																			 ###
####################################################################################

### Estimation and generation from angular density

angular <- function(data, model, n, dep, asy, alpha, beta, df, seed, k, nsim, plot=TRUE, nw=100){
  
  w <- seq(0.00001, .99999, length=nw)
  
  # Simulation of synthetic data
  if(missing(data)){
    if(missing(model) || missing(n) || missing(dep) || missing(seed) ){stop("model, n. dep and missing must be specified when data is not provided")}
    # warning("Data not provided and generated according to the parameters provided")
    
    if(any(model==c("log", "alog", "hr", "neglog", "aneglog", "bilog", "negbilog", "ct", "amix"))){
      set.seed(seed)
      data <- rbvevd(n=n, dep = dep, asy, alpha, beta, model = model, mar1 = c(1,1,1))
      Atrue <- abvevd(w, dep=dep, asy, alpha, beta, model=model) # True pickands dependence function
      htrue <- hbvevd(1-w, dep=dep, asy, alpha, beta, model=model, half=TRUE) # True angular density
      if(any(model==c("alog","aneglog"))){Htrue <- (1-asy)/2}
      if(model=="amix"){Htrue <- c(1-alpha-beta, 1-alpha-2*beta)/2}
    }else if(model=="Extremalt"){
      set.seed(seed)
      data <- rExtDep(n=n, model=model, par=c(dep, df), angular=FALSE, mar=c(1,1,1), num=5e+5)
      Atrue <- rep(NA,nw)
      for(i in 1:nw){Atrue[i] <- index.ExtDep(object="pickands", model="ET", par=c(dep,df), x=c(w[i],1-w[i]))}
      htrue <- dExtDep(x=cbind(w,1-w), method="Parametric", model=model, par=c(dep,df), angular=TRUE, c=0, log=FALSE, vectorial=TRUE)/2
      Htrue <- dExtDep(x=cbind(c(1,0),c(0,1)), method="Parametric", model=model, par=c(dep,df), angular=TRUE, c=0, log=FALSE, vectorial=TRUE)/2
    }
  }else{
    if(!is.matrix(data) || ncol(data)!=2){stop("data must be a matrix with 2 columns")}
    n <- nrow(data)
    model <- dep <- seed <- NULL
    Atrue <- htrue <- 0
    warning("True Pickands function and angular density not computed as data is provided and true model is unsure")
  }
  
  if(missing(k)){stop("k, the polynomial order must be specified")}
  if(missing(nsim)){stop("nsim, the number of of generation from the estimated angular density must be specified")}  
  
  # Compute the angles:
  wdata <- data[,1]/rowSums(data)
  
  # Estimate the pickands function:
  Aest <- beed(data, cbind(w, 1-w), 2, 'md', 'emp', k=k, plot=FALSE)
  beta <- Aest$beta
  
  # Compute the angular density in the interior
  hest <- sapply(1:nw, function(i, w, beta) dh(w[i], beta), w, beta)
  
  # Compute the angular measure
  Hest <- sapply(1:nw, function(i, w, beta) ph(w[i], beta), w, beta)
  
  # Compute the point masses on the corners
  p0 <- (1+k*(beta[2]-1))/2
  p1 <- (1-k*(1-beta[k]))/2
  
  # simulate from the semiparametric angular density
  wsim <- rh(nsim, beta)
  
  if(plot){
    par(mai = c(0.85, 0.85, 0.1, 0.1), mgp = c(2.8, 1, 0), cex.axis=2, cex.lab=2)
    hist(wsim[wsim!=0 & wsim!=1], freq=FALSE, ylim=c(0,3.5), xlim=c(0,1), xlab='w', ylab='h(w)', main="",col="gray")
    lines(w, hest, col=1, lty=2, lwd=2.5)
    if(is.vector(htrue)){ lines(w[2:99], htrue[2:99],col=2, lty=1, lwd=1.5)}
    points(1,sum(wsim==1)/nsim, pch=19, cex=2)
    points(0,sum(wsim==0)/nsim, pch=19, cex=2)
    if(any(model==c('alog','aneglog','amix','Extremalt'))){
      points(0, Htrue[1] , pch=21, cex=2,col=2, lwd=2)
      points(1, Htrue[2], pch=21, cex=2,col=2, lwd=2)
    }
  }
  
  return(list(model=model, n=n, dep=dep, data=data, Aest=Aest, hest=hest, Hest=Hest, p0=p0, p1=p1, wsim=wsim, Atrue=Atrue, htrue=htrue))
  
}