Question 1: Panel Data Models: Fixed and Random Effects

In the STAR experiment, children were randomly assigned within schools into three types of classes: small classes with 13–17 students, regular-sized classes with 22–25 students, and regular-sized classes with a full-time teacher aide to assist the teacher. Student scores on achievement tests were recorded, as was some information about the students, teachers, and schools. Data for the kindergarten classes are contained in the data file star.dat.

Let’s load and view the data’s first few lines as well as summary statistics about it

# load the data
load("./data/star.rda")
kable(head(star))
id schid tchid tchexper absent readscore mathscore totalscore boy white_asian black tchwhite tchmasters freelunch schurban schrural small regular aide
10133 169280 16928003 7 5 427 478 905 1 0 1 1 0 0 0 0 0 0 1
10246 218562 21856202 8 28 450 494 944 0 1 0 1 1 0 1 0 0 0 1
10263 205492 20549204 3 2 483 513 996 0 0 1 1 0 1 0 0 1 0 0
10266 257899 25789904 12 10 456 513 969 1 1 0 1 1 0 0 1 0 1 0
10275 161176 16117602 2 3 411 468 879 1 1 0 1 0 0 0 1 0 0 1
10281 189382 18938204 7 2 443 473 916 1 1 0 1 1 0 0 0 1 0 0 
kable(tidy(star))
## Warning: 'tidy.data.frame' is deprecated.
## See help("Deprecated")
column n mean sd median trimmed mad min max range skew kurtosis se
id 5786 1.559306e+04 2.694317e+03 15380.5 1.543933e+04 1862 10133 21580 11447 0.4033242 2.498993 3.542089e+01
schid 5786 2.110018e+05 3.838193e+04 215533.0 2.138217e+05 29243 112038 264945 152907 -0.5868482 2.365732 5.045887e+02
tchid 5786 2.110019e+07 3.838193e+06 21553306.0 2.138217e+07 2924296 11203801 26494505 15290704 -0.5868482 2.365732 5.045887e+04
tchexper 5766 9.306452e+00 5.767684e+00 9.0 9.003034e+00 4 0 27 27 NA NA 7.595640e-02
absent 5765 1.027511e+01 9.270640e+00 8.0 8.893778e+00 5 0 79 79 NA NA 1.220984e-01
readscore 5786 4.367297e+02 3.171347e+01 433.0 4.338069e+02 19 315 627 312 1.3400597 6.844684 4.169217e-01
mathscore 5786 4.855990e+02 4.769394e+01 484.0 4.831495e+02 30 320 626 306 0.4760213 3.294667 6.270092e-01
totalscore 5786 9.223287e+02 7.374660e+01 915.0 9.179244e+02 46 635 1253 618 0.6301511 3.745125 9.695111e-01
boy 5786 5.134808e-01 4.998614e-01 1.0 5.168467e-01 0 0 1 1 -0.0539429 1.002910 6.571400e-03
white_asian 5786 6.766333e-01 4.678018e-01 1.0 7.207343e-01 0 0 1 1 -0.7552281 1.570369 6.150000e-03
black 5786 3.209471e-01 4.668809e-01 0.0 2.762419e-01 0 0 1 1 0.7670838 1.588418 6.137900e-03
tchwhite 5786 8.354649e-01 3.707925e-01 1.0 9.192225e-01 0 0 1 1 -1.8096048 4.274669 4.874600e-03
tchmasters 5786 3.517110e-01 4.775456e-01 0.0 3.146868e-01 0 0 1 1 0.6211000 1.385765 6.278100e-03
freelunch 5786 4.816799e-01 4.997074e-01 0.0 4.771058e-01 0 0 1 1 0.0733296 1.005377 6.569400e-03
schurban 5786 3.128241e-01 4.636834e-01 0.0 2.660907e-01 0 0 1 1 0.8074134 1.651917 6.095800e-03
schrural 5786 4.709644e-01 4.991994e-01 0.0 4.637149e-01 0 0 1 1 0.1163387 1.013535 6.562700e-03
small 5786 3.003802e-01 4.584629e-01 0.0 2.505400e-01 0 0 1 1 0.8708971 1.758462 6.027200e-03
regular 5786 3.465261e-01 4.759043e-01 0.0 3.082073e-01 0 0 1 1 0.6450337 1.416068 6.256500e-03
aide 5786 3.530937e-01 4.779728e-01 0.0 3.164147e-01 0 0 1 1 0.6147589 1.377928 6.283700e-03

