library(emdbook)

############################# g-h transform#######
gh_trans <- function(u,g,h){
if(g!=0) (exp(g*u)-1)*exp(h/2*u^2)/g else u*exp(h/2*u^2)
}

fstar <- function(zp,g,h){
u <- zp
if(g!=0)
tmp <- -u^2*(1+h)/2-log(exp(g*u)+(exp(g*u)-1)*h*u/g) else tmp <- -u^2*(1+h)/2-log(1+h*u^2)
tmp
}


fstar_hessian <- function(zp,xi,omega,g,h){

obj <- zp_hessian(zp,xi,omega,g,h)

d1Z <- obj$gradient
d2Z <- obj$hessian
if(abs(g)< 10^(-4)) g <- 10^(-4)

t <- exp(g*zp)

tmp1 <- t + g^(-1)*(t-1)*h*zp
tmp2 <- g*t+t*h*zp+g^(-1)*(t-1)*h
tmp3 <- zp*t-g^(-2)*(t-1)*h*zp + g^(-1)*t*h*zp^2
tmp4 <- g^(-1)*(t-1)*zp

tmp1.xi <- tmp2*d1Z[1,]
tmp1.omega <- tmp2*d1Z[2,]
tmp1.g <- tmp3+tmp2*d1Z[3,]
tmp1.h <- tmp4+tmp2*d1Z[4,]


d1f <- d1temp <- rbind(tmp1.xi,tmp1.omega,tmp1.g,tmp1.h)

d2temp <- d2Z

tmp21 <- g^2*t + 2*t*h + g*t*h*zp
tmp22 <- g*zp*t+t+t*h*zp^2-g^(-2)*(t-1)*h+g^(-1)*t*h*zp
tmp23 <- g^(-1)*(t-1)+t*zp
tmp31 <- t+g*zp*t-g^(-1)*t*h*zp-g^(-2)*(t-1)*h+t*h*zp^2+ g^(-1)*t*2*h*zp
tmp2.g <- tmp22+tmp21*d1Z[3,]
tmp2.h <- tmp23+tmp21*d1Z[4,]
tmp3.g <- zp^2*t+2*g^(-3)*(t-1)*h*zp-2*g^(-2)*t*h*zp^2+g^(-1)*t*h*zp^3+tmp31*d1Z[3,]
tmp3.h <- -g^(-2)*(t-1)*zp+g^(-1)*t*zp^2+tmp31*d1Z[4,]
tmp4.h <- (t*zp+g^(-1)*(t-1))*d1Z[4,]


d2temp[1,1,] <- tmp1.xi.xi <- tmp21*d1Z[1,]^2+tmp2*d2Z[1,1,]
d2temp[1,2,] <- tmp1.xi.omega <- tmp21*d1Z[1,]*d1Z[2,]+tmp2*d2Z[1,2,]
d2temp[1,3,] <- tmp1.xi.g <- tmp2.g*d1Z[1,]+tmp2*d2Z[1,3,]
d2temp[1,4,] <- tmp1.xi.h <- tmp2.h*d1Z[1,]+tmp2*d2Z[1,4,]
d2temp[2,2,] <- tmp1.omega.omega <- tmp21*d1Z[2,]^2+tmp2*d2Z[2,2,]
d2temp[2,3,] <- tmp1.omega.g <- tmp2.g*d1Z[2,]+tmp2*d2Z[2,3,]
d2temp[2,4,] <- tmp1.omega.h <- tmp2.h*d1Z[2,]+tmp2*d2Z[2,4,]
d2temp[3,3,] <- tmp1.g.g <- tmp3.g + tmp2.g*d1Z[3,] + tmp2*d2Z[3,3,]
d2temp[3,4,] <- tmp1.g.h <- tmp3.h + tmp2.h*d1Z[3,] + tmp2*d2Z[3,4,]
d2temp[4,4,] <- tmp1.h.h <- tmp4.h + tmp2.h*d1Z[4,] + tmp2*d2Z[4,4,]

d1f[1,] <- -(1+h)*d1Z[1,]*zp-1/tmp1*d1temp[1,]
d1f[2,] <- -(1+h)*d1Z[2,]*zp-1/tmp1*d1temp[2,]-1/omega
d1f[3,] <- -(1+h)*d1Z[3,]*zp-1/tmp1*d1temp[3,]
d1f[4,] <- -(1+h)*d1Z[4,]*zp-1/tmp1*d1temp[4,]-1/2*zp^2


d2f <- d2temp

d2f[1,1,] <- -(1+h)*(d1Z[1,]^2+zp*d2Z[1,1,])+d1temp[1,]^2/tmp1^2-d2temp[1,1,]/tmp1
d2f[2,1,] <- d2f[1,2,] <- -(1+h)*(d1Z[1,]*d1Z[2,]+zp*d2Z[1,2,])+d1temp[1,]*d1temp[2,]/tmp1^2-d2temp[1,2,]/tmp1
d2f[3,1,] <- d2f[1,3,] <- -(1+h)*(d1Z[1,]*d1Z[3,]+zp*d2Z[1,3,])+d1temp[1,]*d1temp[3,]/tmp1^2-d2temp[1,3,]/tmp1
d2f[4,1,] <- d2f[1,4,] <- -zp*d1Z[1,] -(1+h)*(d1Z[1,]*d1Z[4,]+zp*d2Z[1,4,])+d1temp[1,]*d1temp[4,]/tmp1^2-d2temp[1,4,]/tmp1
d2f[2,2,] <- -(1+h)*(d1Z[2,]^2+zp*d2Z[2,2,])+d1temp[2,]^2/tmp1^2-d2temp[2,2,]/tmp1+1/omega^2
d2f[3,2,] <- d2f[2,3,] <- -(1+h)*(d1Z[2,]*d1Z[3,]+zp*d2Z[2,3,])+d1temp[2,]*d1temp[3,]/tmp1^2-d2temp[2,3,]/tmp1
d2f[4,2,] <- d2f[2,4,] <- -zp*d1Z[2,] -(1+h)*(d1Z[2,]*d1Z[4,]+zp*d2Z[2,4,])+d1temp[2,]*d1temp[4,]/tmp1^2-d2temp[2,4,]/tmp1
d2f[3,3,] <- -(1+h)*(d1Z[3,]^2+zp*d2Z[3,3,])+d1temp[3,]^2/tmp1^2-d2temp[3,3,]/tmp1
d2f[4,3,] <- d2f[3,4,] <- -zp*d1Z[3,] -(1+h)*(d1Z[3,]*d1Z[4,]+zp*d2Z[3,4,])+d1temp[3,]*d1temp[4,]/tmp1^2-d2temp[3,4,]/tmp1
d2f[4,4,] <- -2*zp*d1Z[4,] -(1+h)*(d1Z[4,]*d1Z[4,]+zp*d2Z[4,4,])+d1temp[4,]*d1temp[4,]/tmp1^2-d2temp[4,4,]/tmp1


gradient <- c(-rowSums(d1f,na.rm=TRUE))
hessian <- as.matrix(-apply(d2f,c(1,2),sum,na.rm=TRUE))
w1 <- c(1,omega,1,1)
w2 <- diag(gradient*c(0,omega,0,0))
gradient <- gradient*w1
hessian <- hessian*outer(w1,w1)+w2
list(gradient=gradient, hessian=hessian)
}

