#######################################################################
#
# LPcopula Density R routine
#Main Ref: Mukhopadhyay, S. and Parzen, E. (2020) Nonparametric Universal Copula Modeling, 
#     Applied Stochastic Models in Business and Industry, special issue on "Data Science", 36(1), 77-94
# Last edited: June, 2020 by Deep.
#######################################################################

library(nleqslv)

LPcopula <-function(x,y,m=4,ngrid=400,select="AIC"){
	u.grid <- seq(0,1,length=ngrid)
	wd <- diff(u.grid)[1]
	TX <- eLP.poly(x,m)
	TY <- eLP.poly(y,m)
	if(select=="ALL") LP <- cov(TX,TY) else LP <- LP.smooth(cov(TX,TY),length(x),method=select)
	if(sum(LP^2)>0){
	ind.imp <- which(LP!=0,arr.ind=TRUE)
	TX.grid <- matrix(NA,ngrid,ncol(TX))
	for(k in 1:ncol(TX)){
	     A <- approxfun(as.numeric(ecdf(x)(x)),TX[,k],rule=2,method="constant",f=1,ties = "mean")
	     TX.grid[,k] <- A(u.grid)
		}
	TY.grid <- matrix(NA,ngrid,ncol(TY))
	for(k in 1:ncol(TY)){
	     A <- approxfun(as.numeric(ecdf(y)(y)),TY[,k],rule=2,method="constant",f=1,ties = "mean")
	     TY.grid[,k] <- A(u.grid)
		}
		
   #constraints
	ConFun <- function(theta){
		theta0 <- theta[1]
		Theta <- matrix(0,nrow(LP),ncol(LP))
		Theta[which(abs(LP)>0,arr.ind =TRUE)] <- theta[-1]
		P <- wd^2*exp(theta0 + TX.grid%*%Theta%*%t(TY.grid))
		con.moment <- t(TX.grid)%*%P%*%TY.grid
		mat.con <- con.moment[ind.imp]  - LP[ind.imp]
		yy <- c(mat.con)
		y.con <- rep(NA, length(yy)+1) 
		y.con[1] <- sum(P)  - 1
		y.con[2: length(y.con )] <- yy
		return(y.con)
	}
	
	#solver
	theta.l2 <- c(0,c(LP[LP!=0]))
	maxent.coef <- nleqslv(theta.l2,ConFun)$x
	R <- list()
	R$theta0 <- maxent.coef[1]
	Theta <- matrix(0,nrow(LP),ncol(LP))
	Theta[which(abs(LP)>0,arr.ind =TRUE)] <- maxent.coef[-1]
	R$Theta <-  Theta
	R$LP <-  LP
	R$x <- x
	R$y <- y
	}else{
	R <- list()
	R$theta0 <- 0
	R$Theta <-  matrix(0,ncol(TX),ncol(TY))
	R$LP <-  matrix(0,nrow(TX),ncol(TY))
	R$x <- x
	R$y <- y
	}
	#return
	return(R)
 }



Predict.LPcopula <- function(lpcop.obj, eval.points,type="MaxEnt"){
# lpcop.obj: output of the LPcopula R function()
# eval.points: pair of points where you want the density estimate value
# type: L2 or MaxEnt exponential 
	x1 <- lpcop.obj$x 
	x2 <- lpcop.obj$y
	theta0 <- lpcop.obj$theta0
	Theta <- lpcop.obj$Theta
	LP <- lpcop.obj$LP
	u1 <- eval.points[,1]
	u2 <- eval.points[,2]
    S.x1 <- eLP.poly(x1,nrow(Theta))
    S.x2 <- eLP.poly(x2,ncol(Theta))
     
    P.x1 <- matrix(NA,length(u1),ncol(S.x1))
    P.x2 <- matrix(NA,length(u2),ncol(S.x2))

	  for(k in 1:ncol(S.x1)){
	     A <- approxfun(ecdf(x1)(x1),S.x1[,k],rule=2,method="constant",f=1,ties = "mean")
	     P.x1[,k] <- A(u1)
		}
	  for(k in 1:ncol(S.x2)){
	     A <- approxfun(ecdf(x2)(x2),S.x2[,k],rule=2,method="constant",f=1,ties = "mean")
	     P.x2[,k] <- A(u2)
		}

	if(type=="L2") cop <- 1 + P.x1%*%LP%*%t(P.x2) else cop <- exp(theta0 + P.x1%*%Theta%*%t(P.x2))#--estimated copula density
	return(cop)
	}
	





#######################################################################
#
# Helping function
#
#######################################################################

eLP.poly<-function(x,m){ # x: data; m: order of LP-poly
u <- (rank(x,ties.method = c("average")) - .5)/length(x)
m <- min(length(unique(u ))-1, m )
S.mat <- as.matrix(poly(u ,df=m))
return(as.matrix(scale(S.mat)))
}


