LoglikNCSCop.R
#' @title Log-likelihood of a non-central squared copula
#'
#' @description This function computes the log-likelihood vector of a non-central squared copula
#'
#' @param alpha     unconstrained non-centrality parameters a1, a2, and unconstrained copula parameters.
#' @param U         (nx2) data matrix of pseudo-observations.
#' @param family    'Gaussian' , 't' , 'Clayton' , 'Frank' , 'Gumbel'.
#' @param p         number of different non-centrality parameters (0,1,2 default).
#'
#' @author Bouchra R. Nasri, August 14, 2019
#' @return \item{LL}{Vector of log-likelihoods}
#'
#'
#'
#' @examples alpha = c(log(0.2),log(5),log(2),log(12));
#'param = c(0.5,2.5,0.5);
#'data = SimNCSCop('Clayton', 250, param);
#'LL = LoglikNCSCop(alpha,data,'Clayton')
#'
#'
#'@export

LoglikNCSCop <- function(alpha,U,family, p=2)
{

theta = ParamCop(family,alpha[(p+1):(p+2)]);
if(p==1){
a1 = 3/(1+exp(-alpha[1]));
a2 = a1;
}   else if (p==2){
a1 = 3/(1+exp(-alpha[1]));
a2 = 3/(1+exp(-alpha[2]));
}

if(p>0){
d1 = a1^2;
d2 = a2^2;

U1 = U[,1];
U2 = U[,2];

G1 = sqrt(qchisq(U1,1,d1));  #Ga inverse
G2 = sqrt(qchisq(U2,1,d2));

ha = matrix(c((G1-a1),(-G1-a1) , (G2-a2), (-G2-a2)), ncol=4, byrow=FALSE);

htilde = pnorm(ha)+ 1E-10;
w0     = dnorm(ha);

U11 = matrix(c(htilde[,1], htilde[,3]), ncol=2,byrow=FALSE);
U12 = matrix(c(htilde[,1], htilde[,4]), ncol=2,byrow=FALSE);
U21 = matrix(c(htilde[,2], htilde[,3]), ncol=2,byrow=FALSE);
U22 = matrix(c(htilde[,2], htilde[,4]), ncol=2,byrow=FALSE);
w0011  = w0[,1] /(w0[,1]+w0[,2]);
w0012  = 1-w0011;
w0021 = w0[,3] /(w0[,3]+w0[,4]);
w0022  = 1-w0021;

w11 = w0011 * w0021;
w12 = w0011 * w0022;
w21 = w0012 * w0021;
w22 = w0012 * w0022;

switch(family,
"Gaussian" =

{
f11 = copula::dCopula(U11,copula::normalCopula(theta[1], dim = 2))+ 1E-20;
f12 = copula::dCopula(U12,copula::normalCopula(theta[1], dim = 2))+ 1E-20;
f21 = copula::dCopula(U21,copula::normalCopula(theta[1], dim = 2))+ 1E-20;
f22 = copula::dCopula(U22,copula::normalCopula(theta[1], dim = 2))+ 1E-20;

},

"t" = {
f11 = copula::dCopula(U11, copula::tCopula(theta[1], dim = 2,df = theta[2]))+ 1E-20;
f12 = copula::dCopula(U12, copula::tCopula(theta[1], dim = 2,df = theta[2]))+ 1E-20;
f21 = copula::dCopula(U21, copula::tCopula(theta[1], dim = 2,df = theta[2]))+ 1E-20;
f22 = copula::dCopula(U22, copula::tCopula(theta[1], dim = 2,df = theta[2]))+ 1E-20;
},

"Clayton" = {
f11 = copula::dCopula(U11,copula::claytonCopula(theta[1], dim = 2))+ 1E-20;
f12 = copula::dCopula(U12,copula::claytonCopula(theta[1], dim = 2))+ 1E-20;
f21 = copula::dCopula(U21,copula::claytonCopula(theta[1], dim = 2))+ 1E-20;
f22 = copula::dCopula(U22,copula::claytonCopula(theta[1], dim = 2))+ 1E-20;
},

"Frank" = {
f11 = copula::dCopula(U11,copula::frankCopula(theta[1], dim = 2))+ 1E-20;
f12 = copula::dCopula(U12,copula::frankCopula(theta[1], dim = 2))+ 1E-20;
f21 = copula::dCopula(U21,copula::frankCopula(theta[1], dim = 2))+ 1E-20;
f22 = copula::dCopula(U22,copula::frankCopula(theta[1], dim = 2))+ 1E-20;
},

"Gumbel" = {
f11 = copula::dCopula(U11,copula::gumbelCopula(theta[1], dim = 2))+ 1E-20;
f12 = copula::dCopula(U12,copula::gumbelCopula(theta[1], dim = 2))+ 1E-20;
f21 = copula::dCopula(U21,copula::gumbelCopula(theta[1], dim = 2))+ 1E-20;
f22 = copula::dCopula(U22,copula::gumbelCopula(theta[1], dim = 2))+ 1E-20;
}
)

LL = log(w11 * f11 + w12 * f12 + w21 * f21 + w22 * f22);
}
else{

switch(family,
"Gaussian" =

{
LL = log(copula::dCopula(U,copula::normalCopula(theta[1], dim = 2))+ 1E-20);

},

"t" = {
LL = log(copula::dCopula(U, copula::tCopula(theta[1], dim = 2,df = theta[2]))+ 1E-20);
},

"Clayton" = {
LL = log(copula::dCopula(U,copula::claytonCopula(theta[1], dim = 2))+ 1E-20);
},

"Frank" = {
LL = log(copula::dCopula(U,copula::frankCopula(theta[1], dim = 2))+ 1E-20);
},

"Gumbel" = {
LL = log(copula::dCopula(U,copula::gumbelCopula(theta[1], dim = 2))+ 1E-20);
}
)
}
return(LL)
}