1.b) Reestimate the model in part (1.a) with school fixed effects. Compare the results with those in part (1.a). Have any of your conclusions changed?**

# added school fixed effects by "factor(school_id)"
math.fixed <- lm(mathscore ~ small + aide + tchexper + 
                 boy + white_asian + factor(schid)-1, 
                  data = star)

School Fixed effects Model: Classroom Features effects on Math Scores

term estimate std.error statistic p.value
small 9.3495834 1.3970221 6.692509 0.0000000
aide 0.5268855 1.3491170 0.390541 0.6961512
tchexper 0.4201512 0.1083858 3.876441 0.0001072
boy -6.6312110 1.1133604 -5.956033 0.0000000
white_asian 23.6509406 2.3108789 10.234608 0.0000000
factor(schid)112038 433.4416543 6.1915004 70.005915 0.0000000
factor(schid)123056 448.2977456 6.4407403 69.603450 0.0000000
factor(schid)128068 443.4064908 6.2685431 70.735175 0.0000000
factor(schid)128076 443.4736630 6.0680868 73.082947 0.0000000
factor(schid)128079 434.9700466 6.0159497 72.302807 0.0000000
factor(schid)130085 455.0002352 5.7344680 79.344803 0.0000000
factor(schid)159171 485.2103480 4.8297836 100.462131 0.0000000
factor(schid)161176 450.5546086 5.4446595 82.751659 0.0000000
factor(schid)161183 476.6336042 4.8551802 98.170116 0.0000000
factor(schid)162184 460.0258998 5.9868061 76.839952 0.0000000
factor(schid)164198 464.5300240 6.2850936 73.909802 0.0000000
factor(schid)165199 490.5138727 6.4041624 76.592978 0.0000000
factor(schid)166203 446.1648065 5.6974915 78.308991 0.0000000
factor(schid)168211 459.1896239 5.2751913 87.047008 0.0000000
factor(schid)168214 485.3243089 6.2627177 77.494202 0.0000000
factor(schid)169219 484.7343721 6.7015248 72.331953 0.0000000
factor(schid)169229 473.5326526 4.4214784 107.098263 0.0000000
factor(schid)169231 438.6436570 6.0257813 72.794487 0.0000000
factor(schid)169280 462.4427156 6.2918758 73.498386 0.0000000
factor(schid)170295 461.7120027 5.9148825 78.059371 0.0000000
factor(schid)173312 491.5293455 6.1513337 79.906143 0.0000000
factor(schid)176329 479.6244542 5.8925086 81.395631 0.0000000
factor(schid)180344 465.2790091 4.8702851 95.534246 0.0000000
factor(schid)189378 439.0862883 5.7169389 76.804440 0.0000000
factor(schid)189382 450.0979675 6.0424637 74.489147 0.0000000
factor(schid)189396 444.6735027 5.9280641 75.011588 0.0000000
factor(schid)191411 453.1603965 6.8780427 65.885081 0.0000000
factor(schid)193422 473.8613131 6.0915237 77.790277 0.0000000
factor(schid)193423 474.4369684 5.6083034 84.595454 0.0000000
factor(schid)201449 480.0881863 4.6273530 103.750067 0.0000000
factor(schid)203452 446.9405868 5.2891302 84.501718 0.0000000
factor(schid)203457 469.3473431 7.0753889 66.335201 0.0000000
factor(schid)205488 456.0238430 6.2316347 73.178848 0.0000000
factor(schid)205489 462.0905394 6.3582369 72.675893 0.0000000
factor(schid)205490 425.2598669 5.6116565 75.781521 0.0000000
factor(schid)205491 460.2467039 5.7835045 79.579207 0.0000000
factor(schid)205492 494.3397358 5.0780955 97.347467 0.0000000
factor(schid)208501 462.6245100 6.0193562 76.856144 0.0000000
factor(schid)208503 443.5864645 6.0212478 73.670189 0.0000000
factor(schid)209510 441.2641751 5.1770397 85.234844 0.0000000
factor(schid)212522 470.0526927 5.9385167 79.153216 0.0000000
factor(schid)215533 489.4558886 4.5342554 107.946255 0.0000000
factor(schid)216536 462.9346546 4.9934517 92.708347 0.0000000
factor(schid)218562 475.0435181 6.1369286 77.407372 0.0000000
factor(schid)221571 417.6265317 4.7687794 87.575141 0.0000000
factor(schid)221574 446.7331245 5.8582226 76.257452 0.0000000
factor(schid)225585 450.9408926 5.5539973 81.192134 0.0000000
factor(schid)228606 479.7092685 5.4446322 88.106827 0.0000000
factor(schid)230612 496.0533470 6.0411844 82.111936 0.0000000
factor(schid)231616 476.0860030 6.3219494 75.306836 0.0000000
factor(schid)234628 464.1989788 5.0906699 91.186226 0.0000000
factor(schid)244697 441.0988498 4.5819219 96.269395 0.0000000
factor(schid)244708 443.2820927 4.1997491 105.549661 0.0000000
factor(schid)244723 444.5766881 4.3098299 103.154114 0.0000000
factor(schid)244727 474.5772920 4.9480977 95.911058 0.0000000
factor(schid)244728 446.4440489 6.1128748 73.033403 0.0000000
factor(schid)244736 501.1558006 6.0305482 83.102860 0.0000000
factor(schid)244745 492.9133982 5.2290374 94.264654 0.0000000
factor(schid)244746 478.4434703 5.9198110 80.820734 0.0000000
factor(schid)244755 484.7687833 3.9533237 122.623095 0.0000000
factor(schid)244764 471.0968437 7.3331620 64.241979 0.0000000
factor(schid)244774 474.2716313 4.3675713 108.589327 0.0000000
factor(schid)244776 461.6985601 3.8661461 119.420876 0.0000000
factor(schid)244780 545.4808369 5.7576685 94.739882 0.0000000
factor(schid)244796 458.3082722 5.8530133 78.302961 0.0000000
factor(schid)244799 468.4503308 5.4401940 86.109123 0.0000000
factor(schid)244801 473.5923761 5.3249229 88.938822 0.0000000
factor(schid)244806 500.9792528 3.8218885 131.081597 0.0000000
factor(schid)244818 455.1998711 4.9339839 92.258079 0.0000000
factor(schid)244831 462.0976048 5.9883646 77.165910 0.0000000
factor(schid)244839 519.3557376 5.1045488 101.743711 0.0000000
factor(schid)252885 469.8948803 5.6770513 82.770941 0.0000000
factor(schid)253888 466.4495285 7.2087178 64.706310 0.0000000
factor(schid)257899 447.7188038 5.0124125 89.322019 0.0000000
factor(schid)257905 458.8989956 4.6675227 98.317464 0.0000000
factor(schid)259915 461.0594547 6.1924181 74.455478 0.0000000
factor(schid)261927 477.8229310 5.3460476 89.378727 0.0000000
factor(schid)262937 496.4755341 5.9311651 83.706241 0.0000000
factor(schid)264945 481.7808675 5.0198832 95.974517 0.0000000