LP.smooth <- function(CR,n,method){   
  CR.s <- sort(CR^2,decreasing=TRUE,index=TRUE)$x
  aa <- rep(0,length(CR.s))
   if(method=="AIC"){ penalty <- 2}
   if(method=="BIC"){ penalty <- log(n)}
    if(method=="KIC"){ penalty <- log(n)*log(min(nrow(CR), ncol(CR)))}
   aa[1] <- CR.s[1] - penalty/n
   if(aa[1]< 0){ return(rep(0,length(CR)))  }
   if(length(aa)==1){ return(CR) }
   for(i in 2:length(CR.s)){
    aa[i] <- aa[(i-1)] + (CR.s[i] - penalty/n)
    }
  CR[CR^2<CR.s[which(aa==max(aa))]] <- 0
  return(CR)
}









#######################################################################
#
# Demo
#
#######################################################################


#### 1)  Mixed copula: Age vs IQ data: sec 6 of the paper


x <- c(rep(1,3),rep(2,3),rep(3,3),rep(4,3),rep(5,3)) #age
y <-   c(8.62, 9.94, 10.06, 9.85, 10.43, 11.31,9.98, 10.69, 11.4, 9.12, 9.89, 10.57,4.8, 9.18, 9.27) #IQ
u2 <- u1 <- seq(0.0155,1-.015,length=30)
lpcop.obj <- LPcopula(x,y,m=4)
##The LP-co
lpcop.obj$LP 

#cop-density estimate: LP-L2
Z.L2 <- Predict.LPcopula(lpcop.obj, eval.points=cbind(u1,u2),type="L2")
persp(u1,u2,Z.L2,col="rosybrown1",ticktype="detailed",theta=-55,phi=20,expand=.95,main="(b)",cex.axis=1.1,cex.main=1.4,cex.lab=1.2,xlab="u",ylab="v",zlab="")

#cop-density estimate: LP-exponential
Z.L2 <- Predict.LPcopula(lpcop.obj, eval.points=cbind(u1,u2),type="MaxEnt")
persp(u1,u2,Z.L2,col="rosybrown1",ticktype="detailed",theta=-55,phi=20,expand=.95,main="(b)",cex.axis=1.1,cex.main=1.4,cex.lab=1.2,xlab="u",ylab="v",zlab="")



#### 2)  Discrete copula: Fisher's Data 4 x 5 table: sec 3, example 3
# data preparation 
library(MASS)
data(caith)
M<- caith
rownames(M) <-seq(1:nrow(M))
colnames(M) <-seq(1:ncol(M))
library("psych") # convert table to matrix
M <- table2matrix(M) 
x<- M[,1]
y <- M[,2]
u2 <- u1 <- seq(0,1,length=35)

#analysis
u2 <- u1 <- seq(0,1,length=35)
lpcop.obj <- LPcopula(x,y,m=4,select="BIC")
##The LP-co
lpcop.obj$LP 



#cop-density estimate: LP-L2
Z.L2 <- Predict.LPcopula(lpcop.obj, eval.points=cbind(u1,u2),type="L2")
persp(u1,u2,Z.L2,col="rosybrown1",ticktype="detailed",theta=-55,phi=20,expand=.95,main="(b)",cex.axis=1.1,cex.main=1.4,cex.lab=1.2,xlab="u",ylab="v",zlab="")

#cop-density estimate: LP-exponential
Z.L2 <- Predict.LPcopula(lpcop.obj, eval.points=cbind(u1,u2),type="MaxEnt")
persp(u1,u2,Z.L2,col="rosybrown1",ticktype="detailed",theta=-55,phi=20,expand=.95,main="(b)",cex.axis=1.1,cex.main=1.4,cex.lab=1.2,xlab="u",ylab="v",zlab="")



#### 3)  Continuous copula: Insurance data: sec 3, example 1.
# data preparation
library(copula)
data(loss, package="copula")
L<-loss
head(L)
x<-L[,1]
y<-L[,2] 
#analysis
u2 <- u1 <- seq(0.05,1-.05,length=40)
lpcop.obj <- LPcopula(x,y,m=4,select="BIC")
##The LP-co
lpcop.obj$LP 


#cop-density estimate: LP-L2
Z.L2 <- Predict.LPcopula(lpcop.obj, eval.points=cbind(u1,u2),type="L2")
persp(u1,u2,Z.L2,col="rosybrown1",ticktype="detailed",theta=-55,phi=20,expand=.95,main="(b)",cex.axis=1.1,cex.main=1.4,cex.lab=1.2,xlab="u",ylab="v",zlab="")

#cop-density estimate: LP-exponential
Z.L2 <- Predict.LPcopula(lpcop.obj, eval.points=cbind(u1,u2),type="MaxEnt")
persp(u1,u2,Z.L2,col="rosybrown1",ticktype="detailed",theta=-55,phi=20,expand=.95,main="(b)",cex.axis=1.1,cex.main=1.4,cex.lab=1.2,xlab="u",ylab="v",zlab="")