R/hjcondvine.R

In vincenzocoia/copsupp: Vine Copulas made easy

# Assume vine array with 1:d on diagonal.
# Vine functions for mapping of (u[1],...,u[d])  to
# p[1]=u[1], p[2]=C_{2|1}(u[2]|u[1]) , ...
#         p[d]=C_{d|1:(d-1)}(u[d]}u[1:(d-1)])
# and inverse.
# Byproducts :
#  (a) conditional cdf F_{d|1:(d-1)}(x_d| x_1,...,x_{d-1}) after converting
#   x_j to u_j=F_j(x_j) for j=1,...,d
#  (b) conditional quantile F^{-1}_{d|1:(d-1)}(p|x_1,...,x_{d-1})
#  after converting x_j to u_jF_j(x_j) for j=1,...,d-1 and
#  applying F_d^{-1}(u_d) to the output u_d

# Note that if pair-copula families are permutation asymmetric,
# the code below needs to replace pcond (pcondnames)
# with pcond21 (pcond21names) and pcond12 (pcond12names),
# with small changes to the code, replacing pcond by either pcond21 or pcond12

# so below is not valid for pcondgumv pcondgumu etc where
# one of the U(0,1) variables has been reflected to ger negative depc.

#============================================================

# Versions of pcond vectors for R-vine
#   with different pair-copula family for each edge of the vine
# for each row u[1:d] of hypercube data,
# return
# u1, C_{2|1}(u2|u1), C_{3|12}(u3|u1,u2), ...  C_{d|1..d-1}(u_d|u[1:(d-1)])

