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

  print("to learn about your 2SLS Regression")

  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])

}