fstar_hessian_gnull <- function(zp,xi,omega,h){

obj <- zp_hessian.gnull(zp,xi,omega,h)

d1Z <- obj$gradient
d2Z <- obj$hessian


tmp1 <- 1 + h*zp^2
tmp2 <- 2*h*zp
tmp4 <- zp^2

tmp1.xi <- tmp2*d1Z[1,]
tmp1.omega <- tmp2*d1Z[2,]
tmp1.h <- tmp4+tmp2*d1Z[3,]


d1f <- d1temp <- rbind(tmp1.xi,tmp1.omega,tmp1.h)

d2temp <- d2Z

tmp21 <- 2*h 
tmp23 <- 2*zp
tmp2.h <- tmp23+tmp21*d1Z[3,]
tmp4.h <- 2*zp*d1Z[3,]


d2temp[1,1,] <- tmp1.xi.xi <- tmp21*d1Z[1,]^2+tmp2*d2Z[1,1,]
d2temp[1,2,] <- tmp1.xi.omega <- tmp21*d1Z[1,]*d1Z[2,]+tmp2*d2Z[1,2,]
d2temp[1,3,] <- tmp1.xi.h <- tmp2.h*d1Z[1,]+tmp2*d2Z[1,3,]
d2temp[2,2,] <- tmp1.omega.omega <- tmp21*d1Z[2,]^2+tmp2*d2Z[2,2,]
d2temp[2,3,] <- tmp1.omega.h <- tmp2.h*d1Z[2,]+tmp2*d2Z[2,3,]
d2temp[3,3,] <- tmp1.h.h <- tmp4.h + tmp2.h*d1Z[3,] + tmp2*d2Z[3,3,]

d1f[1,] <- -(1+h)*d1Z[1,]*zp-1/tmp1*d1temp[1,]
d1f[2,] <- -(1+h)*d1Z[2,]*zp-1/tmp1*d1temp[2,]-1/omega
d1f[3,] <- -(1+h)*d1Z[3,]*zp-1/tmp1*d1temp[3,]-1/2*zp^2


d2f <- d2temp

d2f[1,1,] <- -(1+h)*(d1Z[1,]^2+zp*d2Z[1,1,])+d1temp[1,]^2/tmp1^2-d2temp[1,1,]/tmp1
d2f[2,1,] <- d2f[1,2,] <- -(1+h)*(d1Z[1,]*d1Z[2,]+zp*d2Z[1,2,])+d1temp[1,]*d1temp[2,]/tmp1^2-d2temp[1,2,]/tmp1
d2f[3,1,] <- d2f[1,3,] <- -zp*d1Z[1,] -(1+h)*(d1Z[1,]*d1Z[3,]+zp*d2Z[1,3,])+d1temp[1,]*d1temp[3,]/tmp1^2-d2temp[1,3,]/tmp1
d2f[2,2,] <- -(1+h)*(d1Z[2,]^2+zp*d2Z[2,2,])+d1temp[2,]^2/tmp1^2-d2temp[2,2,]/tmp1+1/omega^2
d2f[3,2,] <- d2f[2,3,] <- -zp*d1Z[2,] -(1+h)*(d1Z[2,]*d1Z[3,]+zp*d2Z[2,3,])+d1temp[2,]*d1temp[3,]/tmp1^2-d2temp[2,3,]/tmp1
d2f[3,3,] <- -2*zp*d1Z[3,] -(1+h)*(d1Z[3,]*d1Z[3,]+zp*d2Z[3,3,])+d1temp[3,]*d1temp[3,]/tmp1^2-d2temp[3,3,]/tmp1


gradient <- c(-rowSums(d1f,na.rm=TRUE))
hessian <- as.matrix(-apply(d2f,c(1,2),sum,na.rm=TRUE))
w1 <- c(1,omega,1)
w2 <- diag(gradient*c(0,omega,0))
gradient <- gradient*w1
hessian <- hessian*outer(w1,w1)+w2
list(gradient=gradient, hessian=hessian)
}

