condivreg<-function(y1,y2,x1,Z,b_0,alpha){
yall=as.matrix(cbind(y1,y2,x1,Z)) #combine all data into a single matrix
newdata=na.omit(yall) #find all nas in the data and omit them
nall=dim(newdata) #calculates the rows and columns of the new data matrix
nx=dim(as.matrix(x1)) #calculates the dimension of the exogenous x variables
nz=dim(as.matrix(Z)) #calculates the dimension of the IV matrix
k1=nx[2] #number of exognous x variables
k=nz[2] #number of IVs
n=nall[1] #total sample size
y1=newdata[,1] #reform dependent variable
y2=newdata[,2] #reform endogenous variable
ind1=seq(3,3-1+k1) #generate sequence to extract x variables
ind2=seq(3+k1,nall[2]) #generate sequence to extract IVs
x1=cbind(as.matrix((rep(1,n))),as.matrix(newdata[,ind1])) #combine vector of ones with x matrix
X=x1 #define X matrix
Z=as.matrix(newdata[,ind2]) #matrix of IVs
Z_=Z #non-residual based matrix of IVs
Px = X%*%solve(t(X)%*%X)%*%t(X)
Mx = diag(n) - Px
Z = Mx%*%Z
Pz = Z%*%solve(t(Z)%*%Z)%*%t(Z)
Mz = diag(n) - Pz
Y = cbind(y1,y2)
Yt=Mx%*%Y
y1t=Yt[,1]
y2t=Yt[,2]
k2sls=1
b2sls=solve(t(y2t)%*%(diag(n)-k2sls*Mz)%*%y2t)%*%(t(y2t)%*%(diag(n)-k2sls*Mz)%*%y1t) #2sls estimate for beta
A=(t(Yt)%*%Yt)
B=t(Yt)%*%Mz%*%Yt
eigs=eigen(A%*%solve(B))
kliml=eigs$values[2]
bliml=solve(t(y2t)%*%(diag(n)-kliml*Mz)%*%y2t)%*%(t(y2t)%*%(diag(n)-kliml*Mz)%*%y1t) #liml estimate for beta
V = Y-Pz%*%Y-Px%*%Y
omega = (t(V)%*%V)/(n-k-k1-1)
b0=t(cbind(1,-b_0))
a0 = t(cbind(b_0,1))
ZZ=t(Z)%*%Z

#check to see how many IVs we have and then calculate the ZZ^(-1/2) matrix
if (k==1) {ZZ2=ZZ^(-1/2)}  else
if (k>1){z.eig <- eigen(ZZ);z.sqrt <- z.eig$vectors %*% diag(sqrt(z.eig$values)) %*% solve(z.eig$vectors);
ZZ2=solve(z.sqrt)}



#Starting computation of Moreira (2003) conditional likelihood ratio (CLR) test
T = (ZZ2)%*%t(Z)%*%Y%*%(solve(omega))%*%a0%*%((t(a0)%*%(solve(omega))%*%a0)^(-1/2))
S1= (ZZ2)%*%t(Z)%*%Y%*%b0%*%((t(b0)%*%omega%*%b0)^(-1/2))
Qt = t(T)%*%T
lnQtk = log(Qt/k)
Qs = t(S1)%*%S1
Qt= t(T)%*%T
Qst= t(S1)%*%T
Qts= t(T)%*%S1
CLR= 1/2*(Qs-Qt+((Qs-Qt)^2+(4*(Qst^2)))^(1/2));
x=cbind(X,Z_);
SSRur=sum(((y2-x%*%solve(t(x)%*%x)%*%t(x)%*%y2)^2))
SSRr=sum((y2-X%*%solve(t(X)%*%X)%*%t(X)%*%y2)^2)
dimallx=dim(x)
m=CLR
qt=Qt

F=((SSRr-SSRur)/k)/(SSRur/(n-dimallx[2]-1)) #F-stat on excluded IVs

LB=-b2sls*50 #lower bound for null evaluation
UB=b2sls*50 #upper bound for null evaluation
b00=seq(LB,UB,b2sls/5) #generate sequence of nulls
nb=dim(as.matrix(b00)) 
nb=nb[1]
CLRf<-vector('list',nb) #pre-allocate space for CLR evaluated on set of nulls of length nb

#loop for evaluation of CLR across nulls
for(i in 1:nb){
b01=t(cbind(1,-b00[i]))
a01 = t(cbind(b00[i],1))
T2 = (ZZ2)%*%t(Z)%*%Y%*%(solve(omega))%*%a01%*%((t(a01)%*%(solve(omega))%*%a01)^(-1/2))
S12= (ZZ2)%*%t(Z)%*%Y%*%b01%*%((t(b01)%*%omega%*%b01)^(-1/2))
Qt2 = t(T2)%*%T2
Qs2 = t(S12)%*%S12
Qt2= t(T2)%*%T2
Qst2= t(S12)%*%T2
Qts2= t(T2)%*%S12
CLR2=1/2*(Qs2-Qt2+((Qs2-Qt2)^2+(4*(Qst2^2)))^(1/2))
CLRf[[i]]=1/2*(Qs2-Qt2+((Qs2-Qt2)^2+(4*(Qst2^2)))^(1/2))
}



if (k==1) {cv=qchisq(1-alpha,k);pvalue=1-pchisq(m,k)}  else 
if (k>1){K4=gamma(k/2)/((pi^(1/2)*gamma((k-1)/2)));integrand <- function(x) {pchisq((qt+m)/(1+((qt*x^2)/m)),k)*(1-x^2)^((k-3)/2)}
;p1=integrate(integrand, lower = 0, upper = 1);
pvalue=1-2*K4*p1$value;
f <- function (m,a) {p1=integrate(function(x) {pchisq((qt+m)/(1+((qt*x^2)/m)),k)*(1-x^2)^((k-3)/2)}, lower = 0, upper = 1);a-1+2*K4*p1$value};
xmin <- uniroot(f, c(1, 1000), tol = 0.0001, a = alpha);cv=xmin$root}

#Finding confidence intervals
f1<-function(xx){b01=t(cbind(1,-xx));a01 = t(cbind(xx,1));T2 = (ZZ2)%*%t(Z)%*%Y%*%(solve(omega))%*%a01%*%((t(a01)%*%(solve(omega))%*%a01)^(-1/2));S12= (ZZ2)%*%t(Z)%*%Y%*%b01%*%((t(b01)%*%omega%*%b01)^(-1/2));
Qt2 = t(T2)%*%T2;
Qs2 = t(S12)%*%S12;
Qt2= t(T2)%*%T2;
Qst2= t(S12)%*%T2;
Qts2= t(T2)%*%S12;
point=cv-1/2*(Qs2-Qt2+((Qs2-Qt2)^2+(4*(Qst2^2)))^(1/2));}

lb=LB
ub=LB+b2sls/5
ci<-vector('list',1000)
for(ii in 1:10000){checkme<-f1(lb)*f1(ub); if (checkme<0){ci[[ii]]=uniroot(f1, c(lb, ub),tol=.00000001)$root};if (ub>UB){break};lb=ub;ub=ub+b2sls/5}


ci=ci[-(which(sapply(ci,is.null),arr.ind=TRUE))]
ci=do.call(cbind, ci)

if(CLRf[1]<cv && CLRf[nb]<cv && is.numeric(ci)){llb=-Inf;uub=Inf;ci=cbind(-Inf,ci,Inf)}else
if(CLRf[1]>cv && CLRf[nb]>cv){ci=ci}

if (is.numeric(ci)){ci=ci} else ci=cbind(-Inf,Inf)
plot(b00,CLRf,type='l')
lines(b00,rep(cv,nb),col=2)
return(list(TSLS=b2sls,LIML=bliml,CLR=CLR,CI=ci,pvalue=pvalue,First_Stage_Fstat=F,CV=cv,Qt=Qt,N=n))

}