#####____________________________________________________________________________#####
#####                                                                            #####
#####     Restricted Likelihood Ratio Test for AR(1)Predictive Regression Model  ##### 
#####     from "Bias Reduction and Likelihood Based Almost-Exactly Sized         #####
#####     Hypothesis Testing in Predictive Regressions using the Restricted      #####
#####     Likelihood" W. Chen and R. Deo (2009)                                  #####
#####                        March 6, 2009                                       #####
#####____________________________________________________________________________#####
#####                                                                            #####
#####          Likelihood funtion assumes stationarity of process                ##### 
#####          ( instead of fixed initial value)                                 #####
#####____________________________________________________________________________#####
#####                                                                            #####


PredReg<-function(X,Y){
	n<-length(X)-1;n<<-n
	yt<-Y;xt<-X
	
	##########################################
	####	demeaned y_t, x_t, x_t-1      ####
	##########################################
	
	demean<-function(x) {x-mean(x)}
	yc<-demean(yt[-1]); 	### demeaned y_t ###
	xc<-cbind(demean(xt[-1]),demean(xt[-(n+1)]));	### demeaned x_t,x_t-1  ###
	yc<<-yc; xc<<-xc; 	
	ini<-lsfit(xc[,2],xc[,1],intercept=F)$coef;ini<<-c(ini,ini);ini<-rep(min(ini,.99),2)

	############################################
	####         REML estimates for         ####
	#### phi, gamma, alpha, esig2,vsig2 #### 
	############################################

		
	outRes<-res(xc,yc)		### restricted, 1 by 5	 ###		
	outUnRes<-unRes(xc,yc)	### unrestricted, 1 by 5 ###
	
	outNames<-c("phi", "gamma", "alpha","vsig2", "esig2")
	names(outRes)<-outNames; names(outUnRes)<-outNames;
	
	############################################
	####     Computing Liklihood ratio      ####		
	############################################
	LRT<-LR(outRes,outUnRes)	### Liklihood ratio with two estimates  ###
	
	############################################
	####            Print Output            ####		
	############################################
	
	beta<-outUnRes["gamma"]+outUnRes["phi"]*outUnRes["alpha"]  ### get beta from gamma,phi,alpha ###
	
	outUnRes["gamma"]<-beta;names(outUnRes)<-c("phi", "beta", "alpha","vsig2", "esig2")
	
	cat("\n REML Estimates:\n")
	print(round(outUnRes,4))
	LRT<-c(LRT,1-pchisq(LRT,1))
	names(LRT)<-c("Statistic","p-value")
	cat("\n REML LRT:\n")
	print(round(LRT,4))
	
}



#################################################
######		unrestricted estimates         ######
#################################################

unRes<-function(xc,yc){

	####    getting theta and esig2 by OLS 	####
	unResOut<-lsfit(xc,yc,intercept=F)
	unRestheta<-unResOut$coef
	esig2<-1/(n-1)*sum(unResOut$residuals^2)
	
	####   getting alpha and vsig2 by MLE   ####
	#unResLout<-nlminb(ini[2],unResL,lower=-.99999,upper=.99999)
	unResLout<-nlminb(ini[2],unResL,lower=-.99999,upper=1)	
	unResAlpha<-unResLout$par	
	
	vsig2<-1/n*Qx(unResAlpha)

	c(unRestheta,unResAlpha,vsig2,esig2)
}


###### likelihood function ######

unResL<-function(a){	
	
	n*log(Qx(a))-log((1+a)/((n-1)*(1-a)+2))
}



Qx<-function(alpha){
	muhat<-((1-alpha)*sum(xt[2:n])+xt[1]+xt[(n+1)])/((n-1)*(1-alpha)+2)
	x<-xt-muhat
	Bx<-c(x[1]*sqrt(1-alpha^2),x[-1]-alpha*x[-(n+1)])
	sum(Bx^2)
}
####################################################
######		restricted estimates           ######
#####################################################


res<-function(xc,yc){
	#resLout<-nlminb(ini[1],resL,lower=-.99999,upper=.99999)
	resLout<-nlminb(ini[1],resL,lower=-.99999,upper=1)
	resAlpha<-resLout$par
	z<-xc[,1]-resAlpha*xc[,2]
	phihat<-sum(yc*z)/sum(z^2)
	vsig2<-1/n*Qx(resAlpha)
	esig2<-1/(n-1)*R(resAlpha)
	c(phihat,-phihat*resAlpha,resAlpha,vsig2,esig2)
}

###### restricted likelihood function ######

resL<-function(a){
	unResL(a)+(n-1)*log(R(a))
}


###### extra term of RL from URL      ######
R<-function(alpha){
	zt<-xc[,1]-alpha*xc[,2]	  
	sum(yc^2)-sum(yc*zt)^2/sum(zt^2)
}

######################################
####### Likelihood ratio test ########
######################################
####### names(est)<-c("phi", "gamma", "alpha","vsig2", "esig2") ######

LR<-function(est0,est1){
	-2*L(est0)+2*L(est1)
}

L<-function(est){

	a1<--(n-1)/2*log(est["esig2"])
	a2<--1/(2*est["esig2"])*sum((yc-est["phi"]*xc[,1]-est["gamma"]*xc[,2])^2)
	a3<--n/2*log(est["vsig2"])
	a4<-1/2*log((1+est["alpha"])/((n-1)*(1-est["alpha"])+2))
	a5<--1/(2*est["vsig2"])*Qx(est["alpha"])
	a1+a2+a3+a4+a5
}