https://github.com/cran/NCSCopula
Tip revision: 3909ac56a196b79173af915f5dec1582c730b991 authored by Bouchra R. Nasri on 28 November 2019, 15:50:02 UTC
version 1.0.1
version 1.0.1
Tip revision: 3909ac5
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)
}