#' Copula Conditional Distributions
#'
#' These functions compute the conditional distributions of the
#' last variable, given the other variables in the vine.
#'
#' parvec = vector of parameters of pair-copulas
#' udat = nxd matrix with uniform scores
#' A = dxd vine array with 1:d on diagonal
#' ntrunc = truncated level, assume >=1
#' pcondmat = matrix of names of conditional cdfs for trees 1,...,ntrunc
#'  (assuming only one needed for permutation symmetric pair-copulas)
#'   pcondmat is empty for diagonal and lower triangle,
#'    and could have ntrunc rows or be empty for rows ntrunc+1 to d-1
#' np = dxd where np[ell,j] is size for parameter th[ell,j]
#'   for pair-copula in tree ell, variables j and A[ell,j]
#'   np=0 on and below diagonal
#' Output:
#' return C_{2|1}(u2|u1), C_{3|12}(u3|u1,u2), ...  C_{d|1..d-1}(ud|u1,...u[d-1])
#' [This function is modification of rvinellkv.trunc2 in CopulaModel]
#' @author Harry Joe
#' @import CopulaModel
#' @export
#rvinellkv.trunc2=function(parvec,udat,A,ntrunc,logdcopmat,pcondmat,np)
rvinepcond=function(parvec,udat,A,ntrunc,pcondmat,np)
{ d=ncol(A)  # or ncol(udat)
  npmax=max(np); th=array(0,c(npmax,d,d)) # format of parameter for vine
  ii=0;
  for(ell in 1:ntrunc)
  { for(j in (ell+1):d)
    { ipp=1:np[ell,j]
      th[ipp,ell,j]=parvec[ii+ipp]
      ii=ii+np[ell,j]
    }
  }
  out=varray2M(A)
  M=out$mxarray
  icomp=out$icomp
  n=nrow(udat)
  #llkv=rep(0,n)
  vinepcond=matrix(0,n,d)
  vinepcond[,1]=udat[,1]
  v=matrix(0,n,d); vp=matrix(0,n,d); s=matrix(0,n,d);
  nllk=0
  # tree 1
  #for(j in 2:d)
  #{ ipp=1:np[1,j]
    #logdcop=match.fun(logdcopmat[1,j])
    #llkv=llkv+logdcop(udat[,A[1,j]],udat[,j],th[ipp,1,j])
  #}
  # tree 2
  if(ntrunc>=2)
  { for(j in 2:d)
    { ipp=1:np[1,j]
      pcond=match.fun(pcondmat[1,j])
      if(icomp[1,j]==1) vp[,j]=pcond(udat[,A[1,j]],udat[,j],th[ipp,1,j])
      v[,j]=pcond(udat[,j],udat[,A[1,j]],th[ipp,1,j])
    }
    vinepcond[,2]=v[,2]
    for(j in 3:d)
    { if(A[2,j]<M[2,j]) s[,j]=vp[,M[2,j]] else s[,j]=v[,A[2,j]] }
    #for(j in 3:d)
    #{ ipp=1:np[2,j]
      #logdcop=match.fun(logdcopmat[2,j])
      #llkv=llkv+logdcop(s[,j],v[,j],th[ipp,2,j])
    #}
    w=v; wp=vp
  }
  # remaining trees
  if(ntrunc>=3)
  { for(ell in 3:ntrunc)
    { for(j in ell:d)
      { ipp=1:np[ell-1,j]
        pcond=match.fun(pcondmat[ell-1,j])
        if(icomp[ell-1,j]==1) vp[,j]=pcond(s[,j],w[,j],th[ipp,ell-1,j])
        v[,j]=pcond(w[,j],s[,j],th[ipp,ell-1,j])
      }
      vinepcond[,ell]=v[,ell]
      for(j in (ell+1):d)
      { if(A[ell,j]<M[ell,j]) s[,j]=vp[,M[ell,j]] else s[,j]=v[,A[ell,j]] }
      #for(j in (ell+1):d)
      #{ ipp=1:np[ell,j]
        #logdcop=match.fun(logdcopmat[ell,j])
        #llkv=llkv+logdcop(s[,j],v[,j],th[ipp,ell,j])
      #}
      w=v; # wp=vp  # wp is not used
    }
  }
  if(ntrunc>1 & ntrunc<d)
  { ell=ntrunc+1
    for(j in ell:d)
    { ipp=1:np[ell-1,j]
      pcond=match.fun(pcondmat[ell-1,j])
      if(icomp[ell-1,j]==1) vp[,j]=pcond(s[,j],w[,j],th[ipp,ell-1,j])
      v[,j]=pcond(w[,j],s[,j],th[ipp,ell-1,j])
      vinepcond[,j]=v[,j]
    }
  }
  if(ntrunc==1)
  { ell=2
    for(j in ell:d)
    { ipp=1:np[ell-1,j]
      pcond=match.fun(pcondmat[ell-1,j])
      v[,j]=pcond(udat[,j],udat[,A[1,j]],th[ipp,ell-1,j])
      vinepcond[,j]=v[,j]
    }
  }
  #llkv
  vinepcond
}

#============================================================

# modification of previous R-vine simualtion code to get quantile function
# this code is not efficient but should be OK for case of
#  pair-copulas all permutation symmetric

