"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)
}