#
#
#R Code for estimating ZIOP and ZIOPC models
#
#BEB 4/27/2011
#
#
#

#clear memory 
rm( list=ls() ) 

#load two necessary packages
library(mvtnorm)
library(MASS)

#set the seed
set.seed(2)

##################################################################
###########Create a dataset for use in the models below###########
##################################################################

#set the number of observations
  n<-1000

#define true gammas and betas
  G<-rbind(1, -0.25, -1)   
  G<-as.matrix(G)    
  cutpoints<-rbind(-100, 0, 4.5, 5.5, 100)
  B<-rbind(1,1)
  B<-as.matrix(B)

#create some covariates
  x0<-matrix(rep(1,n))                                                           
  x1<-matrix(log(runif(n, min=0, max=100)),n,1)
  x2star<-matrix(runif(n, min=0, max=1), n, 1)      
  x2<-ifelse(x2star>.25,1,0)   
  Z<-cbind(x0,x1,x2)
  Z<-as.matrix(Z)  
  X<-cbind(x1,x2)
  X<-as.matrix(X)

#create a zero inflated dependent variable that
#ranges from 0 to 3, with no correlated disturbances
  phi<-pnorm(Z%*%G) 
  vcv<-matrix(c(1,0,0,1), nrow=2, ncol=2)
  e<-mvrnorm(n, mu=c(0,0), Sigma=vcv)
  ystar<-X%*%B+e[,1]
  y<-matrix(0,n,1)
  for(i in 2:4){y[ystar>cutpoints[i]]<-y[ystar>cutpoints[i]]+1}
  yzero<-matrix(0,n,1)
  error<--1*e[,2]
  flag<-error<qnorm(phi)
  yzero[flag]<-error[flag]
  flag<-yzero==0
  y[flag]<-yzero[flag]

#save completed dataset
  dataset<-cbind(y,x1,x2) 
  colnames(dataset) <- c("y", "x1", "x2")
  dataset<-as.data.frame(dataset)



##################################################################
############################ZIOP Model############################
##################################################################


#This program estimates the ZIOP loglikelihood function

ZIOP<- function(est,Y,X,Z,data) {					      
	n=nrow(data)							      					  
	llik <- matrix(0, nrow=n, ncol = 1)
	y.cat<-nlevels(as.factor(Y))
	y0<-sort(unique(Y))
	V<-matrix(NA,nrow(data),y.cat)
	for(k in 1:y.cat){
		V[,k]<-Y==y0[k]
		}
	tau<-rep(1,max(Y))
	for (i in 1:max(Y)){
		tau[i]<-ifelse(i==1,est[i],tau[i-1]+exp(est[i]))
		}
	gamma<-est[(y.cat):(y.cat+ncol(Z)-1)]
	beta<-est[(y.cat+ncol(Z)):length(est)]
	ZG<-Z%*%gamma
	XB<-X%*%beta
	
	cprobs<-matrix(nrow=length(XB),ncol=y.cat)
	probs<-matrix(nrow=n,ncol=y.cat)
	     
	for(j in 1:(y.cat-1)){
		
		cprobs[,j]<-pnorm(tau[j]-XB)
		
		}
	
	probs[,y.cat]<-(1-cprobs[,(y.cat-1)])*pnorm(ZG) 
	probs[,1]<-(cprobs[,1])*pnorm(ZG)+(1-pnorm(ZG))	
	
	for(j in 2:(y.cat-1)){			
		
		probs[,j]<-(cprobs[,j]-cprobs[,(j-1)])*(pnorm(ZG))
		
		}
	
	llik<--1*sum( log(probs[V]) )
	
	return(llik)
	
	}



#set starting parameters
est<-rbind(.1,.1,.1,.1,.1,.1,.1,.1)

#set data, Y, X, and Z
data<-dataset
Y<-dataset$y
X<-cbind(dataset$x1,dataset$x2)
Z<-cbind(1, dataset$x1, dataset$x2)

#optimize
output.ZIOP<-nlm(f=ZIOP,  p=est, Y=Y,X=X,Z=Z ,iterlim=500, data=dataset, hessian=TRUE)
print(output.ZIOP)
vcv<-solve(output.ZIOP$hessian)
vcv 







####################################################################
############################ZIOPC Model#############################
###################################################################

#

#This program estimates the ZIOPC loglikelihood function