#' n = #replications or sample size
#' pmat = nxd matrix of levels of u1, C_{2|1}(u_2|u_1),...,
#'                          C_{d|1..d-1}(u_d|u[1:(d-1)])
#' parvec = vector of parameters to be optimized in nllk
#' A = dxd vine array with 1:d on diagonal
#' ntrunc = truncated level, assume >=1
#' pcondmat = matrix of names of conditional cdfs for trees 1,...,ntrunc
#'  (assuming only one needed for permutation symmetric pair-copulas)
#' qcondmat = matrix of names of conditional quantile functions for
#'        trees 1,...,ntrunc
#'   qcondmat and pcondmat are empty for diagonal and lower triangle,
#'    and could have ntrunc rows or be empty for rows ntrunc+1 to d-1
#' np = dxd where np[ell,j] is #parameters for parameter th[ell,j]
#'   for pair-copula in tree ell, variables j and A[ell,j]
#'   np=0 on and below diagonal
#' iinv=T to check that this is inverse of rvinepcond()
#'    in this case, columns of pmat come from rvinepcond()
#' iinv=F, get quantiles C_{d|1:(d-1)}(p|u[1:(d-1)]) based on last column of
#'       pmat[,d] where u[1],..,u[d-1] have been previously converted to
#'       u[1], C_{2|1}(u[2]|u[1]), ... C_{d-1|1:(d-2)}(u[d-1]|u[1:(d-2)])
#'       via rvinepcond()
#' @return n x d matrix with values in (0,1) or
#'       quantile C_{d|1...d-1}(p| u[1:(d-1)])
#' [This function is modification of rvinesimvec2 in CopulaModel]
#' @author Harry Joe
#' @import CopulaModel
#' @export
#rvinesimvec2=function(nsim,A,ntrunc,parvec,np,qcondmat,pcondmat,iprint=F)
rvineqcond=function(pmat,A,ntrunc,parvec,np,qcondmat,pcondmat,iinv=F)
{ d=ncol(A)
  # get matrix ip1,ip2 of indices
  ii=0
  ip1=matrix(0,d,d); ip2=matrix(0,d,d)
  for(ell in 1:ntrunc)
  { for(j in (ell+1):d)
    { ip1[ell,j]=ii+1; ip2[ell,j]=ii+np[ell,j]
      ii=ii+np[ell,j]
    }
  }
  #if(iprint) { print(ip1); print(ip2) }
  out=varray2M(A)
  M=out$mxarray
  icomp=out$icomp
  #p=matrix(runif(nsim*d),nsim,d)
  p=pmat
  n=nrow(p)
  qq=array(0,c(n,d,d)); v=array(0,c(n,d,d))
  u=matrix(0,n,d)
  u[,1]=p[,1]; qq[,1,1]=p[,1]; qq[,2,2]=p[,2];
  qcond=match.fun(qcondmat[1,2])
  u[,2]=qcond(p[,2],p[,1],parvec[ip1[1,2]:ip2[1,2]])
  qq[,1,2]=u[,2]
  if(icomp[1,2]==1)
  { #pcond=match.fun(pcondnames[1])
    pcond=match.fun(pcondmat[1,2])
    v[,1,2]=pcond(u[,1],u[,2],parvec[ip1[1,2]:ip2[1,2]])
  }
  # the main loop
  for(j in 3:d)  # variable
  { tt=min(ntrunc,j-1)
    qq[,tt+1,j]=p[,j]
    if(tt>1)
    { for(ell in seq(tt,2)) # tree
      { if(A[ell,j]==M[ell,j]) { s= qq[,ell,A[ell,j]] }
        else { s=v[,ell-1,M[ell,j]] }
        #qcond=match.fun(qcondnames[ell])
        qcond=match.fun(qcondmat[ell,j])
        qq[,ell,j]=qcond(qq[,ell+1,j], s, parvec[ip1[ell,j]:ip2[ell,j]]);
      }
    }
    #qcond=match.fun(qcondnames[1])
    qcond=match.fun(qcondmat[1,j])
    qq[,1,j]=qcond(qq[,2,j],u[,A[1,j]], parvec[ip1[1,j]:ip2[1,j]])
    u[,j]=qq[,1,j]
    # set up for next iteration (not needed for last j=d)
    #pcond=match.fun(pcondnames[1])
    pcond=match.fun(pcondmat[1,j])
    v[,1,j]=pcond(u[,A[1,j]],u[,j],parvec[ip1[1,j]:ip2[1,j]])
    if(tt>1)
    { for(ell in 2:tt)
      { if(A[ell,j]==M[ell,j]) { s=qq[,ell,A[ell,j]] }
        else { s=v[,ell-1,M[ell,j]] }
        if(icomp[ell,j]==1)
        { #pcond=match.fun(pcondnames[ell])
          pcond=match.fun(pcondmat[ell,j])
          v[,ell,j]=pcond(s, qq[,ell,j], parvec[ip1[ell,j]:ip2[ell,j]]);
        }
      }
    }
  }
  if(iinv) u=u[,d]  # last column as quantiles
  u  # nxd matrix to check that this funciton is inverse of rvinepcond
}