Technical comments: The function factor() generates dummy variables for all categories of the variable, taking the first category as the reference. To include the reference in the output, one needs to exclude the constant from the regression model by including the term −1 in the regression formula. When the constant is not excluded, the coefficients of the dummy variables represent, as usual, the difference between the respective category and the benchmark one.

However in R we can use the “plm” (panel) models to build a model object inherent to that idea:

# fixed effects using plm "within" method
math.fixed.within <- plm(mathscore ~ small + aide + 
                           tchexper + boy + white_asian,
                          data = star,
                         model="within",
                         index=c("schid"))

kable(tidy(math.fixed.within))
term estimate std.error statistic p.value
small 9.3495834 1.3970221 6.692509 0.0000000
aide 0.5268855 1.3491170 0.390541 0.6961512
tchexper 0.4201512 0.1083858 3.876441 0.0001072
boy -6.6312110 1.1133604 -5.956033 0.0000000
white_asian 23.6509406 2.3108789 10.234608 0.0000000

This model showed still no significance of AIDE, but the school fixed effects were very (the most possible) significant.
Including the fixed effects (accounting for individual heterogeneity) significantly lowers the marginal effects of the variables. That is to be expected, since the previously the variables needed to explain all (or as much as possible) of the variance between individuals; but now belonging to each particular school can explain much of that variance.

