ivreg = function(Y0,X0,Z10=NULL,Z20,data.df){

Y  = data.df[[Y0]]

X  = data.df[[X0]]

Z1 = NULL

if(length(Z10)>0){

Z1 = data.df[[Z10[1]]]

}

Z2 = data.df[[Z20[1]]]

if(length(Z10)>1){

for(i in 2:length(Z10)){

Z1=cbind(Z1,data.df[[Z10[i]]])

}

}

if(length(Z20)>1){

for(i in 2:length(Z20)){

Z2 =cbind(Z2,data.df[[Z20[i]]])

}

}

## First Stage

Z       = cbind(1,Z1,Z2)

zPz     = t(Z)%*%Z

zPx     = t(Z)%*%X

zPz.inv = solve(zPz)

beta.f  = zPz.inv%*%zPx

xhat    = Z%*%beta.f

SSE.f = sum((X - xhat)^2)

df.d    = length(data.df[,1])-length(Z[1,])

MSE.f   = SSE.f/df.d

RMSE.f  = sqrt(MSE.f)

se.f    = RMSE.f*sqrt(diag(zPz.inv))

t.f     = beta.f/se.f

pval.f  = 2*pnorm(abs(t.f),lower.tail=F)

## Constructing F-test

Z.d       = cbind(rep(1,

length(data.df[,1])),Z1)

zdPzd     = t(Z.d)%*%Z.d

zdPx      = t(Z.d)%*%X

zdPzd.inv = solve(zdPzd)

bhat.d    = zdPzd.inv%*%zdPx

xhat.d    = Z.d%*%bhat.d

SSE.d      = sum((X-xhat.d)^2)

SSEx.z     = SSE.d -SSE.f

MSEx       = SSEx.z/length(Z20)

F.stat     = MSEx/MSE.f

df.n       = length(Z[1,]) -

length(Z.d[1,])

pval.Ftest = pf(F.stat, df.n, df.d ,

lower.tail=FALSE)

ftest  = round(c(F.stat,

pval.Ftest,df.n), 4)

names(ftest) = c("F-stat", "P-value",

"# of instruments")

rn.f = c("constant",Z10,Z20)

first   = data.frame(cbind(beta.f,se.f,

t.f, pval.f),row.names=rn.f)

names(first) = c("Estimate", "Std. Error",

"T-stat", "P-value")

## Second Stage

X.mat   = cbind(1,Z1,xhat)

X.mat2  = cbind(1,Z1,X)

xPx      = t(X.mat)%*%X.mat

xPx.inv  = solve(xPx)

xPy      = t(X.mat)%*%Y

beta2sls = xPx.inv%*%xPy

resid    = Y-X.mat2%*%beta2sls

MSE      = mean(resid^2)

RMSE     = sqrt(MSE)

se2sls   = RMSE*sqrt(diag(xPx.inv))

tstat    = beta2sls/se2sls

pval     = 2*pnorm(abs(tstat),lower.tail=F)

rn = c("constant",Z10,X0)

second   = data.frame(cbind(beta2sls,se2sls,

tstat, pval),row.names=rn)

names(second) = c("Estimate", "Std. Error",

"T-stat", "P-value")

result = list(first,second,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])

}