ZIOPC<- function(est,Y,X,Z,data) {					      
	n=nrow(data)							      					  
	llik <- matrix(0, nrow=n, ncol = 1)
	y.cat<-nlevels(as.factor(Y))
	y0<-sort(unique(Y))
	V<-matrix(NA,nrow(data),y.cat)
	for(k in 1:y.cat){
		V[,k]<-Y==y0[k]
		}
	tau<-rep(1,max(Y))
	for (i in 1:max(Y)){
		tau[i]<-ifelse(i==1,est[i],tau[i-1]+exp(est[i]))
		}

	gamma<-est[((y.cat)):(y.cat+ncol(Z)-1)]
	beta<-est[(y.cat+ncol(Z)):(length(est)-1)]
	rho<-est[length(est)]
	
	sigma<-matrix(c(1,rho,rho,1), nrow=2, ncol=2)
	negsigma<-matrix(c(1,-rho,-rho,1), nrow=2, ncol=2)
	
	
	
	cprobs<-matrix(nrow=n,ncol=y.cat)
	probs<-matrix(nrow=n,ncol=y.cat)
       
	for(j in 1:(y.cat-1)){
		
		for(i in 1:n){
		
			ZG<-Z[i,]%*%gamma
			XB<-X[i,]%*%beta
		
			cprobs[i,j]<-pmvnorm(lower=c(-Inf,-Inf), upper=c(ZG,tau[j]-XB),mean=c(0,0),sigma=negsigma)	
		
			}	
	
		}	
	
	for(i in 1:n){
	
		ZG<-Z[i,]%*%gamma
		XB<-X[i,]%*%beta
	
		probs[i,1]<-pmvnorm(lower=c(-Inf,-Inf), upper=c(ZG,tau[1]-XB),mean=c(0,0),sigma=negsigma)+(1-pnorm(ZG))
	
		}
	
	for(i in 1:n){	
		
		ZG<-Z[i,]%*%gamma
		XB<-X[i,]%*%beta
		
		probs[i,y.cat]<-(pmvnorm(lower=c(-Inf,-Inf), upper=c(ZG,(XB-tau[y.cat-1])),mean=c(0,0),sigma=sigma))  
		
		}
	
	for(j in 2:(y.cat-1)){			
		
		probs[,j]<-(cprobs[,j]-cprobs[,(j-1)])
		
		}
	
	llik<--1*sum( log(probs[V]) )
	
	return(llik)
	
	}


#set starting parameters
est<-rbind(.1,.1,.1,.1,.1,.1,.1,.1,.1)

#set data, Y, X, and Z
data<-dataset
Y<-dataset$y
X<-cbind(dataset$x1,dataset$x2)
Z<-cbind(1, dataset$x1, dataset$x2)

#optimize
output.ZIOPC<-nlm(f=ZIOPC,  p=est, Y=Y,X=X,Z=Z ,iterlim=500, data=dataset, hessian=TRUE)
print(output.ZIOPC)
vcv<-solve(output.ZIOPC$hessian)
vcv 





####################################################################
###################Alternate ZIOPC Model############################
###################################################################

#

#This program estimates a slightly faster version of the ZIOPC model

#For this version, one must set the number of Y categories manually
#(here the categories are set to 0-3, to correspond to the four categories of our Y).

ZIOPC <- function(est,Y,X,Z,data) {					      
	n=nrow(data)							      							      
	llik <- matrix(0, nrow=n, ncol = 1)
	y.cat<-nlevels(as.factor(Y))
	tau<-rep(1,max(Y))
	for (i in 1:max(Y)){
		tau[i]<-ifelse(i==1,est[i],tau[i-1]+exp(est[i]))
		}
	gamma<-est[((y.cat)):(y.cat+ncol(Z)-1)]
	beta<-est[(y.cat+ncol(Z)):(length(est)-1)]
	rho<-est[length(est)]
	for(i in 1:n){
		ZG<-Z[i,]%*%gamma
		XB<-X[i,]%*%beta
		sigma<-matrix(c(1,rho,rho,1), nrow=2, ncol=2)
		negsigma<-matrix(c(1,-rho,-rho,1), nrow=2, ncol=2)
		if(identical(Y[i],0)){llik[i]<-log((1-pnorm(ZG))+ pmvnorm(lower=c(-Inf,-Inf), upper=c(ZG,(tau[1]-XB)),mean=c(0,0),sigma=negsigma))}
		else if(identical(Y[i],1)){llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(ZG,(tau[2]-XB)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(ZG,(tau[1]-XB)),mean=c(0,0),sigma=negsigma))}
		else if(identical(Y[i],2)){llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(ZG,(tau[3]-XB)),mean=c(0,0),sigma=negsigma)-pmvnorm(lower=c(-Inf,-Inf), upper=c(ZG,(tau[2]-XB)),mean=c(0,0),sigma=negsigma))}
		else if(identical(Y[i],3)){llik[i]<-log(pmvnorm(lower=c(-Inf,-Inf), upper=c(ZG,(XB-tau[3])),mean=c(0,0),sigma=sigma))}     
		}
	llik<--1*sum(llik)					      
	return(llik)
	}								      



#set starting parameters
est<-rbind(.1,.1,.1,.1,.1,.1,.1,.1,.1)

#set data, Y, X, and Z
data<-dataset
Y<-dataset$y
X<-cbind(dataset$x1,dataset$x2)
Z<-cbind(1, dataset$x1, dataset$x2)

#optimize
output.ZIOPC<-nlm(f=ZIOPC,  p=est, Y=Y,X=X,Z=Z ,iterlim=500, data=dataset, hessian=TRUE)
print(output.ZIOPC)
vcv<-solve(output.ZIOPC$hessian)
vcv