1.c) Test for the significance of the school fixed effects. Under what conditions would we expect the inclusion of significant fixed effects to have little influence on the coefficient estimates of the remaining variables?

I will test it with the F test for fixed effects:

# Fixed effects test: Ho='No fixed effects'
kable(tidy(pFtest(math.fixed.within, math.pooled)), caption="Fixed effects test: Ho='No fixed effects'")
Fixed effects test: Ho=‘No fixed effects’
df1 df2 statistic p.value method alternative
78 5682 18.06606 0 F test for individual effects significant effects

The fixed effects are all significant. By using the F Test for Individual Effects, with the Null Hypothesis that the fixed effects are zero, the F test shows that the p-value is 0 and the null hypothesis can be rejected, meaning that the fixed effects are indeed significant.
As further verification, we could already detecte that the regression output results for the fixed effects regression model showed that they all have the lowest possible p-value (under 2e-16), which is lower than the other variables in the models.

This means that studying in each particular school accounts for more variance than each of the other variables does, and can better explain the variation in Math outcomes. This stems from the fact that each school is different on many other attributes, including hidden variables that we don’t have data on (such as potentially the demographic composition of the school, area in town, income level distribution of the parents, or even number of hours taught math).
We would expect fixed effects to have little influence on the coefficient estimates of the remaining variables in cases where they don’t help to explain the variance better than the other variables included (although we shouldn’t include other variables which are constant per entity when including fixed effects).

1.d) Reestimate the model in part (a) with school random effects. Compare the results with those from parts (a) and (b). Are there any variables in the equation that might be correlated with the school effects?

The random effects model elaborates on the fixed effects model by recognizing that, since the individuals in the panel are randomly selected, their characteristics, measured by the intercept β1i should also be random. Thus, the random effects model assumes the form of the intercept is this:
\(β_{1i}=\barβ_1+u_i\)
where \(\barβ_1\) stands for the population average and \(u_i\) represents an individual-specific random term. As in the case of fixed effects, random effects are also time-invariant.
So the random effects model is: \(y_{it}=\barβ_1 +β_2 X_{2it}+ν_{it}\)
Where the error termm, \(ν_it\) , incorporates both individual specifics and the initial regression error term:
\(ν_{it}=u_i+e_{it}\)
Hence, the intercept is constant or averaged as opposed to the individual intercepts of the fixed effects model, and the inidividual effect \(e_{it}\) which was part of the intercept at the fixed effects model is now “moved” to the error term and considered part of individual randomness/errors. Testing for random effects amounts to testing the hypothesis that there are no differences among individuals, which implies that the individual-specific random variable has zero variance.

First, let’s test wether a random effects model is appropriate:

# The same function we used for fixed effects can be used for random effects, but setting the argument model= to ‘random’ and selecting a random.method.
mathReTest <- plmtest(math.pooled, effect="individual", type="honda")
# "Random Effects Test for the Math Schore Model"
kable(tidy(mathReTest), caption="Random Effects Lagrange Multiplier Test for the Math Schore Model")
Random Effects Lagrange Multiplier Test for the Math Schore Model
statistic p.value method alternative
Inf 0 Lagrange Multiplier Test - (Honda) for unbalanced panels significant effects

The results shows the p.value of 0, meaning that we reject the null hypothesis of zero variance in individual-specific errors. Hence, heterogeneity among schools should be significant.

# fixed effects using plm "within" method
math.random <- plm(mathscore ~ small + aide + 
                           tchexper + boy + white_asian,
                          #+ factor(schid),
                          data = star,
                         model="random",
                         random.method = "swar",
                         index=c("schid"))

School Fixed effects Model: Classroom Features effects on Math Scores

kable(tidy(math.random), caption="School Fixed effects Model: Classroom Features effects on Math Scores")
School Fixed effects Model: Classroom Features effects on Math Scores
term estimate std.error statistic p.value
(Intercept) 466.5626599 3.0802179 151.4706687 0.0000000
small 9.3011197 1.3964264 6.6606588 0.0000000
aide 0.4852552 1.3483438 0.3598898 0.7189295
tchexper 0.4373365 0.1076016 4.0644060 0.0000482
boy -6.7140946 1.1134201 -6.0301540 0.0000000
white_asian 22.4403265 2.1529340 10.4231372 0.0000000

