Time Trends and FEs

Drew Dimmery

March 14, 2014

Structure

  • Quick Homework Talk
  • Abadie's \(\kappa\)
  • Fixed effects
  • Time Trends

On CIA

  • There seemed to be some confusion on CIA.
  • The "traditional" method of choosing a covariate set is something like:
    • Find the most significant or strongest predictors that you have
    • Include those
    • Don't include the weak predictors
    • Include things that other studies have included
  • This isn't what we want, though.
  • CIA is a story about conditional randomization.
  • You don't justify based on raw correlations.
  • Because independence implies zero correlation
  • But not vice versa
  • Instead, you should think about where you might be seeing random variation, and then adjust your conditioning strategy to ensure that you only work with that variation.

More IV Stuff

  • We're going to be looking at Ananat (2011) in AEJ
  • This study looks at the effect of racial segregation on economic outcomes.
  • Outcome: Poverty rate & Inequality (Gini index)
  • Treatment: Segregation
  • Instrument: "railroad division index"
  • Main covariate of note: railroad length in a town
  • I'm dichotomizing treatment and instrument for simplicity.
  • And my outcomes are for the Black subsample
require(foreign)
d <- read.dta("aej_maindata.dta")
d$herf_b <- with(d, ifelse(herf >= quantile(herf, 0.5), 1, 0))
d$dism1990_b <- with(d, ifelse(dism1990 >= quantile(dism1990, 0.5), 1, 0))
first.stage <- lm(dism1990 ~ herf + lenper, d)
first.stage.b <- lm(dism1990_b ~ herf_b + lenper, d)
require(AER)
gini.iv <- ivreg(lngini_b ~ dism1990 + lenper, ~herf + lenper, d)
gini.iv.b <- ivreg(lngini_b ~ dism1990_b + lenper, ~herf_b + lenper, d)
pov.iv <- ivreg(povrate_b ~ dism1990 + lenper, ~herf + lenper, d)
pov.iv.b <- ivreg(povrate_b ~ dism1990_b + lenper, ~herf_b + lenper, d)

Base Results

round(summary(first.stage)$coefficients[2, ], 3)
##   Estimate Std. Error    t value   Pr(>|t|) 
##      0.357      0.081      4.395      0.000
round(summary(first.stage.b)$coefficients[2, ], 3)
##   Estimate Std. Error    t value   Pr(>|t|) 
##      0.372      0.083      4.481      0.000
round(summary(gini.iv)$coefficients[2, ], 3)
##   Estimate Std. Error    t value   Pr(>|t|) 
##      0.875      0.302      2.895      0.005
round(summary(gini.iv.b)$coefficients[2, ], 3)
##   Estimate Std. Error    t value   Pr(>|t|) 
##      0.211      0.081      2.615      0.010
round(summary(pov.iv)$coefficients[2, ], 3)
##   Estimate Std. Error    t value   Pr(>|t|) 
##      0.258      0.144      1.798      0.075
round(summary(pov.iv.b)$coefficients[2, ], 3)
##   Estimate Std. Error    t value   Pr(>|t|) 
##      0.059      0.039      1.543      0.125

Abadie's \(\kappa\)

  • Recall from the lecture that we can use a weighting scheme to calculate statistics on the compliant population.
  • \(E[g(Y,D,X)|D_1 > D_0] = {1 \over p(D_1>D_0)} E[\kappa g(Y,D,X)]\)
  • \(\kappa = 1 - {D_i (1-Z_i) \over p(Z_i =0|X)} - {(1-D_i)Z_i \over p(Z_i =1|X)}\)
  • \(E[\kappa|X] = E[D_1 -D_0|X] = E[D|X,Z=1] - E[D|X,Z=0]\)
  • Take \(w_i = {\kappa_i \over E[D_{1}-D_{0}|X_i]}\)
  • Use this in calculating any interesting statistics (means, variance, etc)
  • This let's you explore the units composing your LATE.
getKappaWt <- function(D, Z) {
    pz <- mean(Z)
    pcomp <- mean(D[Z == 1]) - mean(D[Z == 0])
    if (pcomp < 0) 
        stop("Assuming p(D|Z) > .5")
    kappa <- 1 - D * (1 - Z)/(1 - pz) - (1 - D) * Z/pz
    # Note that pcomp = mean(kappa)
    kappa/pcomp
}
w <- with(d, getKappaWt(D = dism1990_b, Z = herf_b))
varlist <- c("closeness", "area1910", "ctyliterate1920", "hsdrop_b", "manshr", 
    "ctymanuf_wkrs1920", "ngov62")
samp.stats <- sapply(varlist, function(v) mean(d[, v], na.rm = TRUE))
comp.stats <- sapply(varlist, function(v) weighted.mean(d[, v], w, na.rm = TRUE))

Examine Complier Statistics

summary(w)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   -2.51   -2.43    2.47    1.00    2.47    2.47
rbind(sample = samp.stats, compliers = comp.stats)
##           closeness area1910 ctyliterate1920 hsdrop_b manshr
## sample       -362.4    14626          0.9585   0.2516 0.1892
## compliers    -299.1    18013          0.9515   0.2424 0.2110
##           ctymanuf_wkrs1920 ngov62
## sample               0.4619  55.55
## compliers            0.4266  83.65

