Bivariate Distributions with twopiece and double twopiece marginals

Bivariate Gaussian Copula

The 2-dimensional Gaussian Copula is a distribution on the square \([0,1]\times [0,1]\). For a given correlation matrix \(R = \begin{pmatrix} 1 & \rho \\ \rho & 1\end{pmatrix}\), the Gaussian copula cumulative distribution function (cdf) \[C_{R}(u_1,u_2)=\Phi _{R}\left(\Phi ^{-1}(u_{1}),\Phi ^{-1}(u_{2})\right),\] where \(\Phi_R\) is the bivariate normal cdf with covariance matrix \(R\) (equal to the correlation matrix), and \(\Phi^{-1}\) is the quantile normal function (the inverse of the standard normal cdf). Simulating from this copula is simple as explained in Simulating from a Bivariate Gaussian Copula.

twopiece and double twopiece marginals

The bivariate Gaussian copula can be used to produce bivariate distributions with specific marginals. In many contexts, it is desirable to use flexible marginals that can capture heavy tails and asymmetry (Rubio and Steel 2015). One option to do so are double two piece (DTP, see DTP R package) distributions. An appealing member of this family is the Double two-piece Sinh-Arcsinh distribution (DTPSAS), which can capture asymmetry as well as heavier and lighter tails than those of the normal distribution (allowing for different tail behaviour in each direction).

Of course, other flexible marginals can be used instead. One option is the family of twopiece distributions, which contain fewer parameters but they have the same tail behaviour in each direction. These distributions are implemented in the R package twopiece. An appealing member of this family is the twopiece Sinh-ArcSinh distribution (TPSAS), which can be implemented with the twopiece R package as well.

The following R code shows how to fit a bivariate distribution with DTPSAS and TPSAS marginals and the Gaussian copula in the context of modelling log-returns.

Note that other copulas and marginals could be used as well (see for example: t-copula with t and two piece t marginals).

Modelling bivariate log-returns

Some required functions

#********************************************************************************
# Baseline functions
#********************************************************************************


#######################################################################
# Symmetric sinh-arcsinh (SAS) distribution
#######################################################################

# baseline symmetric SAS pdf
dsas0 <- function(x,delta,log=FALSE){
  logPDF <- dnorm(sinh(delta*asinh(x)),log=T) + log(delta) + log(cosh(delta*asinh(x))) -0.5*log(1+x^2)
  ifelse( is.numeric(logPDF),ifelse( log, return(logPDF), return(exp(logPDF)) ), logPDF )
}

# baseline symmetric SAS cdf
psas0 <- function(x,delta,log.p=FALSE){
  logCDF <- pnorm(sinh(delta*asinh(x)),log.p=T)
  ifelse( is.numeric(logCDF),ifelse( log.p, return(logCDF), return(exp(logCDF)) ), logCDF )
}


#####################################################
# Other functions
#####################################################

# logit function
logit <- Vectorize(function(p){
  val <- log(p) -  log(1-p)
  return(as.vector(val))
})

# expit functon
expit <- Vectorize(function(x){
  a <- exp(x)
  return(as.vector( a/(a + 1) ))          
})

Maximum Likelihood Estimation

# Required packages
library(mnormt)
library(mvtnorm)
library(ghyp)

## Loading required package: numDeriv

## Loading required package: MASS

library(copula)
library(fMultivar)

## Loading required package: timeDate

## Loading required package: timeSeries

## Loading required package: fBasics

# library(devtools)
# github_install("FJRubio67/twopiece")
library(twopiece)
# github_install("FJRubio67/DTP")
library(DTP)


# The SMI stocks data:  "SMI"  vs    "Swiss.Re"
data(smi.stocks)
Y = as.matrix(smi.stocks[ , c(1, 6)])
colnames(smi.stocks[ , c(1, 6)])

## [1] "SMI"      "Swiss.Re"

plot(Y)

##############################################################################################
# bivariate Gaussian-copula with DTP SAS marginals
##############################################################################################

# log likelihood function (reparameterised)
loglikDTP = function(par){
  # Reparameterisation
  mu = par[1:2]; sigma = exp(par[3:4]); 
  epsilon = 2*expit(par[5:6])-1; delta1 = exp(par[7:8]); delta2 <- exp(par[9:10])
  rho = 2*expit(par[11])-1; 
  # Ingredients of the log likelihood
  norm.cop = normalCopula(rho, dim = 2)
  probs =cbind(pdtp(Y[,1],mu[1],sigma[1],epsilon[1],delta1[1],delta1[2],param="eps",F=psas0, f=dsas0),
               pdtp(Y[,2],mu[2],sigma[2],epsilon[2],delta2[1],delta2[2],param="eps",F=psas0, f=dsas0 ))
  val.cop = sum(dCopula(probs, copula = norm.cop, log = TRUE))
  val.marg1 = sum( ddtp(Y[,1],mu[1],sigma[1],epsilon[1],delta1[1],delta1[2],param="eps",f=dsas0,log=T) )
  val.marg2 = sum( ddtp(Y[,2],mu[2],sigma[2],epsilon[2],delta2[1],delta2[2],param="eps",f=dsas0,log=T) )
  
  # output
  out <- -val.cop - val.marg1 - val.marg2
  
  return(out)
}


start2 = c(0,0,log(0.01),log(0.01),0,0,-0.5,-0.5,-0.5,-0.5,0) 

# Optimisation step
OPT.DTP = optim(start2,loglikDTP,control=list(maxit=10000))
OPT.DTP

## $par
##  [1]  0.001318143 -0.001491822 -5.297693393 -5.026813005  0.126564100
##  [6] -0.310798522 -0.476491305 -0.464640560 -0.657931559 -0.513102357
## [11]  1.756384654
## 
## $value
## [1] -10746.84
## 
## $counts
## function gradient 
##     3526       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