1.e) Using the t-test statistic in (15.37) from Principles and a 5% significance level, test whether there are any significant differences between the fixed effects and random effects estimates of the coefficients on SMALL, AIDE, TCHEXPER, and WHITE_ASIAN. What are the implications of the test outcomes? What happens in we apply the test to the fixed and random estimates of the coefficient on BOY?

The t statistic from 15.37 is:
t statistic 15.37

t statistic 15.37

Where the parameter of interest is \(b_k\); denote the fixed effects estimate is\(b_{FE,k}\) and the random effects estimate as \(b_{RE,k}\).

# set the results as data tables for easy indexing
library(data.table)
cfs <- merge(setDT(tidy(math.fixed.within)), setDT(tidy(math.random)), all.y = T, by="term")
kable(cfs)
term estimate.x std.error.x statistic.x p.value.x estimate.y std.error.y statistic.y p.value.y
(Intercept) NA NA NA NA 466.5626599 3.0802179 151.4706687 0.0000000
aide 0.5268855 1.3491170 0.390541 0.6961512 0.4852552 1.3483438 0.3598898 0.7189295
boy -6.6312110 1.1133604 -5.956033 0.0000000 -6.7140946 1.1134201 -6.0301540 0.0000000
small 9.3495834 1.3970221 6.692509 0.0000000 9.3011197 1.3964264 6.6606588 0.0000000
tchexper 0.4201512 0.1083858 3.876441 0.0001072 0.4373365 0.1076016 4.0644060 0.0000482
white_asian 23.6509406 2.3108789 10.234608 0.0000000 22.4403265 2.1529340 10.4231372 0.0000000 
# columns with X are for Fixed Effects, and Y is for Random Effects results.
# use 15.37 to calculate the F statistic per term:
t_stats = (data.frame((cfs$estimate.x - cfs$estimate.y) / (((cfs$std.error.x)^2 - (cfs$std.error.y)^2)^(1/2))))
t_stats = cbind(cfs$term, t_stats)
names(t_stats) = c("term", "t_statistic")
t_stats['t>1.96?'] <- t_stats$t_statistic > 1.96
kable(t_stats)
term t_statistic t>1.96?
(Intercept) NA NA
aide 0.911566 FALSE
boy NaN NA
small 1.188071 FALSE
tchexper -1.320459 FALSE
white_asian 1.441783 FALSE

This shows the t statistics obtained from the calculation given in 15.37. Using the standard 5% large sample critical value of 1.96, we cannot reject the hypothesis that any of the estimators yield identical results.

Our conclusion is that the random effects estimator is inconsistent, and that we should use the fixed effects estimator, or should attempt to improve the model specification. The null hypothesis will be rejected for any reason that makes the two sets of estimates differ, including a misspecified model. There may be nonlinearities in the relationship we have not captured with our model, and other explanatory variables may be relevant.

In fact, more commonly, we would use the Hausman test to test this hypothesis, and I will do this below. The Hausman test tests for the null hypothesis that the individual random effects are exogenous.

Random effects estimator are reliable under the assumption that individual characteristics (heterogeneity) are exogenous, that is, they are independent with respect to the regressors in the random effects equation. The We can therefore use the Hausman test for endogeneity, with the null hypothesis that individual random effects are exogenous. The test function phtest() compares the fixed effects and the random effects models; the next code chunk performs the Hausman endogeneity test.

kable(tidy(phtest(math.fixed.within, math.random)), caption = 
 "Hausman endogeneity test for the random effects Math Scores model")
Hausman endogeneity test for the random effects Math Scores model
statistic p.value parameter method alternative
25.3505 0.0001192 5 Hausman Test one model is inconsistent

The Hausman Test shows a very low (highly significant) p-value of the test, which indicates that the null hypothesis saying that the individual random effects are exogenous is rejected, which makes the random effects equation inconsistent. In this case the fixed effects model is the correct solution.

Question 2: Instrumental Variables