fstar_partial <- function(zp,xi,omega,g,h){
t <- exp(g*zp)
temp <- t + (t-1)*h*zp/g
zp_p <- zp_partial(zp,xi,omega,g,h)
fs_d1 <- fstar_d1(zp,g,h)
partial_xi <- zp_p[1,]*fs_d1
partial_omega <- zp_p[2,]*fs_d1-1/omega
partial_g <- zp_p[3,]*fs_d1-(t*zp-(t-1)*h*zp/g^2+t*h*zp^2/g)/temp
partial_h <- zp_p[4,]*fs_d1-((t-1)*zp/g)/temp -1/2*zp^2
rbind(partial_xi,partial_omega*omega,partial_g,partial_h*h*(1-h))
}


zp_hessian <- function(zp,xi,omega,g,h){
n <- length(zp)
if(abs(g)< 10^(-4)) g <- 10^(-4)
t <- exp(g*zp)
temp <- t + (t-1)*h*zp/g
d1T.z <- exp(h*zp^2/2)*temp*omega
d1T.xi <- 1
d1T.omega <- (t-1)*exp(h*zp^2/2)/g
d1T.g <- -omega*g^(-2)*(t-1)*exp(h*zp^2/2)+omega*t*zp*exp(h*zp^2/2)/g
d1T.h <- omega*(t-1)*exp(h*zp^2/2)*zp^2/2/g

d2T.z.z <- omega*exp(h*zp^2/2)*(temp*h*zp+g*t+t*h*zp+(t-1)*h/g)
d2T.z.xi <- 0
d2T.z.omega <- exp(h*zp^2/2)*temp
d2T.z.g <- omega*exp(h*zp^2/2)*(zp*t-g^(-2)*(t-1)*h*zp+g^(-1)*t*h*zp^2)
d2T.z.h <- omega*exp(h*zp^2/2)*(temp*zp^2/2+g^(-1)*(t-1)*zp)
d2T.xi.xi <- 0
d2T.xi.omega <- 0
d2T.xi.g <- 0
d2T.xi.h <- 0
d2T.omega.omega <- 0
d2T.omega.g <- -g^(-2)*(t-1)*exp(h*zp^2/2)+g^(-1)*t*zp*exp(h*zp^2/2)
d2T.omega.h <- g^(-1)*(t-1)*exp(h*zp^2/2)*zp^2/2
d2T.g.g <- omega*exp(h*zp^2/2)*(2*g^(-3)*(t-1)-2*g^(-2)*t*zp+g^(-1)*t*zp^2)
d2T.g.h <- omega*exp(h*zp^2/2)*(-g^(-2)*(t-1)+g^(-1)*t*zp)*zp^2/2
d2T.h.h <- g^(-1)*(t-1)*omega*exp(h*zp^2/2)*(zp^2/2)^2

d1Z.xi <- -d1T.xi/d1T.z
d1Z.omega <- -d1T.omega/d1T.z
d1Z.g <- -d1T.g/d1T.z
d1Z.h <- -d1T.h/d1T.z

d2Z <- array(NA,dim=c(4,4,n))
d2Z[1,1,] <- d2Z.xi.xi <- (d1T.z)^(-2)*((d2T.z.xi+d2T.z.z*d1Z.xi)*d1T.xi-d1T.z*(d2T.xi.xi+d2T.z.xi*d1Z.xi))
d2Z[2,1,] <- d2Z[1,2,] <- d2Z.xi.omega <- (d1T.z)^(-2)*((d2T.z.omega+d2T.z.z*d1Z.omega)*d1T.xi-d1T.z*(d2T.xi.omega+d2T.z.xi*d1Z.omega))
d2Z[3,1,] <- d2Z[1,3,] <- d2Z.xi.g <- (d1T.z)^(-2)*((d2T.z.g+d2T.z.z*d1Z.g)*d1T.xi-d1T.z*(d2T.xi.g+d2T.z.xi*d1Z.g))
d2Z[4,1,] <- d2Z[1,4,] <- d2Z.xi.h <- (d1T.z)^(-2)*((d2T.z.h+d2T.z.z*d1Z.h)*d1T.xi-d1T.z*(d2T.xi.h+d2T.z.xi*d1Z.h))
d2Z[2,2,] <- d2Z.omega.omega<- (d1T.z)^(-2)*((d2T.z.omega+d2T.z.z*d1Z.omega)*d1T.omega-d1T.z*(d2T.omega.omega+d2T.z.omega*d1Z.omega))
d2Z[3,2,] <- d2Z[2,3,] <- d2Z.omega.g <- (d1T.z)^(-2)*((d2T.z.g+d2T.z.z*d1Z.g)*d1T.omega-d1T.z*(d2T.omega.g+d2T.z.omega*d1Z.g))
d2Z[4,2,] <- d2Z[2,4,] <- d2Z.omega.h <- (d1T.z)^(-2)*((d2T.z.h+d2T.z.z*d1Z.h)*d1T.omega-d1T.z*(d2T.omega.h+d2T.z.omega*d1Z.h))
d2Z[3,3,] <- d2Z.g.g<- (d1T.z)^(-2)*((d2T.z.g+d2T.z.z*d1Z.g)*d1T.g-d1T.z*(d2T.g.g+d2T.z.g*d1Z.g))
d2Z[4,3,] <- d2Z[3,4,] <- d2Z.g.h<- (d1T.z)^(-2)*((d2T.z.h+d2T.z.z*d1Z.h)*d1T.g-d1T.z*(d2T.g.h+d2T.z.g*d1Z.h))
d2Z[4,4,] <- d2Z.h.h <- (d1T.z)^(-2)*((d2T.z.h+d2T.z.z*d1Z.h)*d1T.h-d1T.z*(d2T.h.h+d2T.z.h*d1Z.h))

d1Z <- rbind(d1Z.xi,d1Z.omega,d1Z.g, d1Z.h)

list(gradient=d1Z, hessian=d2Z)
}


