ivregress = function(second, first, data){

s.names = all.vars(second)

f.names  = all.vars(first)

data.names = names(data)

all.names = c(s.names,f.names)

resp   = s.names[1]

endog  = f.names[1]

inst   = f.names[-1]

explan = s.names[-1]

exog   = explan[explan!=endog]

exog.f = paste(exog,collapse="+")

inst.f = paste(inst, collapse="+")

RHS    = paste(exog.f, inst.f, sep="+")

first.form = as.formula( paste(endog, "~", RHS))

first.lm = lm(first.form, data)

ftest    = linearHypothesis(first.lm, inst, rep(0,length(inst)))

x.hat    = fitted(first.lm)

data2    = cbind(data,x.hat)

iname    = paste(endog,".hat",sep="")

names(data2) = c(names(data), iname)

data2.names = names(data2)

RHS2 = paste(exog.f,iname,sep="+")

second.form = as.formula(paste(resp, "~", RHS2))

second.lm = lm(second.form, data2)

Xmat  = data2[c(exog,endog)]

Xmat2  = data2[c(exog,iname)]

z = summary(second.lm)

X     = as.matrix(cbind(1,Xmat))

X2    = as.matrix(cbind(1,Xmat2))

Y     = data[resp]

fit   = X%*%second.lm\$coefficients

res   = Y - fit

xPx     = t(X2)%*%X2

xPx.inv = solve(xPx)

z\$cov.unscaled     = xPx.inv

z\$residuals        = res

z\$sigma            = sqrt(mean(res^2))

varcovmat          = z\$cov.unscaled*z\$sigma

coef               = z\$coefficients[,1]

IV.SE              = z\$sigma*sqrt(diag(xPx.inv))

t.iv               = coef/IV.SE

p.val              = 2*(1-pnorm(abs(t.iv)))

z\$coefficients[,2] = IV.SE

z\$coefficients[,3] = t.iv

z\$coefficients[,4] = p.val

result = list(summary(first.lm),z,ftest)

names(result) = c("first", "second", "ftest")

print("IV object successfully created. Use sum.iv() on object")

return(invisible(result))

}

sum.iv = function(reg.iv, first=FALSE,

second=TRUE, ftest=FALSE) {

x= rep(0,2)

if(first==TRUE) x[1] = 1

if(second==TRUE) x[2]= 2

if(ftest==TRUE) x[3]= 3

print(reg.iv[x])

}