An Application of boot() to IV regression

Bootstrapping standard errors can be a useful technique when obtaining closed form for the standard error formula is difficult or intractable. In this post, I give an example of how to use R to create a bootstrap sampling distribution in the context of IV regression. Specifically, I use boot() to automatically augment a function of mine to resample the indices of my data set with replacement (see the code in the function below).

In my application, I present a function that uses the boot library to report bootstrap standard errors for an instrumental variables regression. My ivboot() function builds on the iv() command I wrote for the tonymisc library.*

Now, onto the ivboot() example. Start by loading the tonymisc library and some data from the library, as well as the boot library.



We can run the following R script to define the function ivboot().



When applied to data and an IV regression model, the ivboot() function creates an IV regression object that -- when tonymisc is loaded -- is compatible with mtable() output. The only difference between an ivboot()-created object and an iv() object is that the ivboot() object has standard errors that are based on the bootstrap distribution of the coefficient estimates (Statistics other than the second-stage coefficient estimates are not bootstrapped).

Here is some code to compare the bootstrap output to the analytical standard error output:



On this small sample (N=100) of simulated data, the ivboot() command took less than a minute to run on my computer (timing of may vary depending on your computer). For much larger data sets, this will be slower. If you have a larger problem or lower standards (or higher standards, your choice), you can use the boots option to ivboot() to specify the number of bootstrap samples. Currently, I have set the default to 500, but you could specify boots = 200 if you want the command to run faster (boots = 10 will make it run even faster, but I don't recommend that!).

Here is the mtable() output, which can easily be ported into LaTeX using the toLatex() command.


*This standard output from an mtable() extension to my iv() command provides quite a bit of information in a convenient format. Another nice feature of iv() is that iv()-created objects have first-stage summary information readily stored in the object for extraction and analysis.

Video Tutorial on IV Regression

Update: [1 May 2011] I am working on a better augmentation of the current IV regression functions (specifically ivreg() in AER) in R. I will post a link here to my new method/function for IV regression when I finish debugging the code.

Update 2: [15 May 2011] The function I wrote here does not work with factors due to the way in which I constructed fitted values for the standard error correction in the second stage.

I recorded a new video tutorial on how to conduct an IV regression. This time I used my new ivregress() command, which has much better syntax than my old one. There are several reasons I needed to do this:

1. Relative to my previous post on IV regression, I have added the ability to conduct overidentifying tests. Now, my ivregress() command has all of the basic functionality of Stata's 2sls option to ivregress.
2. Relative to my previous video tutorial on IV regression, this video uses my newer command, which as much better syntax. As such, I will use this video tutorial to replace the old one.
   
   

Here is the R script file I used to define my ivregress() command and its companion summary command sum.iv().



To get ivregress() and sum.iv() to work, (1) copy the text in the above text box into an R script, (2) run all (or just download and run the R script) and (3) save into your default workspace if you want to use the commands into the future. Then, you're good to go. See my video demonstrating how to use the command here:



Here is the code I used to produce the output in the video.



Finally, here is a link to the synthetic "market share" data I used in the video.

I hope you find this command helpful.

An easier to use IV regression command in R

Update: I have added some functionality to my ivregress() command. Check out my newer post here.

After I posted my last video tutorial on how to use my IV regression function, I received a comment asking why I didn't write the command a different way to make the syntax easier to read.

The answer is that I didn't know how to write an easier to use function a year ago (when I wrote the ivreg() function). After some digging, I figured out how to work with "formula objects" in R and the result is an easier to use IV regression function (called ivregress()).

How to "install" ivregress()

Here's the code you need to run to define ivregress() and its companion summary command sum.iv(). This will provide instrumental variables regression estimates if you have one endogenous regressor with one or more instruments.



You only need to run the above code once to define the function object ivregress() for all future uses (as long as you don't write over it or clear your workspace).

How to use ivregress()

To create an "IV object" simply use the ivregress command as follows (the command relies on the car library):

library(car)
myivobject = ivregress(Y~X1+X2+...+Xk, Xj ~ Z1+Z2+...+Zl, dataframe)

where X1, X2,...Xk are second stage regressors, Xj is an endogenous regressor for which you would like to instrument, and Z1, Z2,...,Zl are instruments.

Then use sum.iv() from the previous posts to produce summary output. Sometime soon, I will post a video tutorial showing how to use ivregress() to perform IV regression.

Video Tutorial on Instrumental Variables in R

Update: I have replaced this video tutorial with a video tutorial on a newer, easier to use IV regression command. Check out that command here.

In this video, I show how to use my instrumental variables function in R, ivreg(), along with its companion summary command, sum.iv().



Here is a link to the synthetic "market share" data I used in the video.

Two-stage least squares

Update: Since this post, I discovered a bunch of better ways to conduct 2SLS regression in R. These methods do not allow the user to look under the hood to see what is going on (so this post may be instructive anyway). If you're interested, ivreg() in AER, tsls() in sem and systemfit() in systemfit are good options. My package tonymisc also has a command, iv() that works when some or all of these fail. I'm working on making that function better (5/18/2011)

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.