zp_hessian.gnull <- function(zp,xi,omega,h){
n <- length(zp)
t <- 1
temp <- 1+h*zp^2
d1T.z <- exp(h*zp^2/2)*temp*omega
d1T.xi <- 1
d1T.omega <- zp*exp(h*zp^2/2)
d1T.h <- omega*exp(h*zp^2/2)*zp^3/2

d2T.z.z <- omega*exp(h*zp^2/2)*(temp*h*zp+2*h*zp)
d2T.z.xi <- 0
d2T.z.omega <- exp(h*zp^2/2)*temp
d2T.z.h <- omega*exp(h*zp^2/2)*(temp*zp^2/2+zp^2)
d2T.xi.xi <- 0
d2T.xi.omega <- 0
d2T.xi.h <- 0
d2T.omega.omega <- 0
d2T.omega.h <- exp(h*zp^2/2)*zp^3/2
d2T.h.h <- zp*omega*exp(h*zp^2/2)*(zp^2/2)^2

d1Z.xi <- -d1T.xi/d1T.z
d1Z.omega <- -d1T.omega/d1T.z
d1Z.h <- -d1T.h/d1T.z

d2Z <- array(NA,dim=c(3,3,n))
d2Z[1,1,] <- d2Z.xi.xi <- (d1T.z)^(-2)*((d2T.z.xi+d2T.z.z*d1Z.xi)*d1T.xi-d1T.z*(d2T.xi.xi+d2T.z.xi*d1Z.xi))
d2Z[2,1,] <- d2Z[1,2,] <- d2Z.xi.omega <- (d1T.z)^(-2)*((d2T.z.omega+d2T.z.z*d1Z.omega)*d1T.xi-d1T.z*(d2T.xi.omega+d2T.z.xi*d1Z.omega))
d2Z[3,1,] <- d2Z[1,3,] <- d2Z.xi.h <- (d1T.z)^(-2)*((d2T.z.h+d2T.z.z*d1Z.h)*d1T.xi-d1T.z*(d2T.xi.h+d2T.z.xi*d1Z.h))
d2Z[2,2,] <- d2Z.omega.omega<- (d1T.z)^(-2)*((d2T.z.omega+d2T.z.z*d1Z.omega)*d1T.omega-d1T.z*(d2T.omega.omega+d2T.z.omega*d1Z.omega))
d2Z[3,2,] <- d2Z[2,3,] <- d2Z.omega.h <- (d1T.z)^(-2)*((d2T.z.h+d2T.z.z*d1Z.h)*d1T.omega-d1T.z*(d2T.omega.h+d2T.z.omega*d1Z.h))
d2Z[3,3,] <- d2Z.h.h <- (d1T.z)^(-2)*((d2T.z.h+d2T.z.z*d1Z.h)*d1T.h-d1T.z*(d2T.h.h+d2T.z.h*d1Z.h))

d1Z <- rbind(d1Z.xi,d1Z.omega, d1Z.h)

list(gradient=d1Z, hessian=d2Z)
}

