R/dcomper_mv.R

In dsfa: Distributional Stochastic Frontier Analysis

#' Multivariate Composed-Error distribution
#'
#' Probablitiy density function, distribution, quantile function and random number generation for the multivariate composed-error distribution
#'
#' @return \code{dcomper_mv} gives the density, \code{pcomper_mv} give the distribution function, and \code{rcomper_mv} generates random numbers, with given parameters.
#' If the derivatives are calculated the output is a \code{derivs} object.
#'
#' @details A bivariate random vector \eqn{(X_1,X_2)=\boldsymbol{X}} follows a multivariate composed-errordistribution \eqn{f_{X_1,X_2}(x_1,x_2)}, which can be rewritten using Sklars' theorem via a copula
#' \deqn{f_{X_1,X_2}(y_1,y_2)=c(F_{X_1}(x_1),F_{X_2}(x_2),\delta) \cdot f_{X_1}(x_1) f_{X_2}(x_2) \qquad,}
#' where \eqn{c(\cdot)} is the density of the copula and \eqn{F_{X_m}(x_m)},\eqn{f_{X_m}(x_m)} are the marginal cdfs and pdfs respectively for \eqn{m \in \{1,2\}}. \eqn{\delta} is the copula parameter.
#'
#'
#' @param x numeric matrix of quantiles. Must have two columns.
#' @param mu numeric matrix of \eqn{\mu}. Must have two columns.
#' @param sigma_v numeric matrix of \eqn{\sigma_V}. Must be positive and have two columns.
#' @param sigma_u numeric matrix of \eqn{\sigma_U}. Must be positive and have two columns.
#' @param s integer vector of length two; each element corresponds to one marginal.
#' @param tri optional; List of objects generated by [trind_generator()].
#' @param log.p logical; if TRUE, probabilities p are given as log(p).
#' @param distr string vector of length three; the first two elements determine the distribution of the marginals. Available are:\cr
#' `normhnorm`, Normal-halfnormal distribution \cr
#' `normexp`, Normal-exponential distribution \cr
#' The last element determines the distribution of the copula:\cr
#' `independent`, Independence copula \cr
#' `normal`, Gaussian copula \cr
#' `clayton`, Clayton copula \cr
#' `gumbel`, Gumbel copula \cr
#' `frank`, Frank copula \cr
#' `joe`, Joe copula \cr
#' `amh`, Ali-Mikhail-Haq copula \cr
#' @inheritParams list2derivs
#' @inheritParams dcop
#'
#' @examples
#' pdf <- dcomper_mv(x=matrix(c(0,10),ncol=2), mu=matrix(c(1,2),ncol=2),
#'                   sigma_v=matrix(c(3,4),ncol=2), sigma_u=matrix(c(5,6),ncol=2),
#'                   delta=c(0.5), s=c(-1,-1), distr=c("normhnorm","normhnorm","normal"),
#'                   deriv=2 ,
#'                   tri=list(trind_generator(3),trind_generator(3),trind_generator(1),
#'                   trind_generator(6),trind_generator(7)), 
#'                   log.p=TRUE)
#' cdf <- pcomper_mv(q=matrix(c(0,10),ncol=2), mu=matrix(c(1,2),ncol=2),
#'                   sigma_v=matrix(c(3,4),ncol=2), sigma_u=matrix(c(5,6),ncol=2),
#'                   delta=c(0.5), s=c(-1,-1), distr=c("normhnorm","normhnorm","normal"))
#' r <- rcomper_mv(n=10, mu=matrix(c(1,2),ncol=2),
#'                 sigma_v=matrix(c(3,4),ncol=2), sigma_u=matrix(c(5,6),ncol=2),
#'                 delta=c(0.5), s=c(-1,-1), distr=c("normhnorm","normhnorm","normal"))
#'
#' @references
#' \itemize{
#' \item \insertRef{aigner1977formulation}{dsfa}
#' \item \insertRef{kumbhakar2015practitioner}{dsfa}
#' \item \insertRef{schmidt2020analytic}{dsfa}
#' \item \insertRef{gradshteyn2014table}{dsfa}
#' \item \insertRef{azzalini2013skew}{dsfa}
#' }
#' 
#' @family distribution
#' 
#' @export
dcomper_mv <- function(x, mu=matrix(c(0,0),ncol=2), sigma_v=matrix(c(1,1),ncol=2), sigma_u=matrix(c(1,1),ncol=2), delta=matrix(0,nrow=1), s=c(-1,-1), distr=c("normhnorm","normhnorm","normal"), rot=0, deriv_order=0, tri=NULL, log.p = FALSE){
  #Density function of the multivariate comper distribution
  
  X<-tryCatch(cbind(x, mu, sigma_v, sigma_u, delta), warning=function(w) {
    stop("Input vectors have incompatible lengths")})

  if(is.null(tri)){
    tri=list(trind_generator(3),trind_generator(3),trind_generator(1),trind_generator(6),trind_generator(7))
  }

  out<-dcomper_mv_cpp(x=X[,1:2, drop=F], m=X[,3:4, drop=F], v=X[,5:6, drop=F], u=X[,7:8, drop=F], delta=X[,9, drop=T], s=s, distr=distr, rot=rot, deriv_order=deriv_order, tri=tri, logp=log.p)

  #Return ouptut
  return(out)
}

