# Companion for OLS Chapter

This page has been moved to https://econ.pages.code.wm.edu/407/notes/docs/index.html and is no longer being maintained here.

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("https://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")  ## 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 "https://rlhick.people.wm.edu/econ407/data/" webuse tobias_koop keep if time==4 regress ln_wage educ pexp pexp2 broken_home  set more off (prefix now "https://rlhick.people.wm.edu/econ407/data") (16,885 observations deleted) Source | SS df MS Number of obs = 1,034 -------------+---------------------------------- F(4, 1029) = 51.36 Model | 37.3778146 4 9.34445366 Prob > F = 0.0000 Residual | 187.21445 1,029 .181938241 R-squared = 0.1664 -------------+---------------------------------- Adj R-squared = 0.1632 Total | 224.592265 1,033 .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 = 1,034 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))  (running regress on estimation sample) Bootstrap replications (1000) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .................................................. 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 = 1,034 Replications = 1,000 Wald chi2(4) = 252.66 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 .0100503 8.48 0.000 .0655743 .1049707 pexp | .2035214 .022688 8.97 0.000 .1590538 .247989 pexp2 | -.0124126 .0022409 -5.54 0.000 -.0168047 -.0080205 broken_home | -.0087254 .0320888 -0.27 0.786 -.0716183 .0541676 _cons | .4603326 .1340845 3.43 0.001 .1975318 .7231333 ------------------------------------------------------------------------------  ## 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."