A consulting firm run by Mr. John Chardonnay is investigating the relative efficiency of wine production at 75 California wineries. John sets up the production function:
\(Q=β_1+β_2MGT+β_3CAP+β_4LAB+e\)
where Q is an index of wine output for a winery, taking into account both quantity and quality, MGT is a variable that reflects the efficiency of management, CAP is an index of capital input, and LAB is an index of labor input. Because he cannot get data on management efficiency, John collects observations on the number of years of experience (XPER) of each winery manager and uses that variable in place of MGT. The 75 observations are stored in chard.dat.

2.a) Estimate the revised equation using least squares and comment on the results.

First let’s look at the data

load("./data/chard.rda")
kable(summary(chard))
q xper cap lab age
Min. : 0.0788 Min. : 3.00 Min. :-0.8358 Min. : 0.0068 Min. :25.00
1st Qu.: 6.5458 1st Qu.:11.00 1st Qu.: 5.4192 1st Qu.: 6.7858 1st Qu.:34.00
Median : 9.8072 Median :14.00 Median : 7.5045 Median : 9.1171 Median :43.00
Mean : 9.6344 Mean :13.88 Mean : 7.8347 Mean :10.0467 Mean :42.95
3rd Qu.:12.3151 3rd Qu.:17.00 3rd Qu.:10.3621 3rd Qu.:13.3645 3rd Qu.:52.00
Max. :19.1759 Max. :27.00 Max. :18.7152 Max. :23.8201 Max. :63.00 
lm1 <- lm(q ~ xper + cap + lab, data=chard)
summary(lm1)
## 
## Call:
## lm(formula = q ~ xper + cap + lab, data = chard)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3447 -1.6842 -0.1289  1.3112  9.4533 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.76226    1.05535   1.670 0.099354 .  
## xper         0.14684    0.06343   2.315 0.023517 *  
## cap          0.43796    0.11756   3.725 0.000388 ***
## lab          0.23916    0.09980   2.396 0.019195 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.756 on 71 degrees of freedom
## Multiple R-squared:  0.5616, Adjusted R-squared:  0.543 
## F-statistic: 30.31 on 3 and 71 DF,  p-value: 9.986e-13

The R squared is 0.56 (or adjusted to 0.543), meaning that this linear model can only explain about half of the variability of the data around its mean, which is considered not a good fit, so we shouldn’t count much on the estimates to reflect reality. Still, higher experience, higher capital and more labor all seem to increase the output (intuitively), and all of significant at the 5% level, whereas captial gains seem to be much more significant even at the 1% level with a p-value of 0.0004.

2.b) Find corresponding interval estimates for wine output at wineries that have the sample average values for labor and capital and have managers with 10, 20 and 30 years of experience.

setDT(chard)
chard.sample <- chard[, .(xper=c(10,20,30), cap=mean(chard$cap), lab=mean(chard$lab))]
predictions <- as.data.frame(cbind(c(10,20,30), predict(lm1, chard.sample, interval = "predict")))
names(predictions)[1] <- "xper"
kable(predictions)
xper fit lwr upr
10 9.064691 3.510797 14.61858
20 10.533072 4.947011 16.11913
30 12.001454 6.105526 17.89738

The intervals seem reasonable, somewhat wide for detecting meaningful differences between 10,20 and 30 years of experience (for example, any quantity between 7-14 could have come from either of them at the 95% confidence level).

2.c) John is concerned that the proxy variable XPER might be correlated with the error term. He decides to do a Hausman test, using the manager’s age (AGE) as an instrument for XPER. Comment on the outcome of the Hausman test.

I am incorporating the Hausman test together with the Insturemental Variables regression at the next step.

2.d) Use the instrumental variables estimator to estimate the equation \(Q=1+2XPER+3CAP+4LAB+e\) with AGE, CAP, and LAB as the instrumental variables. Comment on the results and compare them with those obtained in part (a).

We could use the 2SLS regression method, where we first predict the insterumented variable with the instrument and use those predictions as the instruements for the missing variable; however that would still yeild incorrect standard errors. Luckily, R has a package that estimates instrumental variables with correct errors. To test this, let’s specify an instrumental variables regression

# in R ivyreg function, use the original regresssors | replaced (instrumented) regressors
iv.age.xper <- ivreg(q ~ xper + cap + lab | 
                         age  + cap + lab,
                     data=chard)