New Example Data

  • Now we're going to look at Gentzkow (2006)
  • This study examined how the introduction of television had an effect on voter turnout
  • Outcome: turnout in national elections
  • Treatment: Years since the introduction of TV
  • Identification comes from the idea that there are random variations in the timing of this introduction.
  • It uses local variations in this introduction.
  • Where locality is defined as a within effect.

Rough and Dirty Sweeping Function

  • plm also does this, but I like to be able to use (most of) the functions written for lm objects.
  • This assumes that you'll be clustering at the same level as your swept effects.
sweeplm <- function(formula, dat, ind) {
    newd <- model.matrix(formula, model.frame(formula, dat, na.action = na.pass))
    newd <- newd[, -1]
    newd <- cbind(dat[, as.character(formula[2])], newd)
    ok <- complete.cases(newd)
    newd <- newd[ok, ]
    ind <- ind[ok]
    newd <- apply(newd, 2, function(x) unlist(tapply(x, ind, function(z) z - 
        mean(z, na.rm = TRUE))))
    ind <- sort(ind)
    list(lm(newd[, 1] ~ newd[, -1] - 1, as.data.frame(newd)), newd, as.character(ind))
}

Robust SEs again

robust.se <- function(model, cluster, df.fix = FALSE) {
    require(sandwich)
    require(lmtest)
    M <- length(unique(cluster))
    N <- length(cluster)
    # The df.fix assumes the swept FEs were the clusters
    K <- model$rank + ifelse(df.fix, M, 0)
    dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
    uj <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum))
    rcse.cov <- dfc * sandwich(model, meat. = crossprod(uj)/N)
    rcse.se <- coeftest(model, rcse.cov)
    return(rcse.se)
}

Fixed Effects

  • Let's demonstrate some of the benefits of sweeping out fixed effects.
  • Algebraic equivalence between these various methods.
    • De-mean variables by group
    • Apply the sweep matrix:
      • \(\dot{Q} = (I_T - {1\over T}\iota_T\iota_T') \otimes I_N\)
      • \(\beta = (X_{tv}' \dot{Q} X^{tv})^{-1} X_{tv}' \dot{Q} Y\)
    • Just throw the dummies in the regression.
dat <- read.csv("gentzkow_fixed.csv")
dat$reg5year <- paste0(dat$regions5, dat$year)
dat <- subset(dat, year >= 1930)
base.lm <- lm(turnout ~ yearsoftv + advote, dat)
simple.fe <- lm(turnout ~ yearsoftv + advote + factor(reg5year), dat)
dummy.time <- system.time(lm(turnout ~ yearsoftv + advote + factor(reg5year), 
    dat))
# all.fe <-
# lm(turnout~yearsoftv+advote+factor(reg5year)+factor(stcounty),dat)
simple.swept.fe <- sweeplm(turnout ~ yearsoftv + advote, dat = dat, ind = dat$reg5year)
sweep.time <- system.time(sweeplm(turnout ~ yearsoftv + advote, dat = dat, ind = dat$reg5year))
swept.fe <- sweeplm(turnout ~ yearsoftv + advote + factor(reg5year), dat = dat, 
    ind = dat$stcounty)
rbind(dummy.time, sweep.time)
##            user.self sys.self elapsed user.child sys.child
## dummy.time     0.882    0.098   0.983          0         0
## sweep.time     0.642    0.054   0.732          0         0

FE Results

base.est <- robust.se(base.lm, dat$stcounty[-base.lm$na.action])[2, ]
## Note: no visible binding for global variable 'lmtest'
simple.est <- robust.se(simple.fe, dat$reg5year[-simple.fe$na.action])[2, ]
simple.swept.est <- robust.se(simple.swept.fe[[1]], simple.swept.fe[[3]], df.fix = TRUE)[1, 
    ]
swept.est <- robust.se(swept.fe[[1]], swept.fe[[3]], df.fix = TRUE)[1, ]
all.ests <- rbind(base.est, simple.est, simple.swept.est, swept.est)
rownames(all.ests) <- c("base", "simple", "simple.swept", "full.swept")
all.ests
##              Estimate Std. Error t value  Pr(>|t|)
## base           0.2475    0.01309  18.915 1.497e-79
## simple        -0.1729    0.05990  -2.886 3.899e-03
## simple.swept  -0.1729    0.05990  -2.886 3.899e-03
## full.swept    -0.3863    0.05107  -7.564 3.969e-14

Splines

  • You can also do this with splines, which have a slightly more non-parametric flavor.
require(splines)
spline.trend <- sweeplm(turnout ~ yearsoftv + advote + ns(trend, 3), dat = dat, 
    ind = dat$stcounty)
spline.est <- robust.se(spline.trend[[1]], spline.trend[[3]], df.fix = TRUE)[1, 
    ]
all.ests <- rbind(all.ests, spline.est)
rownames(all.ests)[9] <- "splines"
all.ests
##               Estimate Std. Error t value  Pr(>|t|)
## base            0.2475    0.01309  18.915 1.497e-79
## simple         -0.1729    0.05990  -2.886 3.899e-03
## simple.swept   -0.1729    0.05990  -2.886 3.899e-03
## full.swept     -0.3863    0.05107  -7.564 3.969e-14
## linear          0.1476    0.01587   9.304 1.392e-20
## linear.region   0.1985    0.01503  13.207 9.192e-40
## cubic          -0.2946    0.04593  -6.413 1.434e-10
## cubic.region   -0.1188    0.04229  -2.809 4.966e-03
## splines        -0.4105    0.04949  -8.294 1.122e-16