Here is some code I wrote that allows you to easily perform two-stage least squares regression using R.

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

}

The above function is ivreg(), which takes five arguments: (1) Y0 is the name of the dependent variable, (2) X0 is the name of the endogenous regressor, (3) Z10 is the set of names of exogenous regressors [second stage variables assumed to be exogenous], and (4) Z20 is the set of instruments.

The function will not print any output if run alone. The appropriate syntax is to save the function output into an "iv" object, which you can summarize using the following function:

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

}

Here is an example of how to use the functions to perform a two-stage least squares regression. Imagine that you have a data on prices and quantities of cars sold, and you have some data on attributes that directly affect quantity sold (falls into category Z10) and attributes that only affect quantity sold through the price (falls into category Z20).

You could estimate the model, and summarize its output using the following two commands:

demand.iv = ivreg("quantity", "price", c("horsepower", "mpg"), c("regulations", "cost", "competition"))

sum.iv(demand.iv, first=FALSE, second=TRUE, ftest=FALSE)

Note: The quotes in the ivreg command are necessary to get the function to work properly. This is an unusual artifact of how I wrote the function.

What is nice is that the sum.iv() command allows the user to specify what output he/she wants to see from the two-stage least squares regression. For the sum.iv() command, you could just leave off the logical arguments to get the same output (because these are the defaults), but you can choose to show the first-stage estimates and the F-test on the explanatory power of the instruments by specifying your own options.