"tsls" <-

function(x, z, y, int = T)

{

# two stage least squares

if(int){

       x <- cbind(1, x)

       z <- cbind(1, z)

      }

xhat <- lm(x~z-1)$fitted.values

R <- lm(y ~ xhat -1)

R$residuals <- c(y - x %*% R$coef)

return(R)

}