test3:

#' Estimation of Autoregressive Condtional Mean Models
#' 
#' Estimation of autoregressive conditionl mean models with exogeneous variables.
#' @param y time series of counts.
#' @param order the order of ACM model.
#' @param X matrix of exogenous variables.
#' @param cond.dist conditional distributions. "po" for Poisson, "nb" for negative binomial, "dp" for double Poisson.
#' @param ini initial parameter estimates designed for use in "nb" and "dp".
#' @return ACMx returns a list with components:
#' \item{data}{time series.}
#' \item{X}{matrix of exogenous variables.}
#' \item{estimates}{estimated values.}
#' \item{residuals}{residuals.}
#' \item{sresi}{standardized residuals.}
#' @examples
#' X=matrix(rnorm(100,100,1))
#' y=rpois(100,10)
#' ACMx(y,c(1,1),X,"po")
#' @export
"ACMx" <- function(y,order=c(1,1),X=NULL,cond.dist="po",ini=NULL){
beta=NULL; k=0; withX=FALSE; nT=length(y)
if(!is.null(X)){
withX=TRUE
if(!is.matrix(X))X=as.matrix(X)
T1=dim(X)[1]
if(nT > T1)nT=T1
if(nT < T1){T1=nT; X=X[1:T1,]}
##
k=dim(X)[2]
### Preliminary estimate of regression coefficients
m1=glm(y~X,family=poisson)
m11=summary(m1)
beta=m11$coefficients[-1,1]; se.beta=m11$coefficients[-1,2]
glmresi=m1$residuals
}
### obtain initial estimate of the constant
mm=glm(y[2:nT]~y[1:(nT-1)],family=poisson)
m22=summary(mm)
ome=m22$coefficients[1,1]; se.ome=m22$coefficients[1,2]
###
p=order[1]; q=order[2]; S=1.0*10^(-6); params=NULL; loB=NULL; upB=NULL
###
###### Poisson model
if(cond.dist=="po"){
if(withX){
params=c(params,beta=beta); loB=c(loB,beta=beta-2*abs(beta)); upB=c(upB,beta=beta+2*abs(beta))
}
params=c(params,omega=ome); loB=c(loB,omega=S); upB=c(upB,omega=ome+10*ome)
if(p > 0){
a1=rep(0.05,p); params=c(params,alpha=a1); loB=c(loB,alpha=rep(S,p)); upB=c(upB,alpha=rep(0.5,p))
}
if(q > 0){
b1=rep(0.5,q)
params=c(params,gamma=b1); loB=c(loB,gamma=rep(S,q)); upB=c(upB,gamma=rep(1-S,q))
}
##mm=optim(params,poX,method="L-BFGS-B",hessian=T,lower=loB,upper=upB)
#
cat("Initial estimates: ",params,"\n")
cat("loB: ",loB,"\n")
cat("upB: ",upB,"\n")
fit=nlminb(start = params, objective= poX, lower=loB, upper=upB, PCAxY=y, PCAxX=X, PCAxOrd=order,
control=list(rel.tol=1e-6))
#control=list(trace=3,rel.tol=1e-6))

epsilon = 0.0001 * fit$par
npar=length(params)
Hessian = matrix(0, ncol = npar, nrow = npar)
for (i in 1:npar) {
for (j in 1:npar) {
x1 = x2 = x3 = x4  = fit$par
x1[i] = x1[i] + epsilon[i]; x1[j] = x1[j] + epsilon[j]
x2[i] = x2[i] + epsilon[i]; x2[j] = x2[j] - epsilon[j]
x3[i] = x3[i] - epsilon[i]; x3[j] = x3[j] + epsilon[j]
x4[i] = x4[i] - epsilon[i]; x4[j] = x4[j] - epsilon[j]
Hessian[i, j] = (poX(x1,PCAxY=y,PCAxX=X,PCAxOrd=order)-poX(x2,PCAxY=y,PCAxX=X,PCAxOrd=order)
    -poX(x3,PCAxY=y,PCAxX=X,PCAxOrd=order)+poX(x4,PCAxY=y,PCAxX=X,PCAxOrd=order))/
(4*epsilon[i]*epsilon[j])
}
}
cat("Maximized log-likehood: ",-poX(fit$par,PCAxY=y,PCAxX=X,PCAxOrd=order),"\n")
# Step 6: Create and Print Summary Report:
se.coef = sqrt(diag(solve(Hessian)))
tval = fit$par/se.coef
matcoef = cbind(fit$par, se.coef, tval, 2*(1-pnorm(abs(tval))))
dimnames(matcoef) = list(names(tval), c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
cat("\nCoefficient(s):\n")
printCoefmat(matcoef, digits = 6, signif.stars = TRUE)
###### Compute the residuals
est=fit$par; ist=1
bb=rep(1,length(y)); y1=y
if(withX){beta=est[1:k]; bb=exp(X%*%matrix(beta,k,1))}#; y1=y/bb}
icnt=k+1
ome=est[icnt]
p=order[1]; q=order[2]; nT=length(y)
if(p > 0){a1=est[(icnt+1):(icnt+p)]; icnt=icnt+p
nobe=nT-p; rate=rep(ome,nobe); ist=p+1
for (i in 1:p){
rate=rate+a1[i]*y1[(ist-i):(nT-i)]
}
}
r1=ome
if(q > 0){g1=est[(icnt+1):(icnt+q)]
r1=filter(rate,g1,"r",init=rep(mean(y/bb),q))
}
rate=bb[ist:nT]*r1
resi=y[ist:nT]-rate
sresi=resi/sqrt(rate)
}
#### Negative binomial ######################################################
if(cond.dist=="nb"){
##
if(is.null(ini)){
if(withX){
params=c(params,beta=rep(0.4,k)); loB=c(loB,beta=rep(-3,k)); upB=c(upB,beta=rep(2,k))
}
params=c(params,omega=ome); loB=c(loB,omega=S); upB=c(upB,omega=ome+6*se.ome)
if(p > 0){
a1=rep(0.05,p); params=c(params,alpha=a1); loB=c(loB,alpha=rep(S,p)); upB=c(upB,alpha=rep(0.5,p))
}
if(q > 0){
b1=rep(0.7,q)
params=c(params,gamma=b1); loB=c(loB,gamma=rep(S,q)); upB=c(upB,gamma=rep(1-S,q))
}
## given initial estimates
}
else{
se.ini=abs(ini/10); jst= k
if(withX){
params=c(params,beta=ini[1:k]); loB=c(loB,beta=ini[1:k]-4*se.ini[1:k]); upB=c(upB,beta=ini[1:k]+4*se.ini[1:k])
}
jst=k+1
params=c(params,omega=ini[jst]); loB=c(loB,omega=ini[jst]-3*se.ini[jst]); upB=c(upB,omega=ini[jst]+3*se.ini[jst])
if(p > 0){a1=ini[(jst+1):(jst+p)]; params=c(params,alpha=a1); loB=c(loB,alpha=rep(S,p))
upB=c(upB,alpha=a1+3*se.ini[(jst+1):(jst+p)]); jst=jst+p
}
if(q > 0){
b1=ini[(jst+1):(jst+q)]
params=c(params,gamma=b1); loB=c(loB,gamma=rep(S,q)); upB=c(upB,gamma=rep(1-S,q))
}
}
### size of the negative binomial parameter
meanY=mean(y); varY=var(y); th1=meanY^2/(varY-meanY)*2; if(th1 < 0)th1=meanY*0.1
params=c(params,theta=th1); loB=c(loB,theta=0.5); upB=c(upB,theta=meanY*1.5)
#
cat("initial estimates: ",params,"\n")
cat("loB: ",loB,"\n")
cat("upB: ",upB,"\n")
fit=nlminb(start = params, objective= nbiX, lower=loB, upper=upB, PCAxY=y,PCAxX=X,PCAxOrd=order,
control=list(rel.tol=1e-6))
#control=list(trace=3,rel.tol=1e-6))

epsilon = 0.0001 * fit$par
npar=length(params)
Hessian = matrix(0, ncol = npar, nrow = npar)
for (i in 1:npar) {
for (j in 1:npar) {
x1 = x2 = x3 = x4  = fit$par
x1[i] = x1[i] + epsilon[i]; x1[j] = x1[j] + epsilon[j]
x2[i] = x2[i] + epsilon[i]; x2[j] = x2[j] - epsilon[j]
x3[i] = x3[i] - epsilon[i]; x3[j] = x3[j] + epsilon[j]
x4[i] = x4[i] - epsilon[i]; x4[j] = x4[j] - epsilon[j]
Hessian[i, j] = (nbiX(x1,PCAxY=y,PCAxX=X,PCAxOrd=order)-nbiX(x2,PCAxY=y,PCAxX=X,PCAxOrd=order)
    -nbiX(x3,PCAxY=y,PCAxX=X,PCAxOrd=order)+nbiX(x4,PCAxY=y,PCAxX=X,PCAxOrd=order))/
(4*epsilon[i]*epsilon[j])
}
}
cat("Maximized log-likehood: ",-nbiX(fit$par,PCAxY=y,PCAxX=X,PCAxOrd=order),"\n")
# Step 6: Create and Print Summary Report:
se.coef = sqrt(diag(solve(Hessian)))
tval = fit$par/se.coef
matcoef = cbind(fit$par, se.coef, tval, 2*(1-pnorm(abs(tval))))
dimnames(matcoef) = list(names(tval), c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
cat("\nCoefficient(s):\n")
printCoefmat(matcoef, digits = 6, signif.stars = TRUE)
#### Compute the residuals
est=fit$par
bb=rep(1,length(y)); y1=y
if(withX){beta=est[1:k]; bb=exp(X%*%matrix(beta,k,1))}#; y1=y/bb}
icnt=k+1
ome=est[icnt]
p=order[1]; q=order[2]; nT=length(y)
if(p > 0){a1=est[(icnt+1):(icnt+p)]; icnt=icnt+p
nobe=nT-p; rate=rep(ome,nobe); ist=p+1
for (i in 1:p){
rate=rate+a1[i]*y1[(ist-i):(nT-i)]
}
}
if(q > 0){g1=est[(icnt+1):(icnt+q)]; icnt=icnt+q
r1=filter(rate,g1,"r",init=rep(mean(y/bb),q))
}
rate=bb[ist:nT]*r1
resi=y[ist:nT]-rate
theta=est[icnt+1]
pb=theta/(theta+rate)
v1=theta*(1-pb)/(pb^2)
sresi=resi/sqrt(v1)

}
### Double Poisson model ###################################
if(cond.dist=="dp"){
##
if(withX){
params=c(params,beta=beta); loB=c(loB,beta=beta-2*abs(beta)); upB=c(upB,beta=beta+2*abs(beta))
}
params=c(params,omega=ome*0.3); loB=c(loB,omega=S); upB=c(upB,omega=ome+4*se.ome)
if(p > 0){
a1=rep(0.05,p); params=c(params,alpha=a1); loB=c(loB,alpha=rep(S,p)); upB=c(upB,alpha=rep(1-S,p))
}
if(q > 0){
b1=rep(0.6,q)
params=c(params,gamma=b1); loB=c(loB,gamma=rep(S,q)); upB=c(upB,gamma=rep(1-S,q))
}
### size of the negative binomial parameter
meanY=mean(y); varY = var(y); t1=meanY/varY
params=c(params,theta=t1); loB=c(loB,theta=S); upB=c(upB,theta=2)
#
cat("initial estimates: ",params,"\n")
cat("loB: ",loB,"\n")
cat("upB: ",upB,"\n")
fit=nlminb(start = params, objective= dpX, lower=loB, upper=upB,PCAxY=y,PCAxX=X,PCAxOrd=order,
control=list(rel.tol=1e-6))
#control=list(trace=3,rel.tol=1e-6))

epsilon = 0.0001 * fit$par
npar=length(params)
Hessian = matrix(0, ncol = npar, nrow = npar)
for (i in 1:npar) {
for (j in 1:npar) {
x1 = x2 = x3 = x4  = fit$par
x1[i] = x1[i] + epsilon[i]; x1[j] = x1[j] + epsilon[j]
x2[i] = x2[i] + epsilon[i]; x2[j] = x2[j] - epsilon[j]
x3[i] = x3[i] - epsilon[i]; x3[j] = x3[j] + epsilon[j]
x4[i] = x4[i] - epsilon[i]; x4[j] = x4[j] - epsilon[j]
Hessian[i, j] = (dpX(x1,PCAxY=y,PCAxX=X,PCAxOrd=order)-dpX(x2,PCAxY=y,PCAxX=X,PCAxOrd=order)
   -dpX(x3,PCAxY=y,PCAxX=X,PCAxOrd=order)+dpX(x4,PCAxY=y,PCAxX=X,PCAxOrd=order))/
(4*epsilon[i]*epsilon[j])
}
}
cat("Maximized log-likehood: ",-dpX(fit$par,PCAxY=y,PCAxX=X,PCAxOrd=order),"\n")
# Step 6: Create and Print Summary Report:
se.coef = sqrt(diag(solve(Hessian)))
tval = fit$par/se.coef
matcoef = cbind(fit$par, se.coef, tval, 2*(1-pnorm(abs(tval))))
dimnames(matcoef) = list(names(tval), c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
cat("\nCoefficient(s):\n")
printCoefmat(matcoef, digits = 6, signif.stars = TRUE)
#### Compute the residuals
est=fit$par
bb=rep(1,length(y)); y1=y
if(withX){beta=est[1:k]; bb=exp(X%*%matrix(beta,k,1))}#; y1=y/bb}
icnt=k+1
ome=est[icnt]
p=order[1]; q=order[2]; nT=length(y)
if(p > 0){a1=est[(icnt+1):(icnt+p)]; icnt=icnt+p
nobe=nT-p; rate=rep(ome,nobe); ist=p+1
for (i in 1:p){
rate=rate+a1[i]*y1[(ist-i):(nT-i)]
}
}
if(q > 0){g1=est[(icnt+1):(icnt+q)]; icnt=icnt+q
r1=filter(rate,g1,"r",init=rep(mean(y/bb),q))
}
theta=est[icnt+1]
rate=bb[ist:nT]*r1
resi=y[ist:nT]-rate
v1=rate/theta
sresi=resi/sqrt(rate)
}

#### end of the program
ACMx <- list(data=y,X=X,estimates=est,residuals=resi,sresi=sresi)
}


#' @export
"nbiX" <- function(par,PCAxY,PCAxX,PCAxOrd){
## compute the log-likelihood function of negative binomial distribution
y <- PCAxY; X <- PCAxX; order <- PCAxOrd
withX = F; k = 0
#
if(!is.null(X)){
withX=T; k=dim(X)[2]
}
bb=rep(1,length(y)); y1=y
if(withX){beta=par[1:k]; bb=exp(X%*%matrix(beta,k,1))}#; y1=y/bb}
icnt=k+1
ome=par[icnt]
p=order[1]; q=order[2]; nT=length(y)
if(p > 0){a1=par[(icnt+1):(icnt+p)]; icnt=icnt+p
nobe=nT-p; rate=rep(ome,nobe); ist=p+1
for (i in 1:p){
rate=rate+a1[i]*y1[(ist-i):(nT-i)]
}
}
if(q > 0){g1=par[(icnt+1):(icnt+q)]; icnt=icnt+q
r1=filter(rate,g1,"r",init=rep(mean(y/bb),q))
}
theta=par[icnt+1]
rate=bb[ist:nT]*r1
pb=theta/(theta+rate)
lnNBi=dnbinom(y[ist:nT],size=theta,prob=pb,log=T)
nbiX <- -sum(lnNBi)
}


#' @export
"poX" <- function(par,PCAxY,PCAxX,PCAxOrd){
y <- PCAxY; X <- PCAxX; order <- PCAxOrd
withX = F; k=0
##
if(!is.null(X)){
withX=T; k=dim(X)[2]
}
bb=rep(1,length(y)); y1=y
if(withX){beta=par[1:k]; bb=exp(X%*%matrix(beta,k,1))}#; y1=y/bb}
icnt=k+1; ist=1; nT=length(y); nobe=nT
ome=par[icnt]; rate=rep(ome,nobe)
p=order[1]; q=order[2]
if(p > 0){a1=par[(icnt+1):(icnt+p)]; icnt=icnt+p
nobe=nT-p; rate=rep(ome,nobe); ist=p+1
for (i in 1:p){
rate=rate+a1[i]*y1[(ist-i):(nT-i)]
}
}
r1=rate
if(q > 0){g1=par[(icnt+1):(icnt+q)]
r1=filter(rate,g1,"r",init=rep(mean(y/bb),q))
}

rate=bb[ist:nT]*r1
lnPoi=dpois(y[ist:nT],rate,log=T)
poX <- -sum(lnPoi)
}




#' @export
"dpX" <- function(par,PCAxY,PCAxX,PCAxOrd){
# compute the log-likelihood function of a double Poisson distribution
y <- PCAxY; X <- PCAxX; order <- PCAxOrd
withX=F; k = 0
#
if(!is.null(X)){
withX=TRUE; k=dim(X)[2]
}
bb=rep(1,length(y)); y1=y
if(withX){beta=par[1:k]; bb=exp(X%*%matrix(beta,k,1))}#; y1=y/bb}
icnt=k+1
ome=par[icnt]
p=order[1]; q=order[2]; nT=length(y)
if(p > 0){a1=par[(icnt+1):(icnt+p)]; icnt=icnt+p
nobe=nT-p; rate=rep(ome,nobe); ist=p+1
for (i in 1:p){
rate=rate+a1[i]*y1[(ist-i):(nT-i)]
}
}
#plot(rate,type='l')

if(q > 0){g1=par[(icnt+1):(icnt+q)]; icnt=icnt+q
r1=filter(rate,g1,"r",init=rep(mean(y/bb),q))
}
theta=par[icnt+1]
rate=bb[ist:nT]*r1
yy=y[ist:nT]
rate1 = rate*theta
cinv=1+(1-theta)/(12*rate1)*(1+1/rate1)
lcnt=-log(cinv)
lpd=lcnt+0.5*log(theta)-rate1
idx=c(1:nobe)[yy > 0]
tp=length(idx)
d1=apply(matrix(yy[idx],tp,1),1,factorialOwn)
lpd[idx]=lpd[idx]-d1-yy[idx]+yy[idx]*log(yy[idx])+theta*yy[idx]*(1+log(rate[idx])-log(yy[idx])) 
#plot(lpd,type='l')
#cat("neg-like: ",-sum(lpd),"\n")

dpX <- -sum(lpd)
}


#' @export
"factorialOwn" <- function(n,log=T){
x=c(1:n)
if(log){
x=log(x)
y=cumsum(x)
}
else{
y=cumprod(x)
}
y[n]
}