MLE.DTP <- c(OPT.DTP$par[1:2], exp(OPT.DTP$par[3:4]), 2*expit(OPT.DTP$par[5:6])-1, exp(OPT.DTP$par[7:10]),  2*expit(OPT.DTP$par[11])-1  )
MLE.DTP

##  [1]  0.001318143 -0.001491822  0.005003121  0.006559683  0.063197711
##  [6] -0.154160318  0.620958327  0.628360922  0.517921517  0.598635514
## [11]  0.705512570

AIC.DTP <- 2*OPT.DTP$value + 2*length(MLE.DTP)
BIC.DTP <- 2*OPT.DTP$value + log(length(nrow(Y)))*length(MLE.DTP)

##############################################################################################
# bivariate Gaussian-copula with TPSC SAS marginals
##############################################################################################

# log likelihood function (reparameterised)
loglikTPSC = function(par){
  # Reparameterisation
  mu = par[1:2]; sigma = exp(par[3:4]); 
  epsilon = 2*expit(par[5:6])-1; delta = exp(par[7:8]);
  rho = 2*expit(par[9])-1; 
# Ingredients of the log likelihood
    norm.cop = normalCopula(rho, dim = 2)
    probs =cbind(ptp4(Y[,1],mu[1],sigma[1],epsilon[1],delta[1],param="eps",FUN=psas0),
                 ptp4(Y[,2],mu[2],sigma[2],epsilon[2],delta[2],param="eps",FUN=psas0))
    val.cop = sum(dCopula(probs, copula = norm.cop, log = TRUE))
    val.marg1 = sum( dtp4(Y[,1],mu[1],sigma[1],epsilon[1],delta[1],param="eps",FUN=dsas0,log=T) )
    val.marg2 = sum( dtp4(Y[,2],mu[2],sigma[2],epsilon[2],delta[2],param="eps",FUN=dsas0,log=T) )

    # output
    out <- -val.cop - val.marg1 - val.marg2
    
        return(out)
}


start = c(0,0,log(0.01),log(0.01),0,0,-0.5,-0.5,0) 

# Optimisation step
OPT.TPSC = optim(start,loglikTPSC,control=list(maxit=10000))
OPT.TPSC

## $par
## [1]  0.0013335816 -0.0001225573 -5.2677213114 -4.8030476208  0.1599806742
## [6]  0.0148994028 -0.4496818989 -0.4690818669  1.7432756339
## 
## $value
## [1] -10724.16
## 
## $counts
## function gradient 
##     4366       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

MLE.TPSC <- c(OPT.TPSC$par[1:2], exp(OPT.TPSC$par[3:4]), 2*expit(OPT.TPSC$par[5:6])-1, exp(OPT.TPSC$par[7:8]),  2*expit(OPT.TPSC$par[9])-1  )
MLE.TPSC

## [1]  0.0013335816 -0.0001225573  0.0051553446  0.0082047041  0.0798201678
## [6]  0.0074495636  0.6378310141  0.6255763671  0.7022053081

AIC.TPSC <- 2*OPT.TPSC$value + 2*length(MLE.TPSC)
BIC.TPSC <- 2*OPT.TPSC$value + log(length(nrow(Y)))*length(MLE.TPSC)

# Model comparison (in this case the model with TSPC marginals is favoured)

c(AIC.TPSC, AIC.DTP)

## [1] -21430.32 -21471.67

c(BIC.TPSC, BIC.DTP)

## [1] -21448.32 -21493.67

Contour plots

Let’s now use the selected model (TPSC marginals) to produce contour plots and visualise the fit of this model.

#####################################################################################################################################################
# Contour plots
#####################################################################################################################################################
# Fitted density
OPT <- OPT.TPSC

pred.den <- function(x){
  mu = OPT$par[1:2]
  sigma = exp(OPT$par[3:4])
  epsilon = 2*expit(OPT$par[5:6])-1
  delta = exp(OPT$par[7:8]);
  rho = 2*expit(OPT$par[9])-1; 
  norm.cop = normalCopula(rho, dim = 2)
  probs =cbind(ptp4(x[1],mu[1],sigma[1],epsilon[1],delta[1],param="eps",FUN=psas0),
               ptp4(x[2],mu[2],sigma[2],epsilon[2],delta[2],param="eps",FUN=psas0))
  var = dCopula(probs, copula= norm.cop, log=FALSE)*
    dtp4(x[1],mu[1],sigma[1],epsilon[1],delta[1],param="eps",FUN=dsas0)*
    dtp4(x[2],mu[2],sigma[2],epsilon[2],delta[2],param="eps",FUN=dsas0)
  return(var)
}

# Constructing the grid
x <- (-100:100)/1000
X <- grid2d(x)
z <- apply(cbind(X$x,X$y),1,pred.den)

ZD <- list(x = x, y = x, z = matrix(z, ncol = length(x)))
# Contour Density Plot:
plot(as.vector(Y[,1]),as.vector(Y[,2]),xlim=c(-0.075,0.075),ylim=c(-0.2,0.15), main="Fitted Density Contours",xlab="SMI",ylab="Swiss.Re",cex.axis=1.5,cex.lab=1.5)
par(new=TRUE)
contour(ZD,levels = c(0.75,1,2,4,8,16,32,64,128,256,512,1024,2048),xlim=c(-0.075,0.075),ylim=c(-0.2,0.15),col="red",cex.axis=1.5,cex.lab=1.5)

References

Rubio, FJ, and MFJ Steel. 2015. “Bayesian Modelling of Skewness and Kurtosis with Two-Piece Scale and Shape Distributions.” Electronic Journal of Statistics 9 (2): 1884–1912.