###Evaluating the g-h likelihood function

gh_likefun <- function(y,par,nbreaks=1000){
n <- length(y)
cutoff <- 1/n^2
ymin <- min(y)
ymax <- max(y)

xi <- par[1]
omega <- par[2]
g <- par[3]
h <- par[4]
zr <- numeric(2)

zr[1] <- -10
zr[2] <- 10

tmp <- gh_trans(c(zr[1],zr[2]),g,h)*omega+xi
if(tmp[1]<ymin&&tmp[2]>ymax){
u <- sort((y-xi)/omega)
zp <- seq(zr[1],zr[2],l=100)
ghp <- gh_trans(zp,g,h)

zr1 <- max(zp[(u[1]>ghp)])
zr2 <- min(zp[max(u)<ghp])
zp <- seq(zr1,zr2,l=nbreaks)
#zp <- seq(zr[1],zr[2],l=nbreaks)
zp <- sort(unique(zp))

ghp <- gh_trans(zp,g,h)
dzp <- diff(zp,lag=1)
dfghp <- diff(ghp,lag=1)
slop <- dfghp/dzp
loc <- .bincode(u,breaks=ghp)
uL <- ghp[loc]
zpL <- zp[loc]
zhat <-zpL+(u-uL)/dfghp[loc]*dzp[loc]

ghdq.hat <- fstar(zhat,g,h)
likelihood <- -sum(ghdq.hat)
value <- likelihood+n*log(omega)
}else{
value <- 10^9
}
return(value)
}


#######Check whether the initial values satisfy the condition
Check_initials <- function(y,par0,cutoff){
bnd <- 10
zcutL <- max(bnd,abs(qnorm(cutoff)))
tmp <- gh_trans(c(-zcutL,zcutL),par0[3],par0[4])*par0[2]+par0[1]
ymin <- min(y)
ymax <- max(y)
check <- FALSE
if(tmp[1]<ymin&&tmp[2]>ymax) check=TRUE
check
}

###Get initial values###
Initials <- function(y){
n<-length(y)
cutoff <- 1/n^2
p <- seq(0.05,0.95,l=10)
par0 <- LVII(y,p)
if(par0[4]<0) par0[4]<-0.1
check <- Check_initials(y,par0,cutoff)
while(!check){
par0[4] <- par0[4]+runif(1,0,0.1)
par0[2] <- par0[2]*1.5
check <- Check_initials(y,par0,cutoff)
}
par0
}

###Get initial values when h=0###
Initials_hnull <- function(y){
n<-length(y)
cutoff <- 1/n^2
p <- seq(0.05,0.95,l=10)
par0 <- LVII(y,p)
par0[4] <- 0.1
check <- Check_initials(y,par0,cutoff)
while(!check){
par0[2] <- par0[2]*c(1.5)
check <- Check_initials(y,par0,cutoff)
}
par0
}

###Get initial values when g=0###
Initials_gnull <- function(y){
n<-length(y)
cutoff <- 1/n^2
p <- seq(0.05,0.95,l=10)
par0 <- LVII(y,p)
par0[3] <- 0
if(par0[4]<0) par0[4]<-0.1
check <- Check_initials(y,par0,cutoff)
while(!check){
par0[2] <- par0[2]*c(1.5)
check <- Check_initials(y,par0,cutoff)
}
par0
}

##############LV method for estimating g-h Version II########
LVII <- function(y,p){
K <- length(p)
xi.hat <- median(y)
y.p <- quantile(y,p)
y.p1 <- quantile(y,1-p)
z.p <- qnorm(p)
g.p <- -1/z.p*log((y.p1-xi.hat)/(xi.hat-y.p))
ind <- numeric(K)
g.hat <- median(g.p)
tmp <- log(g.hat*(y.p-y.p1)/(exp(g.hat*z.p)-exp(-g.hat*z.p)))
tmp.z <- z.p^2/2
obj <- lm(tmp~tmp.z)
omega.hat <- exp(coef(obj)[1])
h.hat <- c(coef(obj)[2])

fn <- function(par1){
h <- par1[2]
sum((tmp-par1[1]-h*tmp.z)^2)
}

par0 <- c(log(omega.hat),log(abs(h.hat)))
par1.hat <- optim(par0,fn)$par
omega.hat <- exp(par1.hat[1])
h.hat <- par1.hat[2]
par <- c(xi.hat,omega.hat,g.hat,h.hat)
names(par)=c("xi","omega","g","h")
par
}


