tsls<-function(y1,y2,x,z){
yall=as.matrix(cbind(y1,y2,x,z)) #combine all data into a single matrix
newdata=na.omit(yall) #find all nas in the data and omit them
dall=dim(newdata) #calculates the rows and columns of the new data matrix
n=dall[1]
output=lm(y1~x+y2)
nm=names(output$coefficients)
x=as.matrix(x)
y1=as.matrix(y1)
z=as.matrix(z)
y2=as.matrix(y2)
kx=dim(x)
kx=kx[2]
kz=dim(z)
kz=kz[2]
outputur=lm(y2~x+z)
outputr=lm(y2~x)
uhatur=residuals(outputur)
uhatr=residuals(outputr)
SSRur=sum(uhatur^2)
SSRr=sum(uhatr^2)
F=((SSRr-SSRur)/kz)/(SSRur/(n-(kx+kz)-1))
F=as.matrix(F)
x=as.matrix(cbind(rep(1,n),x))
Z=as.matrix(cbind(x,z))
X=as.matrix(cbind(x,y2))
y2hat=fitted(outputur)
Xhat=as.matrix(cbind(x,y2hat))
B=solve(t(Xhat)%*%X)%*%t(Xhat)%*%y1
uhat=y1-X%*%B
sigma2=sum(uhat^2)/(n-kx-2)
VarCov=solve(t(Xhat)%*%Xhat)*sigma2
SE=as.matrix(diag(VarCov^(1/2)))
Tstats=B/SE
B<-cbind(B,SE,Tstats)
colnames(B) <- c("Coefficients","Standard-Errors","T-stats")
colnames(F) <- c("First-stage F-stat")
rownames(F) <- c("F-stat")
options(scipen=999)
rownames(B) <- c(nm)
return(list(uhat=uhat,B=B,F=F))
}