#' @describeIn dcomper_mv distribution function for the multivariate composed-error distribution.
#' @param q numeric matrix of probabilities.
#' @export
pcomper_mv <- function(q, mu=matrix(c(0,0),ncol=2), sigma_v=matrix(c(1,1),ncol=2), sigma_u=matrix(c(1,1),ncol=2), delta=0, s=c(-1,-1), distr=c("normhnorm","normhnorm","normal"), rot=0, deriv_order=0, tri=NULL, log.p = FALSE){
  #Probability function of the multivariate comper distribution
  if(deriv_order>0){
    stop(paste("No derivatives implemented for the pcomper_mv", "\n", ""))
  }

  X<-tryCatch(cbind(q, mu, sigma_v, sigma_u, delta), warning=function(w) {
    stop("Input vectors have incompatible lengths")})
  
  #Margin 1
  F1<-pcomper(q=X[,1, drop=T], mu=X[,3, drop=T], sigma_v=X[,5, drop=T], sigma_u=X[,7, drop=T], s=s[1], distr=distr[1])

  #Margin 2
  F2<-pcomper(q=X[,2, drop=T], mu=X[,4, drop=T], sigma_v=X[,6, drop=T], sigma_u=X[,8, drop=T], s=s[2], distr=distr[2])

  #Evaluate cdf of copula at probability integral transformed observations
  out<-pcop(W=cbind(F1,F2), delta=X[,9, drop=T], distr=distr[3], rot=rot, log.p = log.p)
  names(out)<-NULL

  #Return ouptut
  return(out)
}


#' @describeIn dcomper_mv random number generation for the multivariate composed-error distribution.
#' @param n positive integer; number of observations.
#' @export
rcomper_mv <- function(n, mu=matrix(c(0,0),ncol=2), sigma_v=matrix(c(1,1),ncol=2), sigma_u=matrix(c(1,1),ncol=2), delta=matrix(0,nrow=1), s=c(-1,-1), distr=c("normhnorm","normhnorm","normal"), rot=0){
  #Function to generate n random numbers from the comper distribution
  X<-tryCatch(cbind(0,0, mu, sigma_v, sigma_u, delta), warning=function(w) {
    stop("Input vectors have incompatible lengths")})
  
  #Generate pseudo observations
  W<-rcop(n=n, delta=X[,9, drop=T], distr=distr[3], rot=rot)

  #Margin 1
  y1<-qcomper(p=W[,1], mu=X[,3, drop=T], sigma_v=X[,5, drop=T], sigma_u=X[,7, drop=T], s=s[1], distr=distr[1])

  #Margin 2
  y2<-qcomper(p=W[,2], mu=X[,4, drop=T], sigma_v=X[,6, drop=T], sigma_u=X[,8, drop=T], s=s[2], distr=distr[2])

  #Combine y1 and y2
  out<-cbind(y1,y2)

  #Return output
  return(out)
}