##minimizing the approximated likelihood function for the case h=0
gh_likest_linear_hnull <- function(y,par0=NULL,nbreaks=1000){

n <- length(y)
ymin <- min(y)
ymax <- max(y)

#initial values for b_n
zr <- numeric(2)
zr[1] <- -10
zr[2] <- 10

if(is.null(par0)){
 par0 <- Initials_hnull(y)
 par0 <- par0[-4]
}

##Approximated likelihood function
fn <- function(par){
xi <- par[1]
omega <- exp(par[2])
g <- par[3]
h <- 0
tmp <- gh_trans(c(zr[1],zr[2]),g,h)*omega+xi
if(tmp[1]<ymin&&tmp[2]>ymax){
u <- sort((y-xi)/omega)
zp <- seq(zr[1],zr[2],l=100)
ghp <- gh_trans(zp,g,h)
####refine the value for b_n##
zr1 <- max(zp[(u[1]>ghp)])
zr2 <- min(zp[max(u)<ghp])
zp <- seq(zr1,zr2,l=nbreaks)
zp <- sort(unique(zp))
ghp <- gh_trans(zp,g,h)
dzp <- diff(zp,lag=1)
dfghp <- diff(ghp,lag=1)
slop <- dfghp/dzp
loc <- .bincode(u,breaks=ghp)
uL <- ghp[loc]
zpL <- zp[loc]
zhat <-zpL+(u-uL)/dfghp[loc]*dzp[loc]
ghdq.hat <- fstar(zhat,g,h)
likelihood <- -sum(ghdq.hat)
value <- likelihood+n*log(omega)
}else{
value <- 10^9
}
return(value)
}

##the gradient of the approximated likelihood function
grn <- function(par){
xi <- par[1]
omega <- exp(par[2])
g <- par[3]
h <- 0
#initial values for b_n
zr <- numeric(2)
zr[1] <- -10
zr[2] <- 10

tmp <- gh_trans(c(zr[1],zr[2]),g,h)*omega+xi
if(tmp[1]<ymin&&tmp[2]>ymax){
u <- sort((y-xi)/omega)
zp <- seq(zr[1],zr[2],l=100)
ghp <- gh_trans(zp,g,h)

####refine the value for b_n##
zr1 <- max(zp[(u[1]>ghp)])
zr2 <- min(zp[max(u)<ghp])

zp <- seq(zr1,zr2,l=nbreaks)
zp <- sort(unique(zp))

ghp <- gh_trans(zp,g,h)
dzp <- diff(zp,lag=1)
dfghp <- diff(ghp,lag=1)
slop <- dfghp/dzp
loc <- .bincode(u,breaks=ghp)
uL <- ghp[loc]
zpL <- zp[loc]
zhat <-zpL+(u-uL)/dfghp[loc]*dzp[loc]

ghdq.hat <- fstar(zhat,g,h)
likelihood <- -sum(ghdq.hat)
value <- likelihood+n*log(omega)
obj2 <- fstar_hessian(zhat,xi,omega,g,h)
partial.hat <- obj2$gradient[-4]
}else{
value <- 10^9
partial.hat <-  rnorm(3)*100
}
return(partial.hat)
}

par0[2] <- log(par0[2])
obj1 <-optim(par0,fn,gr=grn,method="L-BFGS-B",lower = c(-Inf,-10,-Inf),upper=c(Inf,10,Inf),hessian=TRUE) 
par<- obj1$par
gradient <- grn(par)
par[2] <- exp(par[2])
value <- gh_likefun(y,c(par,0),nbreaks=nbreaks)
if(value==10^9) converge <- FALSE else converge <- TRUE
list(par=par,likelihood=value,converge=obj1$converge,gradient=gradient,hessian=obj1$hessian,message=obj1$message)
}

