Companion for OLS Chapter

This companion to the OLS Chapter shows how to implement the models in both Stata and R. Unless stated otherwise, the code shown here is for R. Note, we use a panel dataset (multiple time periods per individual) and arbitrarily restrict our analysis to a cross section dataset by analyzing only records where time is 4. We do this for demonstration purposes only- it is usually a horrible idea to throw away data like this, and later in the course we will learn more about panel data techniques.

# allows us to open stata datasets
library(foreign)
# for heteroskedasticity corrected variance/covariance matrices
library(sandwich)
# for White Heteroskedasticity Test
library(lmtest)
# for bootstrapping
library(boot)
# for column percentiles
library(matrixStats)
tk.df = read.dta("http://rlhick.people.wm.edu/econ407/data/tobias_koop.dta")
tk4.df = subset(tk.df, time == 4)
histwage <- hist(tk4.df$ln_wage,main = "Log Wage",xlab="ln_wage")

histwage.png

The OLS Model

Using ln_wage as our dependent variable, here is the code for running a simple OLS model in stata and R. First in Stata. Since we haven't loaded the data yet, let's do that and run the model:

webuse set "http://rlhick.people.wm.edu/econ407/data/"
(prefix now "http://rlhick.people.wm.edu/econ407/data")
webuse tobias_koop
keep if time==4
(16885 observations deleted)
regress ln_wage educ pexp pexp2 broken_home

      Source |       SS       df       MS              Number of obs =    1034
-------------+------------------------------           F(  4,  1029) =   51.36
       Model |  37.3778146     4  9.34445366           Prob > F      =  0.0000
    Residual |   187.21445  1029  .181938241           R-squared     =  0.1664
-------------+------------------------------           Adj R-squared =  0.1632
       Total |  224.592265  1033  .217417488           Root MSE      =  .42654

------------------------------------------------------------------------------
     ln_wage |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        educ |   .0852725   .0092897     9.18   0.000     .0670437    .1035014
        pexp |   .2035214   .0235859     8.63   0.000     .1572395    .2498033
       pexp2 |  -.0124126   .0022825    -5.44   0.000    -.0168916   -.0079336
 broken_home |  -.0087254   .0357107    -0.24   0.807    -.0787995    .0613488
       _cons |   .4603326    .137294     3.35   0.001     .1909243    .7297408
------------------------------------------------------------------------------

In R, we use

ols.lm = lm(ln_wage~educ + pexp + pexp2 + broken_home, data=tk4.df)
summary(ols.lm)
Call:
lm(formula = ln_wage ~ educ + pexp + pexp2 + broken_home, data = tk4.df)

Residuals:
     Min       1Q   Median       3Q      Max
-1.69490 -0.24818 -0.00079  0.25417  1.51139

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.460333   0.137294   3.353 0.000829 ***
educ         0.085273   0.009290   9.179  < 2e-16 ***
pexp         0.203521   0.023586   8.629  < 2e-16 ***
pexp2       -0.012413   0.002283  -5.438 6.73e-08 ***
broken_home -0.008725   0.035711  -0.244 0.807020
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4265 on 1029 degrees of freedom
Multiple R-squared:  0.1664,	Adjusted R-squared:  0.1632
F-statistic: 51.36 on 4 and 1029 DF,  p-value: < 2.2e-16

Robust Standard Errors