summary(iv.age.xper, diagnostics=TRUE)
## 
## Call:
## ivreg(formula = q ~ xper + cap + lab | age + cap + lab, data = chard)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -7.40413 -2.17750 -0.09044  2.28339 10.62769 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -2.4867     2.7226  -0.913   0.3642  
## xper          0.5121     0.2205   2.323   0.0231 *
## cap           0.3321     0.1545   2.150   0.0349 *
## lab           0.2400     0.1209   1.985   0.0510 .
## 
## Diagnostic tests:
##                  df1 df2 statistic p-value   
## Weak instruments   1  71     9.814 0.00252 **
## Wu-Hausman         1  70     4.830 0.03127 * 
## Sargan             0  NA        NA      NA   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.338 on 71 degrees of freedom
## Multiple R-Squared: 0.3568,  Adjusted R-squared: 0.3296 
## Wald test: 21.24 on 3 and 71 DF,  p-value: 6.307e-10

The results for the quantity IV model are as follows:

  1. Weak instruments test: rejects the null, meaning that at least one instrument is strong
  2. (Wu-)Hausman test for endogeneity: barely rejects the null that the variable of concern is uncorrelated with the error term, indicating that \(xper\) is marginally endogenous at the 5% significance level (but not at the 1%).
  3. Sargan overidentifying restrictions: does not reject the null, meaning that the extra instruments are valid (are uncorrelated with the error term).
Comparing results
stargazer(lm1, iv.age.xper, type="text")
## 
## =================================================================
##                                       Dependent variable:        
##                               -----------------------------------
##                                                q                 
##                                        OLS           instrumental
##                                                        variable  
##                                        (1)               (2)     
## -----------------------------------------------------------------
## xper                                 0.147**           0.512**   
##                                      (0.063)           (0.220)   
##                                                                  
## cap                                  0.438***          0.332**   
##                                      (0.118)           (0.154)   
##                                                                  
## lab                                  0.239**            0.240*   
##                                      (0.100)           (0.121)   
##                                                                  
## Constant                              1.762*            -2.487   
##                                      (1.055)           (2.723)   
##                                                                  
## -----------------------------------------------------------------
## Observations                            75                75     
## R2                                    0.562             0.357    
## Adjusted R2                           0.543             0.330    
## Residual Std. Error (df = 71)         2.756             3.338    
## F Statistic                   30.312*** (df = 3; 71)             
## =================================================================
## Note:                                 *p<0.1; **p<0.05; ***p<0.01

Comparing the results of the OLS versus the Instrumental Variable models, we can see that: * The effect of Experience (xper) is much larger (4 times larger) in the IV model instrumenting it with age * The effect of Capital is smaller in the IV model (~75% of its value before). * The effect of Labor remains about the same. * The standard errors are greater in the IV model; espeically for the instrumented variable and the intercept. * However, the R2 and Adjusted R2 are considerably smaller for the IV model; from ~0.5 to ~0.3.

2.e) (Using the IV equation) Find corresponding interval estimates for wine output at wineries that have the sample average values for labor and capital and have managers with 10, 20 and 30 years of experience.

Using the IV regression model, I build a synthetic small dataset with the desired values to predict for, and predict the the result given those values. Now, to report the confidence intervals, I will use the corrected standard errors given from the ivreg model, as reported in the previous table (comparison table).

setDT(chard)
chard.sample <- chard[, .(xper=c(10,20,30), cap=mean(chard$cap), lab=mean(chard$lab), age=mean(chard$age))]
preds <- predict(iv.age.xper, chard.sample, interval = "predict")
preds.intervals <- as.data.frame(cbind( c(10,20,30), preds, preds - 1.96*iv.age.xper$sigma^2, preds+1.96*iv.age.xper$sigma^2))
names(preds.intervals) <- c("xper", "q_estimate", "lower", "upper")
kable(preds.intervals)
xper q_estimate lower upper
10 7.647467 -14.194997 29.48993
20 12.768487 -9.073976 34.61095
30 17.889508 -3.952956 39.73197

Comparing these interval estimates with those obtained in part (2b), these intervals are MUCH wider - more than twice as wide. The intervals seem much less relevant for interpretation, both because of their mangitude and the resulting nonsensicle result of yielding a negative amount of wine barrels (which makes no since). They present a very wide overlap - between -3 and 29, meaning for example any quantity in the range could have come from either of them at the 95% confidence level).