##minimizing the approximated likelihood function for the case g=0
gh_likest_linear_gnull <- function(y,par0=NULL,nbreaks=1000){
n <- length(y)
cutoff <- 1/n^2
alpha <- 0.1
nbs = floor(nbreaks/3)
p <- seq(alpha,1-alpha,l=nbs)
bnd <- 10
zcutL <- max(bnd,abs(qnorm(cutoff)))
zcutU <- abs(qnorm(alpha))
zp <- seq(-zcutL,zcutL,l=nbreaks)
zp <- sort(unique(zp))
ymin <- min(y)
ymax <- max(y)

if(is.null(par0)){
 par0 <- Initials_gnull(y)
 par0 <- par0[-3]
}

fn <- function(par){
xi <- par[1]
omega <- exp(par[2])
g <- 0
h <- par[3]
zr <- numeric(2)

zr[1] <- -10
zr[2] <- 10
tmp <- gh_trans(c(zr[1],zr[2]),g,h)*omega+xi
if(tmp[1]<ymin&&tmp[2]>ymax){
u <- sort((y-xi)/omega)
zp <- seq(zr[1],zr[2],l=100)
ghp <- gh_trans(zp,g,h)

zr1 <- max(zp[(u[1]>ghp)])
zr2 <- min(zp[max(u)<ghp])
zp <- seq(zr1,zr2,l=nbreaks)

ghp <- gh_trans(zp,g,h)
dzp <- diff(zp,lag=1)
dfghp <- diff(ghp,lag=1)
slop <- dfghp/dzp
loc <- .bincode(u,breaks=ghp)
uL <- ghp[loc]
zpL <- zp[loc]
zhat <-zpL+(u-uL)/dfghp[loc]*dzp[loc]
ghdq.hat <- fstar(zhat,g,h)
likelihood <- -sum(ghdq.hat)
value <- likelihood+n*log(omega)
}else{
value <- 10^9
}
return(value)
}

grn <- function(par){
xi <- par[1]
omega <- exp(par[2])
g <- 0
h <- par[3]
zr <- numeric(2)

zr[1] <- -10
zr[2] <- 10

tmp <- gh_trans(c(zr[1],zr[2]),g,h)*omega+xi
if(tmp[1]<ymin&&tmp[2]>ymax){
u <- sort((y-xi)/omega)

zp <- seq(zr[1],zr[2],l=100)
ghp <- gh_trans(zp,g,h)

zr1 <- max(zp[(u[1]>ghp)])
zr2 <- min(zp[max(u)<ghp])
zp <- seq(zr1,zr2,l=nbreaks)

ghp <- gh_trans(zp,g,h)
dzp <- diff(zp,lag=1)
dfghp <- diff(ghp,lag=1)
slop <- dfghp/dzp
loc <- .bincode(u,breaks=ghp)
uL <- ghp[loc]
zpL <- zp[loc]
zhat <-zpL+(u-uL)/dfghp[loc]*dzp[loc]

ghdq.hat <- fstar(zhat,g,h)
likelihood <- -sum(ghdq.hat)
value <- likelihood+n*log(omega)
obj2 <- fstar_hessian_gnull(zhat,xi,omega,h)
partial.hat <- obj2$gradient
}else{
value <- 10^9
partial.hat <-  rnorm(3)*1000
}
return(partial.hat)
}
names(par0) <- c("xi","omega","h")
par0[2] <- log(par0[2])

obj1 <-optim(par0,fn,gr=grn,method="L-BFGS-B",lower = c(-Inf,-10,0),upper=c(Inf,10,Inf),hessian=TRUE) 
par<- obj1$par
gradient <- grn(par)
par[2] <- exp(par[2])
value <- gh_likefun(y,c(par[1:2],0,par[3]),nbreaks=nbreaks)
if(value==10^9) converge <- FALSE else converge <- TRUE
list(par=par,likelihood=value,converge=obj1$converge,gradient=gradient,hessian=obj1$hessian,message=obj1$message)
}

##minimizing the approximated likelihood function for a general case
gh_likest_linear<- function(y,par0=NULL,nbreaks=1000){
n <- length(y)
ymin <- min(y)
ymax <- max(y)

zr <- numeric(2)
zr[1] <- -10
zr[2] <- 10

if(is.null(par0)) par0 <- Initials(y)


fn <- function(par){
xi <- par[1]
omega <- exp(par[2])
g <- par[3]
h <- max(par[4],0)

tmp <- gh_trans(c(zr[1],zr[2]),g,h)*omega+xi
if(tmp[1]<ymin&&tmp[2]>ymax){
u <- sort((y-xi)/omega)
zp <- seq(zr[1],zr[2],l=100)
ghp <- gh_trans(zp,g,h)

zr1 <- max(zp[(u[1]>ghp)])
zr2 <- min(zp[max(u)<ghp])
zp <- seq(zr1,zr2,l=nbreaks)
zp <- sort(unique(zp))

ghp <- gh_trans(zp,g,h)
dzp <- diff(zp,lag=1)
dfghp <- diff(ghp,lag=1)
slop <- dfghp/dzp
loc <- .bincode(u,breaks=ghp)
uL <- ghp[loc]
zpL <- zp[loc]
zhat <-zpL+(u-uL)/dfghp[loc]*dzp[loc]
ghdq.hat <- fstar(zhat,g,h)
likelihood <- -sum(ghdq.hat)
value <- likelihood+n*log(omega)
}else{

value <- 10^9
}
return(value)
}

grn <- function(par){
xi <- par[1]
omega <- exp(par[2])
g <- par[3]
h <- max(par[4],0)

tmp <- gh_trans(c(zr[1],zr[2]),g,h)*omega+xi
if(tmp[1]<ymin&&tmp[2]>ymax){
u <- sort((y-xi)/omega)

zp <- seq(zr[1],zr[2],l=100)
ghp <- gh_trans(zp,g,h)

zr1 <- max(zp[(u[1]>ghp)])
zr2 <- min(zp[max(u)<ghp])
zp <- seq(zr1,zr2,l=nbreaks)
zp <- sort(unique(zp))

ghp <- gh_trans(zp,g,h)
dzp <- diff(zp,lag=1)
dfghp <- diff(ghp,lag=1)
slop <- dfghp/dzp

loc <- .bincode(u,breaks=ghp)
uL <- ghp[loc]
zpL <- zp[loc]
zhat <-zpL+(u-uL)/dfghp[loc]*dzp[loc]
ghdq.hat <- fstar(zhat,g,h)
likelihood <- -sum(ghdq.hat)
value <- likelihood+n*log(omega)
obj2 <- fstar_hessian(zhat,xi,omega,g,h)
partial.hat <- obj2$gradient
}else{
value <- 10^9
partial.hat <-  rnorm(4)*100
}
return(partial.hat)
}

par0[2] <- log(par0[2])
obj1 <-optim(par0,fn,gr=grn,method="L-BFGS-B",lower = c(-Inf,-10,-Inf,0),upper=c(Inf,10,Inf,Inf),hessian=TRUE) 
par<- obj1$par
gradient <- grn(par)
par[2] <- exp(par[2])
value <- gh_likefun(y,par,nbreaks=nbreaks)
if(value==10^9) converge <- FALSE else converge <- TRUE
list(par=par,likelihood=value,converge=obj1$converge,gradient=gradient,hessian=obj1$hessian,message=obj1$message,niter=obj1$count[2])
}