Note that since the \(\mathbf{b}\)'s are unbiased estimators for \(\beta\), those continue to be the ones to report. However, our variance covariance matrix is $$ VARCOV_{HC1} = \mathbf{(x'x)^{-1}V(x'x)^{-1}} $$

Defining \(\mathbf{e=y-xb}\), $$ \mathbf{V} = diag(\mathbf{ee}') $$

In stata, the robust option is a very easy way to implement the robust standard error correction. Note: this is not robust regression.

regress ln_wage educ pexp pexp2 broken_home, robust

Linear regression                                      Number of obs =    1034
                                                       F(  4,  1029) =   64.82
                                                       Prob > F      =  0.0000
                                                       R-squared     =  0.1664
                                                       Root MSE      =  .42654

------------------------------------------------------------------------------
             |               Robust
     ln_wage |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        educ |   .0852725   .0099294     8.59   0.000     .0657882    .1047568
        pexp |   .2035214   .0226835     8.97   0.000     .1590101    .2480326
       pexp2 |  -.0124126   .0022447    -5.53   0.000    -.0168173   -.0080079
 broken_home |  -.0087254   .0312381    -0.28   0.780     -.070023    .0525723
       _cons |   .4603326   .1315416     3.50   0.000     .2022121     .718453
------------------------------------------------------------------------------

Implementing this in R, we can use the sandwich package's `vcovHC` function with the `HC1` option (stata default):

VARCOV = vcovHC(ols.lm, 'HC1')
se = sqrt(diag(VARCOV))
print(se)
(Intercept)        educ        pexp       pexp2 broken_home
0.131541623 0.009929445 0.022683530 0.002244678 0.031238081

Note that the F-statistic reported by the R OLS result above is inappropriate in the presence of heteroskedasticity. We need to use a different form of the F test. Stata does this automatically. In R, we need to use waldtest for the correct model F statistic:

waldtest(ols.lm, vcov = vcovHC(ols.lm,'HC1'))
Wald test

Model 1: ln_wage ~ educ + pexp + pexp2 + broken_home
Model 2: ln_wage ~ 1
  Res.Df Df     F    Pr(>F)
1   1029
2   1033 -4 64.82 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

So far, we have gotten R to provide us with the correct robust standard error regression results but it isn't in a single table. To put it together for printing, use

# copy summary object
robust.summ = summary(ols.lm)
# overwrite all coefficients with robust counterparts
robust.summ$coefficients = unclass(coeftest(ols.lm, vcov. = VARCOV))
# overwrite f stat with correct robust value
f_robust = waldtest(ols.lm, vcov = VARCOV)
robust.summ$fstatistic[1] = f_robust$F[2]
print(robust.summ)
Call:
lm(formula = ln_wage ~ educ + pexp + pexp2 + broken_home, data = tk4.df)

Residuals:
     Min       1Q   Median       3Q      Max
-1.69490 -0.24818 -0.00079  0.25417  1.51139

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.460333   0.131542   3.500 0.000486 ***
educ         0.085273   0.009929   8.588  < 2e-16 ***
pexp         0.203521   0.022684   8.972  < 2e-16 ***
pexp2       -0.012413   0.002245  -5.530 4.06e-08 ***
broken_home -0.008725   0.031238  -0.279 0.780057
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4265 on 1029 degrees of freedom
Multiple R-squared:  0.1664,	Adjusted R-squared:  0.1632
F-statistic: 64.82 on 4 and 1029 DF,  p-value: < 2.2e-16

The bootstrapping approach to robust standard errors

Given the advent of very fast computation even with laptop computers, we can use boostrapping to correct our standard errrors for any non-iid error process. In R, we run this (and by asking for percentiles) to get non-parametric bootstrap estimates for confidence intervals:

bs_ols <- function(formula, data, indices) {
  d <- data[indices,] # allows boot to select sample
  fit <- lm(formula, data=d) # runs the model for this sample
  return(coef(fit)) # returns coefficients
}

results <- boot(data=tk4.df, statistic=bs_ols,
	R=1000, formula=ln_wage~educ + pexp + pexp2 + broken_home)

From that we can extract the 95\% confidence intervals from the 1000 replicates.

ci_mat = colQuantiles(results$t,probs=c(.025,.975))
ci = data.frame(ci_mat,row.names=c("Intercept","educ","pexp","pexp2","broken_home"))
colnames(ci) = c("2.5%","97.5%")
print(ci)
                   2.5%        97.5%
Intercept    0.19939776  0.714556678
educ         0.06565962  0.104340132
pexp         0.16128174  0.249996019
pexp2       -0.01696095 -0.008148125
broken_home -0.07249099  0.056744122

In Stata, we just run (for 1000 replicates) using the following code. I believe this is doing what is called a parametric bootstrap (confidence intervals assume normality). As you can see, the Stata method is much easier to implement.

regress ln_wage educ pexp pexp2 broken_home, vce(bootstrap,rep(1000))
   50
..................................................   100
..................................................   150
..................................................   200
..................................................   250
..................................................   300
..................................................   350
..................................................   400
..................................................   450
..................................................   500
..................................................   550
..................................................   600
..................................................   650
..................................................   700
..................................................   750
..................................................   800
..................................................   850
..................................................   900
..................................................   950
..................................................  1000

Linear regression                               Number of obs      =      1034
                                                Replications       =      1000
                                                Wald chi2(4)       =    267.45
                                                Prob > chi2        =    0.0000
                                                R-squared          =    0.1664
                                                Adj R-squared      =    0.1632
                                                Root MSE           =    0.4265

------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
     ln_wage |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        educ |   .0852725   .0099813     8.54   0.000     .0657095    .1048355
        pexp |   .2035214   .0224216     9.08   0.000     .1595759    .2474669
       pexp2 |  -.0124126   .0022395    -5.54   0.000     -.016802   -.0080232
 broken_home |  -.0087254   .0323746    -0.27   0.788    -.0721785    .0547277
       _cons |   .4603326   .1313786     3.50   0.000     .2028352    .7178299
------------------------------------------------------------------------------

Testing for Heteroskedasticity

The two dominant heteroskedasticity tests are the White test and the Cook-Weisberg (both Breusch Pagan tests). The White test is particularly useful if you believe the error variances have a potentially non-linear relationship with model covariates. For stata, we need to make sure the basic OLS regression (rather than robust) is run just prior to these commands and then we can use:

quietly regress ln_wage educ pexp pexp2 broken_home
imtest, white

White's test for Ho: homoskedasticity
         against Ha: unrestricted heteroskedasticity

         chi2(12)     =     29.20
         Prob > chi2  =    0.0037

Cameron & Trivedi's decomposition of IM-test

---------------------------------------------------
              Source |       chi2     df      p
---------------------+-----------------------------
  Heteroskedasticity |      29.20     12    0.0037
            Skewness |      14.14      4    0.0068
            Kurtosis |      16.19      1    0.0001
---------------------+-----------------------------
               Total |      59.53     17    0.0000
---------------------------------------------------

In R, we perform the White test by:

bptest(ols.lm,ln_wage ~ educ*pexp + educ*pexp2 + educ*broken_home + pexp*pexp2 + pexp*broken_home +
	 pexp2*broken_home + I(educ^2) + I(pexp^2) + I(pexp2^2) + I(broken_home^2),data=tk4.df)
        studentized Breusch-Pagan test

data:  ols.lm
BP = 29.196, df = 12, p-value = 0.003685

Notice that in R, we need to specify the form of the heteroskedasticity. To mirror the White test, we include interaction and second order terms. Both Stata's imtest and R's bptest allow you to specify more esoteric forms of heteroskedasticity, although in R it is more apparent since we have to specify the form of the heteroskedasticity for the test to run.

A more basic test- one assuming that heteroskedastity is a linear function of the \(\mathbf{x}\)'s- is the Cook-Weisberg Breusch Pagan test.

In stata, we use

hettest

Breusch-Pagan / Cook-Weisberg test for heteroskedasticity
         Ho: Constant variance
         Variables: fitted values of ln_wage

         chi2(1)      =    19.57
         Prob > chi2  =   0.0000

R doesn't implement this test directly. But it can be performed manually following these steps:

  • Recover residuals as above, denote as \(\mathbf{e}\)
  • square the residuals by the sqare root of mean squared residuals:

$$ \mathbf{r} = \frac{diagonal(\mathbf{ee}')}{\sqrt{\mathbf{e'e}/N}} $$

  • Run the regression \(\mathbf{r} = \gamma_0 + \gamma_1 \mathbf{xb} + \mathbf{v}\). Recall that \(\mathbf{xb} = \mathbf{\hat{y}}\).
  • Recover ESS defined as

$$ ESS = (\mathbf{\hat{r}} - \bar{r})'(\hat{\mathbf{r}} - \bar{r}) $$

  • ESS/2 is the test statistic distributed \(\chi^2_{(1)}\)
r = diag(tcrossprod(ols.lm$residuals))/mean(diag(tcrossprod(ols.lm$residuals)))
xb = ols.lm$fitted.values
het_ols_results = lm(r~xb)
mss = crossprod(het_ols_results$fitted.values - mean(r))
cw = mss/2
sprintf("BP / Cook-Weisberg score for heterosk.: %2.3f distributed chi2 with %d df.",cw,1)
[1] "BP / Cook-Weisberg score for heterosk.: 19.572 distributed chi2 with 1 df."