###log-likelihood for normal data

likefun_normfunction<- function(y){
n <- length(y)
xihat <- mean(y)
omegahat <- sqrt(mean((y-mean(y))^2))
u <- (y-xihat)/omegahat
 sum(u^2)/2+n*log(omegahat)
}

###conducting test for g-h distribution########
##Warning: if the algorithm does not converge, you need change values for par0 (initial values for alternative estimation) and par0_null (initial values for null estimation)
gh_test<- function(y,par0=NULL,par0_null=NULL,nbreaks=nbreaks,Test="g=0"){
n <- length(y)
cutoff <- 1/n^2
if(is.null(par0))
par0 <- Initials(y)
par <- gh_likest_linear(y,par0,nbreaks=nbreaks)                        

if(Test=="g=0"){
if(is.null(par0_null))
par0_null <- Initials_gnull(y)
par_null <- gh_likest_linear_gnull(y,par0_null[-3],nbreaks=nbreaks)  
T<- 2*(par_null$likelihood-par$likelihood)   
pvalue= 1-pchisq(T,df=1)
}else if(Test=="h=0"){
if(is.null(par0_null))
par0_null <- Initials_hnull(y)
par_null <- gh_likest_linear_hnull(y,par0_null[-4],nbreaks=nbreaks)   
T<- 2*(par_null$likelihood-par$likelihood)   
pvalue= 1-pchibarsq(T,df=1,mix=0.5)                                        
}else if(Test=="g=0,h=0"){
T <- 2*(likefun_normfunction(y)-par$likelihood)
pvalue= 1-pchibarsq(T,df=2,mix=0.5)
}
list(STAT=T,pvalue=pvalue)
}


#######The approximated inverse g-h transformation
gh_likest_inv <- function(y,xi,omega,g,h,nbreaks=1000){
n <- length(y)
cutoff <- 1/n^2
alpha <- 0.1
nbs <- floor(nbreaks/3)
p <- seq(alpha,1-alpha,l=nbs)
bnd <- 10
zcutL <- max(bnd,abs(qnorm(cutoff)))
zcutU <- abs(qnorm(alpha))
zp <- seq(-zcutL,zcutL,l=nbreaks)
zp <- sort(unique(zp))
u <- (y-xi)/omega
ghp <- gh_trans(zp,g,h)
dfzp <- diff(zp,lag=1)
dfghp <- diff(ghp,lag=1)
slop <- dfghp/dfzp
loc <- .bincode(u,breaks=ghp)
uL <- ghp[loc]
zpL <- zp[loc]
zphat <-zpL+(u-uL)/dfghp[loc]*dfzp[loc] 
zphat 
}

####Regression model with g-h residuals
gh_likest_reg <- function(y,X=NULL,beta0=NULL,par0,nbreaks=1000){
n <- length(y)
ymin <- min(y)
ymax <- max(y)

cutoff <- 1/n^2
bnd <- 10
zcutL <- max(bnd,abs(qnorm(cutoff)))

fn <- function(par){
xi <- par[1]
omega <- exp(par[2])
g <- par[3]
h <- exp(par[4])
h <- h/(1+h)
y.temp  <- y
if(!is.null(X)){
beta <- as.matrix(par[-(1:4)])
y.temp <-  y-X%*%beta
ymin <- min(y.temp)
ymax <- max(y.temp)
}

tmp <- gh_trans(c(-zcutL,zcutL),g,h)*omega+xi

if(tmp[1]<ymin&&tmp[2]>ymax){
zphat <- gh_likest_inv(y.temp,xi,omega,g,h,nbreaks=nbreaks)
zphat <- as.matrix(zphat)
likelihood <- sum(zphat^2)/2+sum(log(abs(exp(g*zphat+h/2*zphat^2)+h*zphat/g*(exp(g*zphat)-1)*exp(h/2*zphat^2))))
value <- likelihood+n*log(omega)
rm(zphat)
}else{
value <- 10^9
}
value
}

par0[2] <- log(par0[2])
par0[4] <- par0[4]/(1-par0[4])
par0[4] <- log(par0[4])
if(!is.null(X)) par0 <- c(par0,beta0)

obj1 <-optim(par0,fn) 
par<- obj1$par
par[4] <- exp(par[4])
par[4] <- par[4]/(1+par[4])
par[2] <- exp(par[2])
value <- obj1$value
if(value==10^9) converge <- FALSE else converge <- TRUE
list(par=par,likelihood=value,converge=converge)
}