Chapter 03: Causal Inference

This note is a high-fidelity Markdown migration of the Causal Inference chapter from the LaTeX source.

Parent map: index Prerequisites: probability-and-mathstats, linear-regression

Concept map

flowchart TD
  A[Potential Outcomes] --> B[Randomized Experiments]
  A --> C[Selection on Observables]
  C --> D[Regression Adjustment]
  C --> E[Matching / IPW / AIPW]
  A --> F[Instrumental Variables]
  F --> G[LATE / Compliers]
  F --> H[MTE]
  A --> I[RD / RDD]
  A --> J[DiD / Panel]
  J --> K[Staggered Adoption]
  J --> L[Synthetic Control]
  A --> M[Decomposition Methods]
  A --> N[Causal DAGs]

Foundations, Experiments

Potential Outcomes

Exposition from

$Y_i$ is the observed outcome, $D_i$ is the treatment with levels $d\in \mathcal{D}$,

potential outcomes denoted $Y_{di},Y_i^d,Y_i(d)$ (interchangeably).

$$Y_i^{\text{obs}}=Y_i(D_i)=\{\begin{array}{ll}{Y}_{1i},& D_i=1,\\ Y_{0i},& D_i=0.\end{array}$$

Equivalently, we have the switching equation

$$Y_i=D_i{Y}_{1i}+(1-D_i)Y_{0i}=Y_{0i}+(Y_{1i}-Y_{0i})D_i.$$

This encodes what is known as the causal-consistency assumption (/ SUTVA).

Generally, define a potential outcome

$Y_i^d=\phi (d,{\mathbf{X}}_i,{\mathbf{U}}_i)$

where ${\mathbf{X}}_i$ is a vector of observed covariates and ${\mathbf{U}}_i$ is a vector of unobservables, and $\phi$ is an unknown measurable function. Typically, we are interested in non-parametric identification of $\phi$ or some features of it.

Given a population of $n$ units, the assignment mechanism is a row-exchangeable function $\mathbf{Pr}(\mathbf{D}\mid \mathbf{X},\mathbf{Y}(0),\mathbf{Y}(1))$ taking values on $[0,1]$ and satisfying

$$\sum _{\mathbf{d}\in \{0,1{\}}^N}\mathbf{Pr}(\mathbf{D}=\mathbf{d}\mid \mathbf{X},\mathbf{Y}(0),\mathbf{Y}(1))=1.$$

A unit level assignment probability for unit $i$ is

$$p_i(\mathbf{X},\mathbf{Y}(0),\mathbf{Y}(1))=\sum _{\mathbf{d}:d_i=1}\mathbf{Pr}(\mathbf{D}=\mathbf{d}\mid \mathbf{X},\mathbf{Y}(0),\mathbf{Y}(1)).$$

A finite population propensity score is

$$e(x)=\frac{1}{N(x)}\sum _{i:{\mathbf{X}}_i=\mathbf{x}}{p}_i(\mathbf{X},\mathbf{Y}(0),\mathbf{Y}(1)).$$

where $N(x)=\#\{i=1,\dots ,N:{\mathbf{X}}_i=\mathbf{x}\}$ is the number of units in each stratum defined by ${\mathbf{X}}_i=\mathbf{x}$.

is a row-exchangeable function of potential outcomes, treatment assignment, and covariates.

$\tau =\tau (\mathbf{Y}(0),\mathbf{Y}(1),\mathbf{X},\mathbf{d})$

$\mathbf{Y}(0),\mathbf{Y}(1)$ are $n-$ vectors of potential outcomes, $\mathbf{X}$ is a $n\times p$ covariate matrix, and $\mathbf{d}$ is an assignment vector.

The most intuitive estimand is a $n-$ vector $\mathbf{\tau }=\mathbf{Y}(1)-\mathbf{Y}(0)$. This is impossible to estimate because of the FPCI, so we instead use summaries, such as its sample average, or subgroup averages.

We never see both potential outcomes for any given unit.

Decompositions of Observed Differences:

$$\underset{\text{observed difference}}{\underbrace{\mathbb{E}[Y_i\mid D_i=1]-\mathbb{E}[Y_i\mid D_i=0]}}=\underset{\text{ATT}}{\underbrace{\mathbb{E}[Y_{1i}\mid D_i=1]-\mathbb{E}[Y_{0i}\mid D_i=1]}}+\underset{\text{selection bias}}{\underbrace{\mathbb{E}[Y_{0i}\mid D_i=1]-\mathbb{E}[Y_{0i}\mid D_i=0]}}.$$ $$\underset{\text{observed difference}}{\underbrace{\mathbb{E}[Y_i\mid D_i=1]-\mathbb{E}[Y_i\mid D_i=0]}}=\underset{\text{ATE}}{\underbrace{\mathbb{E}\left[Y_1\right]-\mathbb{E}\left[Y_0\right]}}+\underset{\text{selection bias}}{\underbrace{\mathbb{E}[Y_{0i}\mid D_i=1]-\mathbb{E}[Y_{0i}\mid D_i=0]}}+\underset{\text{heterogeneous treatment bias}}{\underbrace{(1-\pi )\left(ATT-ATU\right)}}.$$

where $\pi =\mathbb{E}\left[D\right]$ is the share of the sample treated.

$(Y_{1i},Y_{0i})\perp \perp D_i$

This is a Missing Completely at Random (MCAR) assumption on potential outcomes.

Writing outcomes generated by the switching regression assumes that potential outcomes for any unit do not vary with the treatment assigned to other units. In practice, this is equivalent to a no spillovers assumption.

$(Y_{1i},Y_{0i})\perp \perp {\mathbf{D}}_{-i}$

Equivalently, let $\mathbf{D}$ denote a treatment vector for $N$ units, and $\mathbf{Y}(\mathbf{D})$ be the potential outcome vector that would be observed if was based on allocation $\mathbf{D}$. Then, SUTVA requires that for allocations $\mathbf{D},{\mathbf{D}}^{\prime }$,

$Y_i(\mathbf{D})=Y_i({\mathbf{D}}^{\prime })\text{ if}{D}_i=D_i^{\prime }$

Intuitively, SUTVA ensures that the ‘science table’ (Imbens & Rubin 2015) has 2 columns for the two potential outcomes as opposed to $2^n$ (number of potential outcomes with arbitrary interference).

Treatment Effects

Estimands

• ${\tau }_{\text{ATE}}:=\mathbb{E}\left(Y_{1i}-Y_{0i}\right)$

• ${\tau }_{\text{ATT}}:=\mathbb{E}\left(Y_{1i}-Y_{0i}|D_i=1\right)=\mathbb{E}\left[Y_{1i}|D_i=1\right]-\mathbb{E}\left[Y_{0i}|D_i=1\right]$

Under randomisation, ${\tau }_{ATE}={\tau }_{ATT}$, since the treated are a random sample of the population. Under weak(er) assumption of $Y_{0i}\perp \perp D_i$, only ${\tau }_{ATT}$ is identified.

Difference in Means

$$\hat{\tau }=\frac{1}{N_1}\sum _{i=1}^N{D}_i{Y}_i-\frac{1}{N_0}\sum _{i=1}^N(1-D_i)Y_i.$$

Variance of Difference in means estimator is given by

$$\mathrm{Var}({\hat{\tau }}_{\text{DiM}})=\frac{S_0^2}{N_0}+\frac{S_1^2}{N_1}-\frac{S_{01}^2}{N}.$$

where $S_0,S_1$ are sample variances of $Y^0,Y^1$ respectively, and $S_{01}$ is the variance of the unit level treatment effect

$$\frac{1}{N-1}\sum _{i=1}^N(Y_i(1)-Y_i(0)-\tau ).$$

This is not identifiable because of the last term. If the treatment effect is constant in the population, the last term is zero.

A (conservative) variance estimator is given by

$$\hat{\mathbb{V}}({\hat{\tau }}_{\text{DiM}})=\frac{{\hat{\sigma }}_1^2}{N_1}+\frac{{\hat{\sigma }}_0^2}{N_0}.$$

where

$${\hat{\sigma }}_{(d)}^2=\frac{1}{N_d-1}\sum _{i:d_i=d}(Y_i-{\overline{Y}}_d{)}^2,d\in \{0,1\}.$$

These variance estimates can be used to construct 95% confidence intervals

$${\mathcal{C}}_{0.95}(\tau )=\left(\hat{\tau }-1.96\sqrt{\hat{\mathbb{V}}},\hat{\tau }+1.96\sqrt{\hat{\mathbb{V}}}\right).$$

Regression Adjustment

$$Y_i=\alpha +{\tau }_{\text{REG}}{D}_i+{\eta }_i={\overline{Y}}_0+({\overline{Y}}_1-{\overline{Y}}_0)D_i+{\eta }_i.$$

• $\alpha =\mathbb{E}[Y_{0i}]$

• $\tau =\mathbb{E}\left[Y_{1i}-Y_{0i}\right]$

• ${\eta }_i=Y_{0i}-\mathbb{E}(Y_{0i})$ [extra terms above come from allowing for heterogeneous TEs]

Selection bias: $\mathrm{Cov}\left(D_i,{\eta }_i\right)\ne 0$

Suppose $50$ percent of the population gets the treatment. Let ${\mathbf{X}}_i^{\prime }=\left[\begin{array}{cc}{D}_i& 1\end{array}\right]$. Then,

$$\hat{\tau }={\bar{Y}}_T-{\bar{Y}}_C,\hat{\alpha }={\bar{Y}}_C.$$

Generalise to $p$ fraction treated

VCV under homoscedasticity

$$\mathrm{Var}(\hat{\tau })\approx \frac{{\hat{\sigma }}^2}{p(1-p)N},{\hat{\sigma }}^2=\frac{{\sum }_i{\hat{u}}_i^2}{N-2}.$$

VCV under heteroskedasticity

$${\widehat{\mathrm{Var}}}_{HC}(\hat{\tau })=\frac{1}{N^2}\sum _{i=1}^N\left(\frac{D_i{\hat{u}}_i^2}{{\hat{p}}^2}+\frac{(1-D_i){\hat{u}}_i^2}{(1-\hat{p}{)}^2}\right),\hat{p}=\frac{N_1}{N}.$$

Including controls:

$Y_i=\alpha +\tau D_i+{\mathbf{X}}_i^{\prime }\mathbf{\beta }+{\eta }_i$

Corrects for chance covariate imbalances, improves precision by removing variation in outcome accounted for by pre-treatment characteristics.

Freedman (2008) Critique

Regression of the form $Y_i=\alpha +{\tau }_{reg}{D}_i+{\beta }_1{X}_i+{ϵ}_i$

${\hat{\tau }}_{reg}$ is consistent for ATE but has small sample bias (unless model is true); bias is on the order of $1/n$

${\hat{\tau }}_{reg}$ precision does not improve through the inclusion of controls; including controls is harmful to precision if more than $3/4$ units are assigned to one treatment condition

Recommends fitting

$$Y_i=\alpha +{\tau }_{\text{lin}}{D}_i+{\beta }_0^{\top }(X_i-\bar{X})+{\beta }_1^{\top }{D}_i(X_i-\bar{X})+{ϵ}_i$$

Where the two potential outcomes are stipulated to follow

$$\hat{Y}(1)={\bar{Y}}_1+(\bar{X}-{\bar{X}}_1{)}^{\top }({\hat{\beta }}_0+{\hat{\beta }}_1),\hat{Y}(0)={\bar{Y}}_0+(\bar{X}-{\bar{X}}_0{)}^{\top }{\hat{\beta }}_0.$$

which has same small sample bias, but cannot hurt asymptotic precision even if the model is incorrect and will likely increase precision if covariates are predictive of the outcomes.

Randomisation Inference

sharp null: $Y_{1i}=Y_{0i}\forall i$. Implies $H_0:\mathbb{E}\left[Y_i\right]=\mathbb{E}\left[Y_0\right];H_1:\mathbb{E}\left[Y_i\right]\ne \mathbb{E}\left[Y_0\right]$.

To test sharp null, set $Y_1=Y_0$ for all units and re-randomize treatment. Complete randomisation of $2N$ units with $N$ treated. $\left(\genfrac{}{}{0ex}{}{2N}{N}\right)$ assignment vectors. $P$ value can be as small as $1/\left(\genfrac{}{}{0ex}{}{2N}{N}\right)$.

$\Omega$ is the full set of randomisation realisations, and $\omega$ is an element in the set (drawn either under complete randomization or binomial randomization), with associated probability $1/\left(\genfrac{}{}{0ex}{}{2N}{N}\right)$

One sided P-value : $\mathbf{Pr}\left(\left(\hat{\alpha }(\omega )\ge {\hat{\tau }}_{\text{ATE}}\right)\right)$

Blocking

Stratify randomisation to ensure that groups start out with identical observable characteristics on blocked factors.

$V\left[{\hat{\tau }}_{BR}\right]<V\left[{\hat{\tau }}_{CR}\right]$ if $\frac{SS{R}_{{\hat{\epsilon }}^{\ast }}}{n-k-1}<\frac{SS{R}_{\hat{\epsilon }}}{n-2}$ where $\widehat{ϵ}$ and $\widehat{{ϵ}^{\ast }}$ are errors from specification omitting and including block dummies respectively.

For $J$ blocks,

Point estimate ${\hat{\tau }}_B={\sum }_{j=1}^J\frac{N_j}{N}{\hat{\tau }}_j$

Variance Randomisations within each block are independent, so the variances are simple means (with squared weights). $\mathrm{Var}\left({\hat{\tau }}_B\right)={\sum }_{j=1}^J{\left(\frac{N_j}{N}\right)}^2\mathrm{Var}\left({\hat{\tau }}_j\right)$

Regression Formulation

$$y_i=\tau D_i+\sum _{j=2}^J{\beta }_j{B}_{ij}+{ϵ}_i$$

If treatment probabilities vary by block, then weight by

$$w_{ij}=\frac{D_i}{p_{ij}}+\frac{1-D_i}{1-p_{ij}}.$$

Efficiency Gains from Blocking

Complete Randomisation : $Y_i = \alpha + \tau_{CR} D_i

• \epsilon_i$

Block Randomisation: $Y_i=\alpha +{\tau }_{BR}{D}_i+{\sum }_{j=2}^J{\beta }_j{B}_{ij}+{ϵ}_i^{\ast }$

$$\mathrm{Var}({\hat{\tau }}_{CR})=\frac{{\sigma }_{\epsilon }^2}{{\sum }_{i=1}^n{\left(D_i-\overline{D}\right)}^2},{\hat{\sigma }}_{\epsilon }^2=\frac{{\sum }_{i=1}^n{\hat{\epsilon }}_i^2}{n-2}=\frac{SS{R}_{\hat{\epsilon }}}{n-2}.$$ $$\mathrm{Var}({\hat{\tau }}_{BR})=\frac{{\sigma }_{{\epsilon }^{\ast }}^2}{{\sum }_{i=1}^n{\left(D_i-\overline{D}\right)}^2\left(1-R_j^2\right)},{\hat{\sigma }}_{{\epsilon }^{\ast }}^2=\frac{{\sum }_{i=1}^n{\left({\hat{\epsilon }}_i^{\ast }\right)}^2}{n-k-1}=\frac{SS{R}_{{\hat{\epsilon }}^{\ast }}}{n-k-1}.$$

Where $R_j^2$ is the fit from regressing $D$ on all $B_j$ dummies. Since $R_j^2\approx 0$ by randomisation,

$$\mathrm{Var}({\hat{\tau }}_{BR})<\mathrm{Var}({\hat{\tau }}_{CR})⟺\frac{SS{R}_{{\hat{\epsilon }}^{\ast }}}{n-k-1}<\frac{SS{R}_{\hat{\epsilon }}}{n-2}.$$

Power Calculations

Basic idea: With large enough samples, $\mathbb{V}\left[{\bar{Y}}_1-{\bar{Y}}_0\right]\approx \frac{{\sigma }_1^2}{pN}+\frac{{\sigma }_0^2}{(1-p)N}$[where $p=N_1/N$ is the share of sample treated]. Set $p$ to minimise overall variance. Yields $p^{\ast }=\frac{{\sigma }_1}{{\sigma }_1+{\sigma }_0}$. With homoskedasticity, this is $\frac{1}{2}$ Treatment, $\frac{1}{2}$ control.

$\tau ={\mu }_1-{\mu }_0$ (effect size)

Test for $\tau >(t_{1-\kappa }+t_{\alpha /2}SE(\hat{\beta })$.

For common variance $\sigma$,

$$\pi =\mathbf{Pr}\left(|t|>1.96\right)=\Phi \left(-1.96-\frac{\tau \sqrt{N}}{2\sigma }\right)+\left(1-\Phi \left(1.96-\frac{\tau \sqrt{N}}{2\sigma }\right)\right).$$

General formula for Power with unequal variances

$$\pi =\Phi \left(-1.96+\frac{\tau }{\sqrt{\frac{{\sigma }_1^2}{pN}+\frac{{\sigma }_0^2}{(1-p)N}}}\right)+\left(1-\Phi \left(1.96-\frac{\tau }{\sqrt{\frac{{\sigma }_1^2}{pN}+\frac{{\sigma }_0^2}{(1-p)N}}}\right)\right).$$

This yields

Common variance (assumed)

$$\mathrm{MDE}(\tau )=M_{n-2}\sqrt{\frac{{\sigma }^2}{Np(1-p)}}.$$

where $M_{n-2}=t_{(1-\alpha /2)}+t_{1-\kappa }$ = Critical t-value to reject null + t-value for alternative (where $1-\kappa$) is power.

MDES (Minimum Detectable Effect Size in Standard Deviation Units):

$\text{MDES}(\tau )=M_{n-2}\sqrt{\frac{1}{Np(1-p)}}$

Multiplier $M_{n-2}$ simplifies to $1.96+0.84\approx 2.8$

$$\text{MDE}\approx (0.84+1.96)\text{SE}(\hat{\tau })\approx 2.8\text{SE}(\hat{\tau }).$$

Rearrange to get necessary sample size for any given hypothesised MDE and expected variance.

$$N=(z_{1-\kappa }+z_{\alpha /2}{)}^2\cdot \frac{1}{p(1-p)}\cdot \frac{{\sigma }^2}{{\text{MDE}}^2}.$$

MDES for Blocking

$$\text{MDES}({\tau }_{BR})=M_{n-k-1}\sqrt{\frac{1-R_B^2}{Np(1-p)}}.$$

where $R_B^2$ is the R-squared from regressing $Y$ on block dummies.

To test $H_0:\mathbb{E}\left[Y_i(1)-Y_i(0)\right]=0$ against the alternative, we look at the T Statistic

$$T=\frac{{\overline{Y}}_t^{obs}-{\overline{Y}}_c^{obs}}{\sqrt{S_y^2/N_t+S_y^2/N_c}}\approx \mathcal{N}\left(\frac{\tau }{\sqrt{{\sigma }^2/N_t+{\sigma }^2/N_c}},1\right).$$

Inverting this for size $\alpha /2$ gives us a required sample size

$$\text{Required Sample Size}=N=\frac{{\left({\Phi }^{-1}(\beta )+{\Phi }^{-1}(1-\alpha /2)\right)}^2}{(\tau /\sigma {)}^2\cdot \gamma \cdot (1-\gamma )}.$$

typically, $\beta =0.8$, $\alpha =0.05,\gamma =0.5$, so by substitution:

$$N=\frac{{\left({\Phi }^{-1}(0.8)+{\Phi }^{-1}(0.975)\right)}^2}{(\tau /\sigma {)}^2\cdot 0.5^2}.$$

Selection On Observables

typology

Regression estimators: rely on consistent estimation of ${\mu }_0(\mathbf{x}),{\mu }_1(\mathbf{x})$ Matching estimators Propensity score estimators: rely on estimation of $\pi (\mathbf{x})$ Combination methods (augmented IPW, bias-corrected Matching, etc)

Regression Anatomy / FWL

$${\beta }_k=\frac{\mathrm{Cov}\left(Y_i,{\tilde{x}}_{ki}\right)}{\mathbb{V}\left[{\tilde{x}}_{ki}\right]}.$$

where ${\tilde{x}}_{ki}$ is the residual from a regression of $x_{ki}$ on all other covariates.

If structural (long) equation is $Y_i=\alpha +\tau D_i+W_i^{\prime }\gamma +{ϵ}_i$, with $W_i$ vector of unobserved, and we estimate short $Y_i=\alpha +\rho D_i+{ϵ}_i$, then we can write the specification as $y=\tau D_i+\underset{{\nu }_i}{\underbrace{W_i^{\prime }\gamma +ϵ}}$

$\rho =\frac{\mathrm{Cov}\left(Y_i,D_i\right)}{\mathbb{V}\left[D\right]}=\tau +{\gamma }^{\prime }{\delta }_{WD}$

equivalently,

$$\mathrm{plim}{\hat{\tau }}_{\text{OLS}}=\tau +\delta \gamma =\tau +\underset{\text{Omitted Variables Bias}}{\underbrace{\mathrm{plim}\left[{\left(N^{-1}{D}^{\prime }D\right)}^{-1}{N}^{-1}{D}^{\prime }W\right]\gamma }}.$$

Coefficient in Short Regression = Coefficient in long regression + effect of omitted $\times$ regression of omitted on included. This bias can be arbitrarily large.

Identification of Treatment Effects under Unconfoundedness

• Unconfoundedness / Selection on Observables / Ignorability / Conditional Independence Assumption: $(Y_0,Y_1)\perp \perp D|\mathbf{X}$ In terms of densities, this is equivalent to the validity of the following density factorisation

$$f_{Y(d),D\mid \mathbf{X}}(y,d\mid \mathbf{x})=f_{Y(d)\mid \mathbf{X}}(y\mid \mathbf{x})f_{D\mid \mathbf{X}}(d\mid x)=f_{Y\mid D,\mathbf{X}}(y\mid d,\mathbf{x})f_{D\mid \mathbf{X}}(d\mid x).$$

• common support $0<\mathbf{Pr}\left(D=1|X\right)<1$

$$\mathbb{E}\left[Y_d\right]=\int \mathbb{E}\left[Y^d\mid \mathbf{X}=x\right]d{P}_{\mathbf{x}}(x)=\int \mathbb{E}\left[Y^d\mid D=d,\mathbf{X}=x\right]d{P}_{\mathbf{x}}(x)=\int \mathbb{E}\left[Y\mid D=d,\mathbf{X}=x\right]d{P}_{\mathbf{x}}(x).$$

The third quantity is estimable using observed data.

Estimators:

Discrete Case: $\mathbf{x}$ has finite values indexed by $k=1,\dots ,K$ with generic entry ${\mathbf{x}}_k$

$${\tau }_{\text{ATE}}=\sum _{k=1}^K\left(\mathbb{E}\left[Y\mid D=1,\mathbf{X}=x_k\right]-\mathbb{E}\left[Y\mid D=0,\mathbf{X}=x_k\right]\right)\mathbf{Pr}\left(\mathbf{X}=x_k\right).$$ $${\tau }_{\text{ATT}}=\sum _{k=1}^K\left(\mathbb{E}\left[Y\mid D=1,\mathbf{X}=x_k\right]-\mathbb{E}\left[Y\mid D=0,\mathbf{X}=x_k\right]\right)\mathbf{Pr}\left(\mathbf{X}=x_k\mid D=1\right).$$

Multi-valued and Continuous Treatments

Treatment values: $\mathcal{D}$ finite if multi-valued / $\subset \mathbb{R}$ for continuous, with corresponding dose-responses $Y_i(d)$. We are interested in dose-response function $\mu (d)=\mathbb{E}\left[Y_i(d)\right]$, and contrasts.

First define Generalised propensity score : $R:=r(d,\mathbf{x})={\mathsf{f}}_{D|\mathbf{X}}\left(d|\mathbf{x}\right)$

Assumptions: Weak unconfoundedness: $Y(d)\perp \perp D|\mathbf{X}=\mathbf{x}\forall D\in \mathcal{D}$ Conditional density overlap: $\mathsf{f}\left(D=d|\mathbf{X}=\mathbf{x}\right)>0$

Bias removal using the generalised propensity score: Estimate the conditional expectation of the outcome as a function of treatment level $d$ and GPS $R$ as $\beta (d,r)=\mathbb{E}\left[Y(d)|r(d,\mathbf{X})=r\right]=\mathbb{E}\left[Y|D=d|R=r\right]$ Estimate the dose-response function of the treatment by averaging the conditional expectation at that particular level of treatment $\mu (d)=\mathbb{E}\left[\beta (d,r(d,\mathbf{X}))\right]$

Then compute contrasts to get first derivative (MTE)

$\frac{\partial }{\partial d}\mathbb{E}\left[\mu (d)\right]$

Estimators of $\mathbb{E}\left[Y^d\right]$

which can be used to construct estimators of ATE(${\hat{\gamma }}_1-{\hat{\gamma }}_0$), ATT($({\gamma }_1-{\hat{\gamma }}_0|D=1)$, and other estimands. reference: , David Childers’ lecture notes.

Regression Adjustment

Estimate ${\mu }_d(\mathbf{x})=\mathbb{E}\left[Y\mid D=d,\mathbf{X}=x\right]$ by a nonparametric regression estimator ${\hat{\mu }}_d(\mathbf{x})$.

$${\hat{\gamma }}_d^{\text{reg}}:=\frac{1}{n}\sum _{i=1}^n{\hat{\mu }}_d({\mathbf{x}}_i).$$

Since the average predicted treated outcome for the treated equals the average predicted control outcome for controls, one ATE form is

$${\hat{\tau }}_{\text{reg}}^{\text{ATE}}=\frac{1}{n}\sum _{i=1}^n\left[D_i\left(Y_i-{\hat{\mu }}_0({\mathbf{X}}_i)\right)+(1-D_i)\left({\hat{\mu }}_1({\mathbf{X}}_i)-Y_i\right)\right].$$

SATT only requires imputation of one potential outcome:

$${\hat{\tau }}_{\text{reg}}^{\text{ATT}}=\frac{1}{N_t}\sum _{i=1}^n{D}_i\left[Y_i-{\hat{\mu }}_0({\mathbf{X}}_i)\right].$$

Inverse Propensity Weighting

Estimate propensity score $\pi (\mathbf{x})=\mathbf{Pr}\left(D=d\mid \mathbf{X}=x\right)$ by $\hat{\pi }(d\mid \mathbf{x})$. Then

$${\hat{\gamma }}_d^{\text{IPW}}:=\frac{1}{n}\sum _{i=1}^n\frac{Y_i{1}_{D_i=d}}{\hat{\pi }(d\mid {\mathbf{x}}_i)}.$$

Augmented Inverse Propensity Weighting / Combination Methods

Estimate ${\mu }_d(\mathbf{x})$ and $\pi (\mathbf{x})$, then average

$${\hat{\gamma }}_d^{\text{AIPW}}:=\frac{1}{n}\sum _{i=1}^n\left[{\hat{\mu }}_d({\mathbf{x}}_i)+\frac{\left(Y_i-{\hat{\mu }}_d({\mathbf{x}}_i)\right)1_{D_i=d}}{\hat{\pi }(d\mid {\mathbf{x}}_i)}\right].$$

Normalized outcome regression:

$${\tilde{\mu }}_1=\frac{{\hat{\mu }}_1(\mathbf{x})}{\hat{\pi }(\mathbf{x})},{\tilde{\mu }}_0=\frac{{\hat{\mu }}_0(\mathbf{x})}{1-\hat{\pi }(\mathbf{x})}.$$

Subclassification / Blocking

Weighted combination of $K$ subclasses of covariate values, which partition the population

$${\hat{\tau }}^{\text{ATE}}=\sum _{k=1}^K\left({\overline{Y}}_1^k-{\overline{Y}}_0^k\right)\cdot \left(\frac{N^k}{N}\right).$$ $${\hat{\tau }}_{\text{ATT}}=\sum _{k=1}^K\left({\overline{Y}}_1^k-{\overline{Y}}_0^k\right)\cdot \left(\frac{N_1^k}{N_1}\right).$$

Regression Adjustment

A single regression with controls $X$ is potentially problematic because of Simpson’s paradox. To account for this in a parametric setup, assume a set of iid subjects $i=1,\dots n$ we observe a tuple $(X_i,Y_i,D_i)$, comprised of

feature vector ${\mathbf{X}}_i\in {\mathbb{R}}^p$

response $Y_i\in \mathbb{R}$

treatment assignment $D_i\in \left\{0,1\right\}$

Define conditional response surfaces as

${\mu }_{(d)}(x):=\mathbb{E}\left[Y_i|{\mathbf{X}}_{\mathbf{i}}=\mathbf{x},D_i=d\right]$

First pass regression adjustment estimator (using OLS)

$$\hat{\tau }=\frac{1}{n}\sum _{i=1}^n\left[{\hat{\mu }}_{(1)}({\mathbf{X}}_i)-{\hat{\mu }}_{(0)}({\mathbf{X}}_i)\right].$$

where ${\hat{\mu }}_{(d)}(\mathbf{x})$ is obtained via OLS. This generically doesn’t work for regularised regression.

With known propensity score $\pi (\mathbf{X})$ (as in case of regression), an efficient estimator weights all estimated treatment effects ${\hat{\mu }}_1({\mathbf{X}}_i)-{\hat{\mu }}_0({\mathbf{X}}_i)$ by the propensity score:

$${\tilde{\tau }}_{\text{reg}}^{\text{ATT}}=\frac{{\sum }_{i=1}^n\pi ({\mathbf{X}}_i)\left[{\hat{\mu }}_1({\mathbf{X}}_i)-{\hat{\mu }}_0({\mathbf{X}}_i)\right]}{{\sum }_{i=1}^n\pi ({\mathbf{X}}_i)}.$$

Additional Assumptions for consistent estimate of ATE from OLS:

1. Constant treatment effects

2. Outcomes linear in X

$⟹\tau$ will provide unbiased and consistent estimates of ATE.

• 1. fails - ${\tau }_{OLS}$ is Best Linear Approximation of average causal response function $\mathbb{E}\left[Y|D=1,X\right]-\mathbb{E}\left[Y|D=0,X\right]$.
• 1. fails - ${\tau }_{OLS}$ is conditional variance weighted average of underlying $\tau$s.

Pretend there are $m$ strata of $\mathbf{X}$. Then, OLS estimates

$${\tau }_{\text{OLS}}=\sum _{k=1}^m\left(\mathbb{E}[Y\mid X=x_k,D=1]-\mathbb{E}[Y\mid X=x_k,D=0]\right){\omega }_k.$$

where the weight

$${\omega }_k=\frac{\mathbb{V}\left[D\mid X=x_k\right]\mathbf{Pr}\left(X=x_k\right)}{{\sum }_{r=1}^m\mathbb{V}\left[D\mid X=x_r\right]\mathbf{Pr}\left(X=x_r\right)}.$$

${\tau }_{\text{OLS}}$ weighs up groups where the size of the treated and untreated population are roughly equal, and weighs down groups with large imbalances in the size of these two groups.

${\tau }_{\text{OLS}}$ is true effect IFF constant treatment effects holds.

Matching

Regression estimators impute missing potential outcomes by imputing it using ${\hat{\mu }}_d({\mathbf{X}}_i)$. Matching estimators proceed by imputing the potential outcome using the observed outcome from the ‘closest’ control unit.

Define ${\ell }_m(i)$ as the index that satisfies

$$\sum _{j:d_j\ne d_i}1\left\{\Vert {\mathbf{X}}_j-{\mathbf{X}}_i\Vert \le \Vert {\mathbf{X}}_{{\ell }_m(i)}-{\mathbf{X}}_i\Vert \right\}=m.$$

So, ${\ell }_m(i)$ is the index of the unit in the opposite treatment group that is $m-$th closest to unit $i$ in terms of covariate values in terms of the norm $\Vert \cdot \Vert$. Let ${\mathcal{J}}_M(i):=\left\{{\ell }_1(i),\dots ,{\ell }_M(i)\right\}$ denote the indices of the first $M$ matches for unit $i$. Then, impute potential outcomes as

$${\hat{Y}}_i(0)=\{\begin{array}{ll}{Y}_i& \text{if }{D}_i=0,\\ \frac{1}{M}{\sum }_{j\in {\mathcal{J}}_M(i)}{Y}_j& \text{if }{D}_i=1.\end{array}$$ $${\hat{Y}}_i(1)=\{\begin{array}{ll}\frac{1}{M}{\sum }_{j\in {\mathcal{J}}_M(i)}{Y}_j& \text{if }{D}_i=0,\\ Y_i& \text{if }{D}_i=1.\end{array}$$

then, the simple matching (with replacement) estimator for ATE is

$${\hat{\tau }}_{\text{Match}}^{\text{ATE}}=\frac{1}{n}\sum _{i=1}^n\left[{\hat{Y}}_i(1)-{\hat{Y}}_i(0)\right].$$

and corresponding ATT

$${\hat{\tau }}_{\text{Match}}^{\text{ATT}}=\frac{1}{N_1}\sum _{i:D_i=1}\left(Y_i-Y_{j(i)}\right).$$

where $M=1$ corresponds with one-to-one matching and $M>1$ is many-to-one. Many-to-one matching is not $\sqrt{n}$ consistent (Abadie and Imbens (2006)) and has a bias of $O(N^{-1/k})$ where $k$ is the number of continuous covariates.

Bias-corrected (Abadie-Imbens)

$${\hat{\tau }}_{\text{ATT}}=\frac{1}{N_1}\sum _{i:D_i=1}\left[\left(Y_i-Y_{j(i)}\right)-\left({\hat{\mu }}_0(X_i)-{\hat{\mu }}_0(X_{j(i)})\right)\right].$$

Where ${\mu }_0(x)=E[Y|X=x,D=0]$ is the regression function under the control.

Metrics

• Euclidian Distance

$$\Vert X_i-X_j\Vert =ED(X_i,X_j)=\sqrt{(X_i-X_j{)}^{\prime }(X_i-X_j)}.$$

• Stata diagonal distance

$$\mathrm{StataD}(X_i,X_j)=\sqrt{(X_i-X_j{)}^{\prime }\mathrm{diag}{\left({\hat{\Sigma }}_X\right)}^{-1}(X_i-X_j)}.$$

where the normalisation factor is the diagonal element of $\hat{\Sigma}$, the estimated variance covariance matrix.

• Mahalanobis distance (scale-invariant)

$$\text{MD}(X_i,X_j)=\sqrt{(X_i-X_j{)}^{\prime }{\Sigma }^{-1}(X_i-X_j)}.$$

Where $\hat{\Sigma }$ is the variance-covariance matrix.

Matching estimators have a normal distribution in large samples provided that bias is small.

For matching without replacement,

$${\hat{\sigma }}_{\text{ATT}}^2=\frac{1}{N_T}\sum _{D_i=1}{\left(Y_i-\frac{1}{M}\sum _{m=1}^M{Y}_{jm(i)}-{\hat{\delta }}_{\text{ATT}}\right)}^2.$$

For matching with replacement,

$${\hat{\sigma }}_{\text{ATT}}^2=\frac{1}{N_T}\sum _{D_i=1}{\left(Y_i-\frac{1}{M}\sum _{m=1}^M{Y}_{jm(i)}-{\hat{\delta }}_{\text{ATT}}\right)}^2+\frac{1}{N_T}\sum _{D_i=0}\left(\frac{K_i(K_i-1)}{M^2}\right)\mathbb{V}\left[ϵ\mid X_i,D_i=0\right].$$

where $K_i$ is the number of times observation $i$ is used in a match, and the last error variance term is estimated by matching also. the bootstrap doesn’t work for matching.

PScore is a balancing score: conditioning on propensity score is equivalent to conditioning on covariates.

$$\mathbf{Pr}\left(D=1\mid Y_0,Y_1,\pi (X)\right)=\mathbf{Pr}\left(D=1\mid \pi (X)\right)=\pi (X).$$ $$\{Y_i(0),Y_i(1)\}\perp \perp D_i\mid {\mathbf{X}}_i⟺\{Y_i(0),Y_i(1)\}\perp \perp D_i\mid \pi ({\mathbf{X}}_i).$$

defines the semiparametric Efficiency Bound for ATE: the asymptotic variance of any regular estimator of $\tau$ of the population ATE obeys

$\sqrt{n}(\hat{\tau }+{\tau }^P)\stackrel{d}{\to }\mathcal{N}\left(0,\mathbb{V}\right)$

where

$$\mathbb{V}\ge V_{\text{eff}}^{\text{PATE}}:=\mathbb{E}\left[\frac{{\sigma }_1^2(\mathbf{X})}{\pi (\mathbf{X})}+\frac{{\sigma }_0^2(\mathbf{X})}{1-\pi (\mathbf{X})}+{\left(\tau (\mathbf{X})-\tau \right)}^2\right].$$

and for PATE ($\gamma$)

$$V_{\text{eff}}^{\text{PATT}}=\mathbb{E}\left[\frac{\pi (\mathbf{X}){\sigma }_1^2(\mathbf{X})}{p^2}+\frac{\pi (\mathbf{X}{)}^2{\sigma }_0^2(\mathbf{X})}{p^2\left(1-\pi (\mathbf{X})\right)}+\frac{{\left(\tau (\mathbf{X})-\gamma \right)}^2\pi (\mathbf{X})}{p^2}\right].$$

where ${\sigma }_d^2(\mathbf{X})=\mathbb{V}\left[Y^d|\mathbf{X}\right]$, $\tau (X):=\mathbb{E}\left[Y^1-Y^0|\mathbf{X}\right]$, and $p:=\mathbb{E}\left[\pi (\mathbf{X})\right]$.

Any regular estimator whose asymptotic variance achieves this efficiency bound is equal to $\frac{1}{n}{\sum }_{i=1}^n{\psi }_i(\mu )+O_P(\sqrt{n})$, where

$${\psi }_i(\mu )=\mu (1,{\mathbf{X}}_i)-\mu (0,{\mathbf{X}}_i)+\frac{D_i\left(Y_i-\mu (1,{\mathbf{X}}_i)\right)}{\pi ({\mathbf{X}}_i)}-\frac{(1-D_i)\left(Y_i-\mu (0,{\mathbf{X}}_i)\right)}{1-\pi ({\mathbf{X}}_i)}.$$

is the Efficient Influence Function for estimating $\tau$.

shows that

$$V_{\text{eff}}^{\text{SATE}}=V_{\text{eff}}^{\text{PATE}}-\underset{\text{Variance of treatment effect}}{\underbrace{\mathbb{E}\left[{\left[(Y_1-Y_0)-{\tau }^P\right]}^2\right]}}.$$

Estimators in this section try to attain the SPEB.

$${\tau }_{\text{ATE}}^{\text{IPW}}=\mathbb{E}\left[Y\cdot \frac{D-\pi (\mathbf{X})}{\pi (\mathbf{X})(1-\pi (\mathbf{X}))}\right]=\mathbb{E}\left[\frac{YD}{\pi (\mathbf{X})}-\frac{Y(1-D)}{1-\pi (\mathbf{X})}\right].$$

and

$${\tau }_{\text{ATT}}^{\text{IPW}}=\frac{1}{\mathbf{Pr}\left(D=1\right)}\mathbb{E}\left[Y\cdot \frac{D-\pi (\mathbf{X})}{1-\pi (\mathbf{X})}\right].$$

The counterfactual mean $\mathbb{E}\left[Y^0|D=1\right]={\mu }_0^1$ can be identified as

$\mathbb{E}\left[\frac{\pi ({\mathbf{X}}_i)}{\rho }\frac{1-D_i}{1-\pi ({\mathbf{X}}_i)}{Y}_i\right]$

where $\rho =\mathbf{Pr}\left(D=1\right)$.

$${\hat{\tau }}_{\text{ipw}}^{\text{ate}}=\frac{1}{n}\sum _{i=1}^n\left[\frac{Y_i{D}_i}{\hat{\pi }({\mathbf{X}}_i)}-\frac{Y_i(1-D_i)}{1-\hat{\pi }({\mathbf{X}}_i)}\right]=\frac{1}{n}\sum _{i=1}^n{Y}_i\left[\frac{D_i}{\hat{\pi }({\mathbf{X}}_i)}-\frac{1-D_i}{1-\hat{\pi }({\mathbf{X}}_i)}\right].$$

normalise both pieces using a Hajek-style adjustment, since extreme values of $\hat{\pi }$ makes variance explode. Often advisable to trim or use Hajek weights, which introduces limited bias at the cost of large decreases in variance.

$${\hat{\tau }}_{\text{ipw2}}^{\text{ate}}=\frac{{\sum }_{i=1}^n\frac{Y_i{D}_i}{\hat{\pi }({\mathbf{X}}_i)}}{{\sum }_{i=1}^n\frac{D_i}{\hat{\pi }({\mathbf{X}}_i)}}-\frac{{\sum }_{i=1}^n\frac{Y_i(1-D_i)}{1-\hat{\pi }({\mathbf{X}}_i)}}{{\sum }_{i=1}^n\frac{1-D_i}{1-\hat{\pi }({\mathbf{X}}_i)}}.$$

Similarly, for the effect on the treated

$${\hat{\tau }}_{\text{ipw}}^{\text{att}}=\frac{1}{n}\sum _{i=1}^n\left[\frac{Y_i}{\mathbf{Pr}\left(D_i=1\right)}\cdot \frac{D_i-\hat{\pi }(X_i)}{1-\hat{\pi }(X_i)}\right].$$ $${\hat{\tau }}_{\text{ipw2}}^{\text{att}}=\left[\frac{1}{N_1}\sum _{i:D_i=1}{Y}_i\right]-\left[\frac{{\sum }_{i:D_i=0}{Y}_i\cdot \frac{\hat{\pi }({\mathbf{X}}_i)}{1-\hat{\pi }({\mathbf{X}}_i)}}{{\sum }_{i:D_i=0}\frac{\hat{\pi }({\mathbf{X}}_i)}{1-\hat{\pi }({\mathbf{X}}_i)}}\right].$$

Horvitz-Thompson Estimator as Regression

$Y_i=\alpha +\tau D_i+{ϵ}_i$ with IPW weights ${\lambda }_i=\sqrt{\frac{D_i}{\pi (X_i)}+\frac{1-D_i}{1-\pi (X_i)}}$

define the Weighted ATE (WATE) as

$${\tau }_{\text{ATE}}=\frac{\int \mathbb{E}\left[Y^1-Y^0\mid \mathbf{X}=x\right]g(\mathbf{x})d\mathbb{F}(\mathbf{x})}{\int g(\mathbf{x})d\mathbb{F}(\mathbf{x})}.$$

where $g(\mathbf{x})$ is a weighting function. ATT is constructed when $g(\mathbf{x})=\pi (\mathbf{x})$

the corresponding estimator is

$${\hat{\tau }}_{\text{WATE}}=\frac{{\sum }_{i=1}^{n}g({\mathbf{x}}_i)\left[\frac{Y_i{D}_i}{\hat{\pi }({\mathbf{x}}_i)}-\frac{Y_i(1-D_i)}{1-\hat{\pi }({\mathbf{x}}_i)}\right]}{{\sum }_{i=1}^{n}g({\mathbf{x}}_i)}.$$

Sample drawn from $f(\mathbf{X})$, and can represent a target population as $g(\mathbf{X})\propto f(\mathbf{X})h(\mathbf{X})$ where $h(\cdot )$ is the tilting function.

Define $f_d(\mathbf{x})=\mathbf{Pr}\left(\mathbf{X}=\mathbf{x}|D=d\right)$, which gives $f_1(\mathbf{x})\propto f(\mathbf{x})\pi (\mathbf{x});f_0(\mathbf{x})\propto f(\mathbf{x})(1-\pi (\mathbf{x})$

For a given tilting function, to estimate ${\tau }_h$, weight $f_d(x)$

$$w_1(\mathbf{x})=\frac{h(\mathbf{x})}{\pi (\mathbf{x})},w_0(\mathbf{x})=\frac{h(\mathbf{x})}{1-\pi (\mathbf{x})}.$$

Target $h(x)$ Estimand $w_1,w_0$

---

Combined 1 ATE $\left(\frac{1}{\pi (x)},\frac{1}{1-\pi (x)}\right)$ [IPW] Treated $\pi (x)$ ATT $\left(1,\frac{\pi (x)}{1-\pi (x)}\right)$ Control $1-\pi (x)$ ATC $\left(\frac{1-\pi (x)}{\pi (x)},1\right)$ Overlap $\pi (x)(1-\pi (x))$ ATO $(1-\pi (x),\pi (x))$

Overlap weights are defined by choosing $h(x)$ that minimises asymptotic variance of ${\hat{\tau }}^h$. The achieve exact balance on covariates included in the propensity score estimation.

$$\hat{\tau }={\hat{\mu }}_1^h-{\hat{\mu }}_0^h=\frac{{\sum }_{i=1}^N{w}_1({\mathbf{x}}_i)D_i{Y}_i}{{\sum }_{i=1}^N{w}_1({\mathbf{x}}_i)D_i}-\frac{{\sum }_{i=1}^N{w}_0({\mathbf{x}}_i)(1-D_i)Y_i}{{\sum }_{i=1}^N{w}_0({\mathbf{x}}_i)(1-D_i)}.$$

${\tau }^{OW}$ can be interpreted as treatment effect among population that have good balance on observables.

Implemented in PSweight.

Entropy weights $w_i$ for each control unit are chosen by a reweighting scheme

${\max }_{w_i}H(w)=-{\sum }_{i:D=0}{w}_i\log (w_i)$

subject to balance/moment-condition and normalising constraints

$$\sum _{i:D_i=0}{w}_i{c}_{ri}({\mathbf{X}}_i)=m_r,r\in \{1,\dots ,R\}.$$ $$\sum _{i:D_i=0}{w}_i=1,w_i\ge 0\forall i\text{ such that }{D}_i=0.$$

The above problem is convex but has dimensionality of $n_0$ (nonnegativity) + $p$ (moment conditions) + $1$ (normalisation). The dual, on the other hand, only has dimensionality $p+1$ and unconstrained, which is considerably easier to solve using Newton-Raphson.

propose CBPS, which is a method that involves modifying an initial propensity score estimate (e.g. by changing coefficients from a logistic model) iteratively until a balance criterion is reached.

Their basic insight is that when we use a logistic regression to estimate a propensity score, we assert that the pscore takes the form ${\pi }_{\beta }({\mathbf{x}}_i)=\Lambda ({\mathbf{x}}_i^{\top }\mathbf{\beta })=\frac{\exp ({\mathbf{x}}_i^{\top }\mathbf{\beta })}{1+\exp ({\mathbf{x}}_i^{\top }\mathbf{\beta })}$, and maximise the bernoulli log likelihood

$$\sum _{i=1}^n\left[d_i\log \left({\pi }_{\beta }({\mathbf{x}}_i)\right)+(1-d_i)\log \left(1-{\pi }_{\beta }({\mathbf{x}}_i)\right)\right].$$

which is then solved by the corresponding score

$$\frac{1}{n}\sum _{i=1}^n\left[\frac{d_i{\pi }_{\beta }^{\prime }({\mathbf{x}}_i)}{{\pi }_{\beta }({\mathbf{x}}_i)}-\frac{(1-d_i){\pi }_{\beta }^{\prime }({\mathbf{x}}_i)}{1-{\pi }_{\beta }({\mathbf{x}}_i)}\right]=0.$$

this score balances a particular function of covariates: ${\pi }_{\beta }^{\prime }({\mathbf{x}}_i)$. Alternatively, we could choose that function by specifying a moment condition

$$\mathbb{E}\left[\frac{d_{i}f({\mathbf{x}}_i)}{{\pi }_{\beta }({\mathbf{x}}_i)}-\frac{(1-d_i)f({\mathbf{x}}_i)}{1-{\pi }_{\beta }({\mathbf{x}}_i)}\right]=0.$$

Analogously for ATT, this moment condition is

$$\mathbb{E}\left[d_{i}f({\mathbf{x}}_i)-\frac{{\pi }_{\beta }({\mathbf{x}}_i)}{1-{\pi }_{\beta }({\mathbf{x}}_i)}(1-d_i)f({\mathbf{x}}_i)\right]=0.$$

When this balance condition is solved independently, the problem is just-identified. When it is used in conjunction with the conventional bernoulli likelihood, the problem is over-identified. Implemented in CBPS::CBPS as well as balance.

The estimand is ${\mu }_1=\mathbb{E}\left[Y(1)\right]$ (with ${\mu }_0$ defined analogously). The estimator for this quantity is written

${\hat{\mu }}_1=\frac{1}{n}{\sum }_{i=1}^n{D}_i\hat{\gamma }({\mathbf{X}}_i)Y_i$

where the weights $\gamma (\cdot )$ are chosen to satisfy the sample balance property

$$\frac{1}{n}\sum _{i=1}^n{D}_i\gamma ({\mathbf{X}}_i)f({\mathbf{X}}_i)\approx \frac{1}{n}\sum _{i=1}^{n}f({\mathbf{X}}_i).$$

for any bounded function $f$.

in words: for every function $f(x)$, the weighting function equates weighted averages of $f$ over the treated units to unweighted averages over the study population.

The weights are solved by solving an optimisation problem to trade off imbalance and some measure of complexity

$$\hat{\gamma }=\underset{\gamma }{\mathrm{argmin}}\left\{\underset{\zeta (\cdot )}{\underbrace{{\text{imbalance}}_{\mathcal{M}}^2(\gamma )}}+\underset{\chi (\gamma \{X_i\})}{\underbrace{\frac{{\sigma }^2}{n^2}\sum _{D_i}\gamma ({\mathbf{X}}_i{)}^2}}\right\}.$$

with convex $\zeta ,\chi$ functions.

A common imbalance measure is

$${\text{imbalance}}_{\mathcal{M}}^2(\gamma )=\underset{j=1,\dots ,p}{\max }\left|\frac{1}{n}\sum _{i=1}^n{X}_{ij}-\frac{1}{n}\sum _{i=1}^n{D}_i\gamma ({\mathbf{X}}_i)X_{ij}\right|.$$

for $\mathcal{M}=\left\{\beta \cdot x:{\Vert \beta \Vert }_1\le 1\right\}$

Hybrid Estimators

A doubly-robust estimator is consistent if one gets either the propensity score $\hat{\pi }$ or the regression $\hat{\mu }$ right.

Oracle AIPW

$${\hat{\tau }}_{\text{AIPW}}^{\ast }=\frac{1}{n}\sum _{i=1}^n\left[{\mu }_{(1)}({\mathbf{X}}_i)-{\mu }_{(0)}({\mathbf{X}}_i)+\frac{D_i\left(y_i-{\mu }_{(1)}({\mathbf{X}}_i)\right)}{e({\mathbf{X}}_i)}+\frac{(1-D_i)\left(y_i-{\mu }_{(0)}({\mathbf{X}}_i)\right)}{1-e({\mathbf{X}}_i)}\right].$$

Feasible AIPW

$${\hat{\tau }}_{\text{AIPW}}=\frac{1}{N}\sum _{i=1}^n\left[\frac{D_i\left(Y_i-{\hat{\mu }}_1(X_i)\right)}{\hat{\pi }(X_i)}-\frac{(1-D_i)\left(Y_i-{\hat{\mu }}_0(X_i)\right)}{1-\hat{\pi }(X_i)}+\left({\hat{\mu }}_1(X_i)-{\hat{\mu }}_0(X_i)\right)\right].$$

This is the Augmented-Inverse-Propensity Weighting Estimator (AIPW) introduced by . Additional overviews:. General double-robustness property also shared by targeted maximum-likelihood estimators(TMLE) - due to .

Similarly, analogous estimator for ATT

$${\hat{\tau }}_{\text{aipw}}^{\text{att}}=\frac{1}{n}\sum _{i=1}^n\left[\frac{D_i\left(Y_i-{\hat{\mu }}_0({\mathbf{X}}_i)\right)}{\hat{\rho }}-\frac{\hat{\pi }({\mathbf{X}}_i)(1-D_i)\left(Y_i-{\hat{\mu }}_0({\mathbf{X}}_i)\right)}{\hat{\rho }\left(1-\hat{\pi }({\mathbf{X}}_i)\right)}\right].$$

where $\rho =\mathbf{Pr}\left(D_i=1\right)$ and $\hat{\rho }$ is its empirical analogue.

The Cross-fit version can be stated as

$${\hat{\tau }}_{\text{IPW}}=\frac{1}{n}\sum _{i=1}^n\left[{\hat{\mu }}_{(1)}^{-k(i)}({\mathbf{X}}_i)-{\hat{\mu }}_{(0)}^{-k(i)}({\mathbf{X}}_i)+D_i\frac{y_i-{\hat{\mu }}_{(1)}^{-k(i)}({\mathbf{X}}_i)}{{\hat{\pi }}^{-k(i)}({\mathbf{X}}_i)}-(1-D_i)\frac{y_i-{\hat{\mu }}_{(0)}^{-k(i)}({\mathbf{X}}_i)}{1-{\hat{\pi }}^{-k(i)}({\mathbf{X}}_i)}\right].$$

where $k(i)$ is a mapping that takes an observation and puts it into one of the $k$ folds. ${\hat{\mu }}_{(1)}^{-k(i)}$ is an estimator excluding the $k^{th}$ fold.

Define individual treatment effect score as

$${\hat{\Gamma }}_i={\hat{\mu }}_{(1)}({\mathbf{X}}_i)-{\hat{\mu }}_{(0)}({\mathbf{X}}_i)+\frac{D_i\left(Y_i-{\hat{\mu }}_{(1)}({\mathbf{X}}_i)\right)}{\hat{\pi }({\mathbf{X}}_i)}-\frac{(1-D_i)\left(Y_i-{\hat{\mu }}_{(0)}({\mathbf{X}}_i)\right)}{1-\hat{\pi }({\mathbf{X}}_i)}.$$

Then, $\hat{\tau }=\frac{1}{n}{\sum }_{i=1}^n{\hat{\Gamma }}_i$

We can form level-$\alpha$ CIs ${\mathcal{I}}_{\alpha }$ :

$${\mathcal{I}}_{\alpha }=\hat{\tau }\pm z_{1-\alpha /2}{\hat{V}}^{1/2},\hat{V}=\frac{1}{n(n-1)}\sum _{i=1}^n{\left({\hat{\Gamma }}_i-\hat{\tau }\right)}^2.$$

grf has a forest-based implementation of AIPW

cf = causal_forest(X, Y, D)
ate_hat = average_treatment_effect(cf)

partially-linear setup

$$y_i=d_i\tau +g({\mathbf{x}}_i)+{\epsilon }_i,\mathbb{E}\left[{\epsilon }_i\mid {\mathbf{x}}_i,d_i\right]=0.$$ $$d_i=m({\mathbf{x}}_i)+{\eta }_i,\mathbb{E}\left[{\eta }_i\mid {\mathbf{x}}_i\right]=0.$$

where $d_i$ is a scalar treatment indicator. Observations are independent but not necessarily identically distributed. We are interested in inference about $\tau$ that is robust to mistakes in model-selection.

Approximate $g$ and $m$ with linear combinations of control terms ${\mathbf{c}}_i=P({\mathbf{x}}_i)$, which may contain interactions and non-linear transformations.

Assume approximate sparsity ($:=$ there are only a small number of relevant controls, and irrelevant controls have a high probability of being small).

Naive (incorrect) approach: use LASSO on an eqn of the form

$$y_i=\tau D_i+{\mathbf{x}}_i^{\prime }\beta +{\epsilon }_i,h(\mathbf{\beta })={\Vert \mathbf{\beta }\Vert }_1=\sum _j\left|{\beta }_j\right|.$$

where the treatment $\tau$ is not penalised. This will mean we drop any control that is highly correlated with the treatment if the control is moderately correlated with the outcome. Then, if we use a post-LASSO selection to estimate the treatment effect, the effect will be contaminated with an omitted variable bias.

recommended two-step approach

Estimate $y_i={\mathbf{c}}_i^{\prime }\mathbf{\beta }+{\nu }_i$ with LASSO, select predictive variables (i.e. those with nonzero coefficients) in $\mathcal{A}$

Estimate $d_i={\mathbf{c}}_i^{\prime }\mathbf{\beta }+{\nu }_i$ with LASSO, select predictive variables (i.e. those with nonzero coefficients) in $\mathcal{B}$

Estimate $y_i=\tau D_i+{\mathbf{e}}_i^{\prime }\mathbf{\kappa }+{\upsilon }_i$ where ${\mathbf{c}}_i:=\mathcal{A}\cup \mathcal{B}$ [i.e. control for variables that are selected in either the first or second regression]

Let ${\hat{l}}_1,{\hat{l}}_2\subset \left\{1,\dots ,p\right\}$ be the indices of the selected controls for the outcome and treatment respectively.

The post-double-selection estimator is

$$(\check{\tau },\check{\mathbf{\beta }})=\underset{\tau \in \mathbb{R},\mathbf{\beta }\in {\mathbb{R}}^k}{\mathrm{argmin}}{\mathbb{E}}_n\left[(y_i-d_i\tau -{\mathbf{x}}_i^{\prime }\mathbf{\beta }{)}^2\right]\text{s.t.}{\beta }_j=0\forall j\notin {\hat{l}}_1\cup {\hat{l}}_2.$$

Can use plugin estimator for variance based on residuals

$${\sigma }_n^{-1}\sqrt{\check{\tau }-\tau }\stackrel{d}{\to }\mathcal{N}(0,1),{\sigma }_n^2=\frac{\mathbb{E}\left[v_i^2{\psi }_i^2\right]}{\mathbb{E}{\left[v_i^2\right]}^2}.$$

where $\begin{array}{rl}{\hat{\psi }}_i& =(y_i-d_i\check{\tau }-{\mathbf{x}}_i^{\prime }\check{\mathbf{\beta }})\sqrt{\frac{n}{n-\hat{s}-1}}\\ {\hat{v}}_i& =d_i-{\mathbf{x}}_i^{\prime }\hat{\mathbf{\beta }}\\ \hat{\beta }& =\underset{\mathbf{\beta }\in {\mathbb{R}}^p}{\mathrm{argmin}}\left\{{\mathbb{E}}_n\left[(d_i-{\mathbf{x}}_i^{\prime }\mathbf{\beta }{)}^2\right]:{\beta }_j=0\forall j\notin \hat{l}\right\}\end{array}$

Implemented in hdm::rlassoEffect(., ’double selection’)

Let the target parameter ${\tau }_0$ solve the equation $\mathbb{E}\left[m(Z_i,{\tau }_0,{\beta }_0)\right]=0$ for known score function $m$, vector of observables $Z_i:={\left\{{\mathbf{X}}_i,D_i,Y_i\right\}}_{i=1}^n$, and nuisance parameter ${\beta }_0$. In fully parametric models, $m$ is simply the score function [derivative of the log-likelihood]. For ATE, $m(Z_i,\tau ,\beta ):=(Y_i-D_i\tau -{\mathbf{X}}_i^{\prime }\beta )D_i$.

In naive double-ML settings, $\mathbb{E}\left[{\partial }_{\beta }m(Z_i,{\tau }_0,{\beta }_0)={\pi }_0{\mu }_1\ne 0\right]$. So, we replace $m$ with the Neyman-orthogonal score $\psi$ s.t.

$\mathbb{E}\left[{\partial }_{\eta }\psi (Z_i,{\tau }_0,{\eta }_0)=0\right].$

which yields the Orthogonalised Moment Condition $\mathbb{E}\left[\psi (Z_i,{\tau }_0,{\eta }_0)\right]=0$ for some real-valued condition $\psi (.)$.

Using a Neyman-orthogonal score eliminates first-order biases arising from the replacement of ${\eta }_0$ with ${\hat{\eta }}_0$.

Reference:

Consider data $\mathcal{W}:=(Y,D,X)$ with $D\in \left\{0,1\right\}$

Partial linear setup $Y=D{\theta }_0+g_0(X)+U;D=m_0(X)+V$.

Score function is

$$\psi(\Sett{W}; \theta, \eta) = (Y - \Ubr{l(X)}{$\mathbb{E}\left[Y|X\right]$} - \theta(D

• \Ubr{m(X)}{$\mathbb{E}\left[D|X\right]$})) (D - m(X))$$

Partially Linear IV

$$Y-D{\theta }_0=g_0(X)+\zeta ,\mathbb{E}(\zeta \mid Z,X)=0.$$ $$Z=m_0(X)+V,\mathbb{E}(V\mid X)=0.$$

Score is

$\psi (\mathcal{W},\theta ,\eta ):=(Y-\underset{\text{E[Y∣X]}}{\underbrace{l(X)}}-\theta (D-\underset{\text{E[D∣X]}}{\underbrace{r(X)}}))(Z-\underset{\text{E[Z∣X]}}{\underbrace{m(X)}})$

Interactive Regression $\begin{array}{rl}Y& =\underset{\text{E[Y∣D,X]}}{\underbrace{g_0(D,X)}}+\epsilon ,\mathbb{E}\left[\epsilon |D,X\right]=0\\ D& =\underset{\text{PrD=1∣X}}{\underbrace{m_0(X)}}+\xi ,\mathbb{E}\left[\xi |X\right]=0\end{array}$

Here, the estimands are

$${\theta }_0^{ATE}=\mathbb{E}\left[g_0(1,X)-g_0(0,X)\right].$$ $${\theta }_0^{ATT}=\mathbb{E}\left[g_0(1,X)-g_0(0,X)\mid D=1\right].$$

The score function for ATE (Hahn (1998))

$${\psi }^{ATE}(Z_i,\theta ,\eta )=\left(g(1,X)-g(0,X)\right)+\frac{D\left(Y-g(1,X)\right)}{m(X)}-\frac{(1-D)\left(Y-g(0,X)\right)}{1-m(X)}-\theta .$$

The nuisance parameter true value is ${\eta }_0=(g_0,m_0)$. For ATET,

$${\psi }^{ATT}(Z_i,\theta ,\eta )=\frac{1}{{\pi }_0}\left(D-\frac{m(X)}{1-m(X)}(1-D)\right)\left(Y-g(0,X)\right)-\frac{D}{{\pi }_0}\theta .$$

Take a $K$-fold random partition $(I_k{)}_{k=1,\dots ,K}$ of observation indices $\left\{1,\dots ,n\right\}$ s.t. each fold $I_k$ has size $n/K$. For each $k$, define $I_k^C:=\left\{1,\dots ,n\right\}{I}_k$ as the complement / auxiliary sample.

For each $k\in \left\{1,\dots ,K\right\}$, construct a ML estimator of ${\eta }_0$ using only the auxiliary sample $I_k^C$; ${\hat{\eta }}_k=\hat{\eta }((Z_i{)}_{i\in I_k^C})$

For each $k\in \left\{1,\dots ,K\right\}$, using the main sample $I_k$, construct the estimator ${\check{\tau }}_k$ as the solution of $\frac{1}{n/K}{\sum }_{i\in I_k}\psi (Z_i,{\check{\tau }}_k,{\hat{\eta }}_k)=0$

Aggregate the estimators ${\check{\tau }}_k$ on each main sample $\check{\tau }=\frac{1}{K}{\sum }_{k=1}^K{\check{\tau }}_k$

Simple implementation of Cross-fitting for Treatment effects

Partition the data in two, such that each fold $I_1,I_2$ has size $n/2$.

Using only sample $I_1$, construct a ML estimator of $g(0,X)$ and $m(X)$,e.g. a feedforward nnet of $Y_i$ on $X_i$, denoted as $\widehat{g_{I_1}}(x)$, and logit-lasso of $D_i$ on $X_i$, denoted by $\widehat{m_{I_1}}(x)$.

Use the estimators on the hold-out sample $I_2$ to compute the T.E ${\check{\tau }}_{I_2}=\frac{1}{{\sum }_{i\in I_2}{D}_i}\left[D_i-\frac{{\hat{m}}_{I_1}(X_i)}{1-{\hat{m}}_{I_1}(X_i)}(1-D_i)\right](Y_i-{\hat{g}}_{I_1}(X_i))$

Repeat (2,3) swapping the roles of $I_1$ and $I_2$ to get ${\check{\tau }}_{I_1}$

Aggregate the estimators: $\check{\tau }=\frac{{\check{\tau }}_{I_1}+{\check{\tau }}_{I_2}}{2}$

Implemented in DoubleML

Augmented Balancing

Loosely: AIPW without the (potentially fraught) inversion of the propensity score step. Exposition based on Bruns-Smith et al (2023)

setup: Covariates $X\in \mathcal{X}\subseteq {\mathbb{R}}^k$ , $Y\in \mathbb{R}$ outcome, two populations $p$ and $q$ that are distributions over $(X,Y)$ $p$ is ‘source’, $q$ is ‘target’ (e.g. treatment group and overall sample) Estimand is ${\mathbb{E}}_q\left[Y\right]$ Identification Assumptions Conditional Mean Ignorability: ${\mathbb{E}}_p\left[Y|X\right]={\mathbb{E}}_q\left[Y\mid X\right]$ Population Overlap: $q(x)$ is absolutely continuous w.r.t. $p(x)$

Effect Functionals

Regression Functional

$${\mathbb{E}}_q\left[{\mathbb{E}}_p\left[Y\mid X\right]\right]={\mathbb{E}}_q\left[{\mathbb{E}}_q\left[Y\mid X\right]\right]={\mathbb{E}}_q\left[Y\right].$$

Weighting Functional

$${\mathbb{E}}_p\left[\frac{dq}{dp}(X)Y\right]={\mathbb{E}}_p\left[\frac{dq}{dp}(X){\mathbb{E}}_p\left[Y\mid X\right]\right]={\mathbb{E}}_q\left[{\mathbb{E}}_p\left[Y\mid X\right]\right]={\mathbb{E}}_q\left[Y\right].$$

Doubly-Robust Functional

$${\mathbb{E}}_q\left[{\mathbb{E}}_p\left[Y\mid X\right]\right]+{\mathbb{E}}_p\left[\frac{dq}{dp}(X)\left(Y-{\mathbb{E}}_p\left[Y\mid X\right]\right)\right].$$

Balancing Weights: Rationale

$\frac{dq}{dp}(X)$ is difficult to estimate using plug-in estimation Alternative: weighting for balance $\equiv$ automatic estimation of the Riesz representer

Weighting to minimise covariate imbalance

$$\underset{w}{\min }\left\{\underset{\text{Imbalance}}{\underbrace{\underset{f\in \mathcal{F}}{\sup }\left({\mathbb{E}}_p\left[w(X)f(X)\right]-{\mathbb{E}}_q\left[f(X)\right]\right)}}+\delta {\Vert w\Vert }_2^2\right\}.$$

Direct estimation of the density ratio

$$\underset{f\in \mathcal{F}}{\min }{\mathbb{E}}_p\left[{\left(f(X)-\frac{dq}{dp}(X)\right)}^2\right].$$

Minimum variance weights that balance $\mathcal{F}$ are also guaranteed to balance all other measurable functions in $\mathcal{F}$.

In linear setting, relevant imbalance is captured entirely by feature mean imbalance $X_p,Y_p$ are $n$ iid draws from $p$, $X_q$ are m draws from $q$ Define feature map $\varphi :\mathcal{X}\to {\mathbb{R}}^d$; construct gram matrices ${\Phi }_p:=\varphi (X_p)$ ${\overline{\Phi }}_q:={\hat{\mathbb{E}}}_q[{\Phi }_q]$ Let $\mathcal{F}=\left\{f(x)={\theta }^{\top }\varphi (x):\Vert \theta \Vert \le r\right\}$ Let ${\Vert \cdot \Vert }_{\ast }$ denote dual norm [${\ell }_2:{\ell }_2$, ${\ell }_1:{\ell }_{\infty }$]

${\widehat{\text{Imbalance}}}_{\mathcal{F}(w)}={\Vert w{\Phi }_p-{\overline{\Phi }}_q\Vert }_{\ast }$

Three Equivalent Representations

Penalised form

$$\underset{w\in {\mathbb{R}}^n}{\min }\left\{{\Vert w{\Phi }_p-{\overline{\Phi }}_q\Vert }_{\ast }^2+{\delta }_1{\Vert w\Vert }_2^2\right\}.$$

Constrained form

$$\underset{w\in {\mathbb{R}}^n}{\min }{\Vert w\Vert }_2^2\text{s.t.}\Vert w{\Phi }_p-{\overline{\Phi }}_q\Vert \le {\delta }_2.$$

Automatic form

$$\underset{\theta \in {\mathbb{R}}^d}{\min }\left\{{\theta }^{\top }({\Phi }_p^{\top }{\Phi }_p)\theta -2{\theta }^{\top }{\overline{\Phi }}_q+{\delta }_3\Vert \theta \Vert \right\}.$$

OLS is equivalent to a weighting estimator that exactly balances the feature means. Let ${\hat{\beta }}_{\text{OLS}}=({\Phi }_p^{\top }{\Phi }_p{)}^{-1}{\Phi }_p^{\top }{Y}_p$ be the linear regression fit on $p$ (source sample). Then,

$${\hat{\mathbb{E}}}_q[{\Phi }_q{\hat{\beta }}_{\text{OLS}}]={\hat{\mathbb{E}}}_p[{\hat{w}}_{\text{exact}}\circ Y_p].$$ $${\hat{\mathbb{E}}}_q[{\Phi }_q({\Phi }_p^{\top }{\Phi }_p{)}^{-1}{\Phi }_p^{\top }{Y}_p]={\hat{\mathbb{E}}}_p[{\overline{\Phi }}_q({\Phi }_p^{\top }{\Phi }_p{)}^{-1}{\Phi }_p^{\top }{Y}_p].$$

Analogue for Ridge

$${\hat{\mathbb{E}}}_q[{\Phi }_q({\Phi }_p^{\top }{\Phi }_p+\delta I{)}^{-1}{\Phi }_p^{\top }{Y}_p]={\hat{\mathbb{E}}}_p[{\overline{\Phi }}_q({\Phi }_p^{\top }{\Phi }_p+\delta I{)}^{-1}{\Phi }_p^{\top }{Y}_p].$$

$\forall {\hat{\beta }}_{reg}^{\lambda }\in {\mathbb{R}}^d$, and any linear balancing weight estimator with estimated coefficients ${\theta }^d\in {\mathbb{R}}^d$, ${\hat{w}}^{\delta }=\hat{\theta }{\Phi }_p^{\top }$, and ${\hat{\Phi }}_q^{\delta }=\hat{w}{\Phi }_p$

$${\hat{\mathbb{E}}}_q[{\Phi }_q{\hat{\beta }}_{reg}^{\lambda }]+{\hat{\mathbb{E}}}_p[{\hat{w}}^{\delta }\circ (Y_p-{\Phi }_p{\hat{\beta }}_{reg}^{\lambda })]={\hat{\mathbb{E}}}_p[{\hat{w}}^{\delta }\circ Y_p]+{\hat{\mathbb{E}}}_q\left[({\overline{\Phi }}_q-{\hat{\Phi }}_q^{\delta }){\hat{\beta }}_{reg}^{\lambda }\right]={\hat{\mathbb{E}}}_q\left[{\hat{\Phi }}_q^{\delta }{\hat{\beta }}_{OLS}+({\overline{\Phi }}_q-{\hat{\Phi }}_q^{\delta }){\hat{\beta }}_{reg}^{\lambda }\right]={\hat{\mathbb{E}}}_q[{\Phi }_q{\hat{\beta }}_{aug}].$$ $${\hat{\beta }}_{aug,j}:=(1-a_j^{\delta }){\hat{\beta }}_{reg,j}^{\lambda }+a_j^{\delta }{\hat{\beta }}_{OLS,j},a_j^{\delta }:=\frac{{\hat{\Phi }}_{q,j}^{\delta }-{\overline{\Phi }}_{p,j}}{{\overline{\Phi }}_{q,j}-{\overline{\Phi }}_{p,j}}.$$

In words: when both outcome and weighting models are linear, the augmented estimator is equivalent to a linear model with coefficients that are element-wise affine combinations of base learner ${\hat{\beta }}_{reg}^{\lambda }$ and coefs ${\hat{\beta }}_{OLS}$ from regressing $Y_p$ on ${\Phi }_p$

Heterogeneous Treatment Effects with selection on observables

Conditional Average Treatment Effects (CATEs) ($\tau (x)=\mathbb{E}\left[Y^1-Y^0|\mathbf{X}=\mathbf{x}\right]$) are often of great policy interest for targeting those who have largest potential gains. However, conventional methods are prone to a severe risk of fishing from researchers (cf ‘conditional effects’ in most published work in the social sciences).

Instead, recent work proposes to use nonparametric estimators to find subgroups, use sample-splitting for honesty.

transformed outcome regression use outcome transformed w pscore $H=\frac{DY}{p(\mathbf{X})}-\frac{(1-D)Y}{1-p(\mathbf{X})}$

conditional mean regression use the fact that under SOO $\tau (\mathbf{x})=\mathbb{E}\left[Y_1|\mathbf{X}=\mathbf{x}\right]-\mathbb{E}\left[Y_0|\mathbf{X}=\mathbf{x}\right]={\mu }_1(\mathbf{x})-{\mu }_0(\mathbf{x})$

(1) typically inefficient because of pscore in denominator, so most focus is on (2). Random forests are a flexible method that is widely liked.

Consider a model for $\tau (x)$ where

$Y_i(d)=f({\mathbf{X}}_i)+d\cdot \tau ({\mathbf{X}}_i)+\epsilon (d),\mathbf{Pr}{d}_i|{\mathbf{X}}_i=e(x)$

where $\tau (\mathbf{x})=\psi (\mathbf{x})\mathbf{\beta }$ for some pre-determined set of basis functions $\psi :\mathcal{X}\to {\mathbb{R}}^k$. We allow for non-parametric relationships between ${\mathbf{X}}_i,y_i,d_i$, but the treatment effect function itself is parametrised by $\beta \in {\mathbb{R}}^k$. showed that under unconfoundedness, we can rewrite the semiparametric setup above as

$$Y_i-m({\mathbf{X}}_i)=(d_i-e({\mathbf{X}}_i))\psi ({\mathbf{X}}_i)\cdot \beta +{\epsilon }_i.$$

where

$$m(\mathbf{x})=\mathbb{E}\left[Y_i\mid {\mathbf{X}}_i=\mathbf{x}\right]=f({\mathbf{X}}_i)+e({\mathbf{X}}_i)\tau ({\mathbf{X}}_i).$$

The oracle algorithm for estimating $\mathbf{\beta }$ is (1) define ${\tilde{Y}}_i^{\ast }=Y_i-m({\mathbf{X}}_i)$ and ${\tilde{Z}}_i^{\ast }=(d_i-e({\mathbf{X}}_i)\psi ({\mathbf{X}}_i)$, then estimate residuals-on-residual regression. This procedure is $\sqrt{n}$-consistent and asymptotically normal.

Use cross-fitting to emulate the Oracle.

Run non-parametric regressions $Y\sim \mathbf{X}$ and $D\sim \mathbf{X}$ to get $\hat{m}(\mathbf{x}),\hat{e}(\mathbf{x})$

define transformed features ${\tilde{Y}}_i=Y_i-{\hat{m}}^{-k(i)}({\mathbf{X}}_i)$, $\tilde{Z}=(D_i-{\hat{e}}^{-k(i)}({\mathbf{X}}_{\mathbf{i}})\psi ({\mathbf{X}}_{\mathbf{i}})$

Estimate ${\hat{\zeta }}_b$ by regressing ${\tilde{Y}}_i\sim {\tilde{Z}}_i$

To define R-Loss , under more general setup restate unconfoundedness as follows

$$\mathbb{E}\left[{\epsilon }_i(d_i)\mid {\mathbf{X}}_i,d_i\right]=0.$$

where

$${\epsilon }_i(d):=Y_i(d)-\left({\mu }_{(0)}({\mathbf{X}}_i)+d\tau ({\mathbf{X}}_i)\right).$$

and follow Robinson’s approach to write

$Y_i-m({\mathbf{X}}_i)=(D_i-e({\mathbf{X}}_i))\tau ({\mathbf{X}}_i)+{\epsilon }_i$

R-loss is then written

$$\tau (\cdot )=\underset{{\tau }^{\prime }}{\mathrm{argmin}}\mathbb{E}\left[{\left(Y_i-m({\mathbf{X}}_i)-(d_i-e({\mathbf{X}}_i)){\tau }^{\prime }({\mathbf{X}}_i)\right)}^2\right].$$

Define $e(x)=\mathbf{Pr}\left(D=1|X=x\right)$ and $m(x)=\mathbb{E}\left[Y|X=x\right]$ The R-learner consists of the following steps

Use any method to estimate the response functions $\hat{e}(x),\hat{m}(x)$

Minimise R-loss using cross-fitting for nuisance components

$$\hat{\tau }(\cdot )=\underset{\tau }{\mathrm{argmin}}\left[{\left((Y_i-\hat{m}({\mathbf{X}}_i))-\tau ({\mathbf{X}}_i)(D_i-\hat{e}({\mathbf{X}}_i))\right)}^2+{\Lambda }_n(\tau (\cdot ))\right].$$

where ${\Lambda }_n$ is a regulariser.

Causal forest as implemented by grf starts by fitting two separate trees to estimate $\hat{m},\hat{e}$, makes out-of-bag predictions [using cross-fitting] using the two first-stage forests, then grow causal forest via

$$\hat{\tau }(x)=\frac{{\sum }_{i=1}^n{\alpha }_i(x)\left(Y_i-\hat{m}({\mathbf{X}}_i)\right)\left(D_i-\hat{e}({\mathbf{X}}_i)\right)}{{\sum }_{i=1}^n{\alpha }_i(x){\left(D_i-\hat{e}({\mathbf{X}}_i)\right)}^2}.$$

where ${\alpha }_i(x)=\frac{1}{B}{\sum }_b\frac{1_{{\mathbf{X}}_{\mathbf{i}}\in {\mathcal{L}}_b(x),i\in B}}{|i:X_i\in {\mathcal{L}}_b(x),i\in B|}$

are the learned adaptive weights.

Draw a subsample of size $s$ from the sample with replacement and divide it into disjoint sets $\mathcal{I},\mathcal{J};\left|\mathcal{I}\right|=\left|\mathcal{J}\right|=n/2$.

Grow a tree via recursive partitioning, with splits chosen from $\mathcal{J}$ (i.e. without using $Y$ observations from $\mathcal{I}$ sample)

Estimate leaf responses using only $\mathcal{I}$ sample

Finally, aggregate all trees over subsamples of size $s$

$$\hat{\mu }(\mathbf{x},{\mathbf{Z}}_1,\dots ,{\mathbf{Z}}_n)={\left(\genfrac{}{}{0ex}{}{n}{s}\right)}^{-1}\sum _{1\le i_1<\cdots <i_s\le n}{\mathbb{E}}_{\xi \in \Xi }\left[T(\mathbf{x},\xi ,{\mathbf{Z}}_{i_1},\dots ,{\mathbf{Z}}_{i_s})\right]\approx \frac{1}{B}\sum _{b=1}^{B}T(\mathbf{x},{\xi }_b^{\ast },{\mathbf{Z}}_{b,1}^{\ast },\dots ,{\mathbf{Z}}_{b,s}^{\ast }).$$

where $\xi$ summarises randomness in the selection of the variable when growing the tree, ${\mathbf{Z}}_i:=(D_i,{\mathbf{X}}_i,Y_i)$ is shorthand for a training sample.

where the base learner $T(\mathbf{x};{\xi }_b^{\ast },{\mathbf{Z}}_{b,1}^{\ast },\dots ,{\mathbf{Z}}_{b,s}^{\ast })={\sum }_{i\in \left\{i_{b,1},\dots i_{b,s}\right\}}{\alpha }_{i,b}^{\ast }(\mathbf{x})Y_{i,b}^{\ast };{\alpha }_{i,b}^{\ast }(\mathbf{x})=\frac{1_{{\mathbf{X}}_{i,b}^{\ast }\in {\mathcal{L}}_b^{\ast }(\mathbf{x})}}{|i:{\mathbf{X}}_{i,b}^{\ast }\in {\mathcal{L}}_b^{\ast }(\mathbf{x})|}$

the ‘honesty’ property is making ${\alpha }_{i,b}^{\ast }(\mathbf{x})$ independent of $Y_{i,b}^{\ast }$, i.e. do not use the same data to select partition (splits) and make predictions.

Implemented in causalForest and grf.

Multi-action policy learning

$i=1,\dots ,N$ units, to be assigned to $J+1$ actions $A_i\in \{0,1,\dots ,J\}=:\mathcal{A}$, which has have corresponding rewards $\{Y_i^{(0)},Y_i^{(1)},\dots ,Y_i^{(J)}\}$. Each observation has covariate ${\mathbf{X}}_i\in \mathcal{X}\subseteq {\mathbb{R}}^d$. Define a policy function

$\pi :\mathbf{x}\to \mathcal{A}$

A given policy assigns each unit to a treatment level. Each policy has a corresponding value function

$V(\pi ):=\mathbb{E}\left[Y(\pi (\mathbf{x}))\right]$

An optimal policy ${\pi }^{\ast }\in \Pi$ is defined as

${\pi }^{\ast }={\mathrm{argmax}}_{\pi \in \Pi }\mathbb{E}\left[Y(\pi )\right]$

Deviations from this optimum is called regret

$R(\pi )=\mathbb{E}\left[Y({\pi }^{\ast })\right]-\mathbb{E}\left[Y(\pi )\right]=V({\pi }^{\ast })-V(\pi )$

Define a CEF as

${\mu }_i(a,{\mathbf{x}}_i)=\mathbb{E}\left[Y_i^{(a)}\mid {\mathbf{x}}_i\right]$

The first-best optimal rule is

${\pi }_i({\mathbf{x}}_i)={\mathrm{argmax}}_{a\in \mathcal{J}}\left\{{\mu }_i(a,{\mathbf{x}}_i)\right\}$

In the binary action case, this simplifies to $\pi ({\mathbf{x}}_i)=1\left\{\mu (1,{\mathbf{x}}_i)\ge \mu (0,{\mathbf{x}}_i)\right\}=1\left\{\tau (x)>0\right\}$ which is the conditional empirical success (CES) rule of Manski (2004).

Under unconfoundedness and Overlap, we can estimate $\hat{\mu }$s and construct an empirical analogue of the value function for a policy $\pi$ using the following familiar estimators

$${\hat{V}}_{RA}(\pi )=\frac{1}{n}\sum _{i=1}^n{\hat{\mu }}_i(\pi ({\mathbf{x}}_i),{\mathbf{x}}_i).$$ $${\hat{V}}_{IPW}(\pi )=\frac{1}{n}\sum _{i=1}^n\frac{1[A_i=\pi ({\mathbf{x}}_i)]}{{\hat{p}}_{A_i}({\mathbf{X}}_i)}{Y}_i.$$ $${\hat{V}}_{AIPW}(\pi )=\frac{1}{n}\sum _{i=1}^n\left[{\hat{\mu }}_i(\pi ({\mathbf{x}}_i),{\mathbf{x}}_i)+\frac{1[A_i=\pi ({\mathbf{x}}_i)]}{{\hat{p}}_{A_i}({\mathbf{x}}_i)}\left(Y_i-{\hat{\mu }}_i(\pi ({\mathbf{x}}_i),{\mathbf{x}}_i)\right)\right].$$

A $\sqrt{n}$-convergent estimator of the value function is the Cross-fit Augmented Inverse Probability Weighted Learning (CAIPWL) estimator of , which is constructed as a cross-fit analogue of the AIPW estimator.

Sensitivity Analysis

Check balance by computing SDiff for observable confounders

$\text{Standardised Difference}=\frac{{\overline{X}}_t-{\overline{X}}_c}{\sqrt{(s_t^2+s_c^2)/2}}$

Three valued treatment indicator: $T_i\in \{-1,0,1\}$ corresponding with ineligibles, eligible nonparticipants, and participants. We can test unconfoundedness by comparing ineligibles with eligible nonparticipants, i.e. test

$Y_i\perp \perp 1\left\{T_i=0\right\}|X_i,T_i\in \{-1,0\}$

Placebo Outcomes

Covariates included lagged outcomes $Y_{i,-1},\dots Y_{i,-T}$. Test $Y_{i,-1}\perp \perp D_i|Y_{i,-2},Y_{i,-T},X_i$ e.g. Earnings in 1975 in Lalonde

$U$ is a nuisance parameter.

$Y_1,Y_0\perp \perp D|X,U$

Where $U\sim \mathcal{B}(\pi =0.5)$, and $U\perp \perp X$. $P(U=1)=P(U=0)=0.5$.

Propensity score is Logistic:

$P(D=1|X,U)=\frac{\exp (X\theta +\gamma U)}{1+\exp (X\theta +\gamma U)}$

$\gamma$ indicates strength of relationship between $U$ and $D|X$.

Y is conditionally normal $Y|X,U\sim \mathcal{N}(\alpha D+X\beta +\delta U,{\sigma }^2)$

$\delta$ indicates strength of relationship between $U$ and $Y|X$.

MLE setup

Construct grid of $(\gamma ,\delta )$ and calculate the MLE for $\hat{\alpha }(\gamma ,\delta )$ by maximising $l(\alpha ,\beta ,\theta ,\gamma ,\delta )$ over $(\gamma ,\delta )$.

Use 2 partial $R^2$s:

• $R_{Y,par}^2(\delta )$: Residual variation in outcome explained by $U$ (after partialling out $X$).

• $R_{D,par}^2(\gamma )$: Residual variation in treatment assignment explained by $U$ (after partialling out $X$).

Draw threshold contours, should expect most covariates to be clustered around origin.

Rosenbaum (2002)[]

Tuning parameter $\Gamma \ge 1$ that measures departure from zero hidden bias.

For any two observations $i$ and $j$ with identical covariate values $X_i=X_j$, under unconfoundedness, probability of assignment into treatment should be identical $\pi (X_i)=\pi (X_j)$

Treatment assignment probability may differ due to unobserved binary confounder $U$. We can bound this by the ratio:

$$\frac{1}{\Gamma }\le \frac{{\hat{\pi }}_i(1-{\hat{\pi }}_j)}{(1-{\hat{\pi }}_i){\hat{\pi }}_j}\le \Gamma .$$

$\gamma =1⟹$ No bias. $\Gamma =2⟹$ $i$ is twice as likely to be treated than $j$ despite identical $x$.

$\Gamma$ is assumed to satisfy

$$\frac{1}{\Gamma }\le \frac{\mathbf{Pr}\left(D=d\mid X=x\right)/(1-\mathbf{Pr}\left(D=d\mid X=x\right))}{\mathbf{Pr}\left(D=d\mid X=x,Y(d)=y\right)/(1-\mathbf{Pr}\left(D=d\mid X=x,Y(d)=y\right))}\le \Gamma .$$

For any given candidate $\Gamma >1$, estimates of the treatment effect can be computed. Implemented in rbounds::hlsens.

Altonji, Elder, Taber (2005)

Only informative if selection on observables is informative about selection on unobservables.

How much does treatment effect move when controls are added? Estimate model with and without controls:

• $Y_i={\alpha }^F{D}_i+X\beta +ϵ$

• $Y_i={\alpha }^R{D}_i+ϵ$

AET ratio: $\rho =\frac{{\hat{\alpha }}^F}{{\hat{\alpha }}^R-{\hat{\alpha }}^F}$

Want $\rho$ to be as big as possible (i.e. ${\hat{\alpha }}^R-{\hat{\alpha }}^F\to 0$ under unconfoundedness).

Define proportional selection coefficient

$$\delta =\frac{\mathrm{Cov}\left(ϵ,D\right)/\mathbb{V}\left[ϵ\right]}{\mathrm{Cov}\left(X^{\prime }\gamma \right)/\mathbb{V}\left[X^{\prime }\gamma \right]}.$$

Then,

$${\beta }^{\ast }\approx \tilde{\beta }-\delta (\dot{\beta }-\tilde{\beta })\frac{R_{max}-\tilde{R}}{\tilde{R}-\dot{R}}\stackrel{p}{\to }\beta .$$

where

$\dot{\beta },\dot{R}$ are from a univariate regression of $Y$ on $T$

$\tilde{\beta },\tilde{R}$ are from a regression including controls

$R_{max}$ is maximum achievable $R^2$

True model is $Y=\tau D+\mathbf{X}\mathbf{\beta }+\gamma Z+\epsilon$, but we don’t observe $Z$. We would like to quantify how biased the coefficient from the short regression ${\hat{\tau }}_s$ is for the long regression coefficient $\tau$. From OVB FOrmula, we know ${\hat{\tau }}_s=\hat{\tau }+\hat{\gamma }+\hat{\delta }$ where $\hat{\gamma }$ is the conditional association between the omitted $Z$ and $Y$ (‘impact’) and $\hat{\delta }$ is the coefficient from regressing $Z$ on $D$ (‘imbalance’).

The bias from this omission is

$$\left|\widehat{\text{Bias}}\right|=\sqrt{\frac{R_{Y\sim Z\mid D,\mathbf{X}}^2{R}_{D\sim Z\mid \mathbf{X}}^2}{1-R_{D\sim Z\mid \mathbf{X}}^2}}\cdot \frac{\text{sd}\left(Y^{\perp (\mathbf{X},D)}\right)}{\text{sd}\left(D^{\perp \mathbf{X}}\right)}.$$

They then define

${\text{RV}}_q=\frac{1}{2}\left[\sqrt{f_q^4+4f_q^2}-f_q^2\right]$

where $f_q:=q\left|f_{Y\sim D\mid \mathbf{X}}\right|$ where $f_{Y\sim D\mid \mathbf{X}}$ is the partial Cohen’s $f$ of the treatment with the outcome, and $q$ is the proportion of reduction on the treatment coefficient $\tau$ that would be deemed problematic.

Partial Identification

the ATE can be decomposed as

$$\text{ATE}=\mathbb{E}\left[Y(1)\right]-\mathbb{E}\left[Y(0)\right]=\mathbb{E}\left[Y_i(1)\mid D_i=1\right]\mathbf{Pr}\left(D_i=1\right)+\mathbb{E}\left[Y_i(1)\mid D_i=0\right]\mathbf{Pr}\left(D_i=0\right)-\mathbb{E}\left[Y_i(0)\mid D_i=1\right]\mathbf{Pr}\left(D_i=1\right)-\mathbb{E}\left[Y_i(0)\mid D_i=0\right]\mathbf{Pr}\left(D_i=0\right).$$

The terms in red are counterfactual outcomes for which the data contains no information. Bounding approaches involve estimators for these missing quantities.

Suppose all we know is $Y^d\in [0,1]$

w.l.o.g. given bounded support $[\underline{Y},\overline{Y}]$, we can always min-max rescale to $\frac{Y-\underline{Y}}{\overline{Y}-\underline{Y}}$

$$\mathbb{E}\left[Y^1-Y^0\right]\in \left[\left\{\mathbb{E}\left[Y\mid D=1\right]\mathbf{Pr}\left(D=1\right)-\mathbb{E}\left[Y\mid D=0\right]\left(1-\mathbf{Pr}\left(D=1\right)\right)\right\}-\mathbf{Pr}\left(D=1\right),\left\{\mathbb{E}\left[Y\mid D=1\right]\mathbf{Pr}\left(D=1\right)-\mathbb{E}\left[Y\mid D=0\right]\left(1-\mathbf{Pr}\left(D=1\right)\right)\right\}+\left(1-\mathbf{Pr}\left(D=1\right)\right)\right].$$

Width of possible interval learnable from data is $[0,1]$ at largest, $[-1,0]$ at smallest, so worst case interval always contains 0. Need theory/assumptions to even get the sign right.

Assume bounded support for the outcome. Replace missing values with maximum ($y^{UB}$) or minimum ($y^{LB}$) of support. These are worst-case bounds and yield intervals that are basically uninformative.

$$\mathbb{E}{\left[Y(1)\right]}^{UB}=\mathbb{E}\left[Y\mid D=1\right]\mathbf{Pr}\left(D=1\right)+y^{UB}\mathbf{Pr}\left(D=0\right).$$ $$\mathbb{E}{\left[Y(1)\right]}^{LB}=\mathbb{E}\left[Y\mid D=1\right]\mathbf{Pr}\left(D=1\right)+y^{LB}\mathbf{Pr}\left(D=0\right).$$ $$\mathbb{E}{\left[Y(0)\right]}^{UB}=y^{UB}\mathbf{Pr}\left(D=1\right)+\mathbb{E}\left[Y\mid D=0\right]\mathbf{Pr}\left(D=0\right).$$ $$\mathbb{E}{\left[Y(0)\right]}^{LB}=y^{LB}\mathbf{Pr}\left(D=1\right)+\mathbb{E}\left[Y\mid D=0\right]\mathbf{Pr}\left(D=0\right).$$

And denote ${\Delta }^{UB}:=\mathbb{E}{\left[Y(1)\right]}^{UB}-\mathbb{E}{\left[Y(0)\right]}^{LB}$ ${\Delta }^{LB}:=\mathbb{E}{\left[Y(1)\right]}^{LB}-\mathbb{E}{\left[Y(0)\right]}^{UB}$

Monotone Treatment Response: assume mean potential outcome under treatment cannot be lower than under control $\mathbb{E}\left[Y(1)\right]\ge \mathbb{E}\left[Y(0)\right]=\Delta \ge 0$. Then

${\Delta }^{LB}=\max (\mathbb{E}{\left[Y(1)\right]}^{LB}-\mathbb{E}{\left[Y(0)\right]}^{UB},0)$

Monotone Treatment Selection: subjects select themselves into treatment in a way the mean potential outcomes of the treatment and control groups can be ordered. Positive MTS implies $\mathbb{E}\left[Y(1)|D=1\right]\ge \mathbb{E}\left[Y(1)|D=0\right]$ and $\mathbb{E}\left[Y(0)|D=1\right]\ge \mathbb{E}\left[Y(0)|D=0\right]$. This implies $\mathbb{E}{\left[Y(0)\right]}^{LB}=\mathbb{E}\left[Y|D=0\right]$ and $\mathbb{E}{\left[Y(1)\right]}^{UB}=\mathbb{E}\left[Y|D=1\right]$

Let ${\tau }_i:=Y_{1i}-y_{0i}$ denote the treatment effect and $\mathbb{F}$ denote its distribution, and let ${\mathbb{F}}_1,{\mathbb{F}}_0$ denote the distributions of outcomes for the two potential outcomes. Then, ${\mathbb{F}}^L(b)\le \mathbb{F}(b)\le {\mathbb{F}}^U(b)$ where

$${\mathbb{F}}^L(b)=\max \left\{\underset{y}{\max }{\mathbb{F}}_1(y)-{\mathbb{F}}_0(y-b),0\right\}$$ $${\mathbb{F}}^U(b)=1+\min \left\{\underset{y}{\min }{\mathbb{F}}_1(y)-{\mathbb{F}}_0(y-b),0\right\}$$

Instrumental Variables

SOO Fails/$\mathbb{E}\left[X_i{ϵ}_i\right]\ne 0$ because of OVB, then ${\hat{\beta }}_{OLS}$ is no longer consistent. Use $Z$ as instrument for $D$ which isolates variation unrelated to the omitted variable.

Traditional IV Framework (Constant Treatment Effects)

Setup

Second Stage: $Y={\alpha }_0+{\alpha }_{1}D+u_2$

First Stage: $D={\pi }_0+{\pi }_{1}Z+u_1$

Reduced Form: $$\begin{aligned} Y & = \gamma_0 + \gamma_1 Z + u_3 \ & = \alpha_0 + \alpha_1 ( \pi_0 + \pi_1 Z + u_1) + u_2 \ & = (\alpha_0 + \alpha_1 \pi_0) + \underbrace{\hlred{(\alpha_1 \pi_1)}}_{\gamma_1} Z + (\alpha_1 u_1 + u_2)

\end{aligned}$$

Exogeneity (as good as random conditional on covariates): $\mathrm{Cov}\left(u_1,Z\right)=0$

Exclusion Restriction: $\mathrm{Cov}\left(u_2,D\right)=0$, $Z$ has no effect on $Y$ except through $D$.

Relevance: $Z$ affects $D$

With the above assumptions, we can write

${\hat{\beta }}_{IV}={\left({\mathbf{Z}}^{\prime }\mathbf{X}\right)}^{-1}{\mathbf{Z}}^{\prime }\mathbf{y}$

This is equivalent to

With binary treatment and binary instrument, one can write the IV effect as

$${\alpha }_1=\frac{{\gamma }_1}{{\pi }_1}=\frac{\mathrm{Cov}\left(Y,Z\right)}{\mathrm{Cov}\left(D,Z\right)}=\frac{\mathbb{E}\left[Y\mid Z=1\right]-\mathbb{E}\left[Y\mid Z=0\right]}{\mathbb{E}\left[D\mid Z=1\right]-\mathbb{E}\left[D\mid Z=0\right]}.$$

With multiple instruments or endogenous variables,

$${\hat{\alpha }}_{2SLS}={\left({\mathbf{X}}^{\prime }{\mathbf{P}}_z\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\prime }{\mathbf{P}}_z\mathbf{y}$$

where ${\mathbf{P}}_z=\mathbf{Z}{\left({\mathbf{Z}}^{\prime }\mathbf{Z}\right)}^{-1}{\mathbf{Z}}^{\prime }$ is $\mathbf{X}$ projected in the column space of $\mathbf{Z}.$

$${\hat{\alpha }}_k={\left({\mathbf{X}}^{\prime }(\mathbf{I}-k{\mathbf{P}}_z)\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\prime }(\mathbf{I}-k{\mathbf{P}}_z)\mathbf{y}$$

which nests 2SLS, LIML, and Fuller’s estimator as special cases. Specifically,

$k=0⟹{\hat{\alpha }}_k$ is OLS

$k=1⟹{\hat{\alpha }}_k$ is 2SLS

$k=k_{\text{LIML}}⟹{\hat{\alpha }}_k$ is LIML

$k=k_{\text{LIML}}-\frac{b}{n-L-p};b>0⟹{\hat{\alpha }}_k$ is Fuller’s estimator

here, $k_{\text{LIML}}$ is the minimum value of $k$ that satisfies

$$\det \left(\begin{array}{cc}{\mathbf{y}}^{\top }(\mathbf{I}-k{\mathbf{P}}_z)\mathbf{y}& {\mathbf{y}}^{\top }(\mathbf{I}-k{\mathbf{P}}_z)\mathbf{X}\\ {\mathbf{X}}^{\top }(\mathbf{I}-k{\mathbf{P}}_z)\mathbf{y}& {\mathbf{X}}^{\top }(\mathbf{I}-k{\mathbf{P}}_z)\mathbf{X}\end{array}\right)=0$$

Implemented in ivmodel, which takes model fits from AER::ivreg and computes LIML / k-class estimates.

Asymptotically, all $k-$ class estimators are consistent for $\alpha$ when $k\to 1,n\to \infty$.

Inference

Under homoscedasticity, $\mathbb{V}\left[{\hat{\alpha }}_{2SLS}\right]={\sigma }^2{\left({\mathbf{X}}^{\prime }{\mathbf{P}}_z\mathbf{X}\right)}^{-1}$

Under heteroskedasticity, $\mathbb{V}({\hat{\beta }}_{IV})={\left({\mathbf{Z}}^{\prime }{\mathbf{P}}_z\mathbf{X}\right)}^{-1}{\mathbf{P}}_z{\mathbf{X}}^{\prime }\hat{\Omega }{\mathbf{P}}_z\mathbf{X}{\left({\mathbf{X}}^{\prime }{\mathbf{P}}_z\mathbf{X}\right)}^{-1};\hat{\Omega }=\text{Diag}[{\hat{u}}_i^2]$

Test statistic and null distribution $H:=\frac{({\hat{\beta }}_{2sls}-{\hat{\beta }}_{ols}{)}^2}{\hat{V}({\hat{\beta }}_{2sls})-\hat{V}({\hat{\beta }}_{ols})}\sim {\chi }_1^2$

Equivalently, Assuming the instrument $Z$ is valid, we can test for whether $\mathbf{x}$ is endogenous by estimating the following regression

$y_i={\mathbf{Z}}_i^{\prime }\mathbf{\pi }+x_i{\alpha }_1+{\hat{v}}_i{\rho }_1+{ϵ}_i$

where $\hat{v}$ are the (fitted) residuals from estimating the first stage regression $x_i={\mathbf{Z}}_i^{\prime }\mathbf{\psi }+v_i$. A standard t-test for $\rho$ tests whether $x$ is exogenous assuming ${\mathbf{Z}}_i$ is a valid set of instruments. [means this test is not that useful in practice]

Weak Instruments

$$\mathrm{plim}{\alpha }_{IV}=\frac{\mathrm{Cov}\left(Y,Z\right)}{\mathrm{Cov}\left(Z,D\right)}+\frac{\mathrm{Cov}\left(Z,u_2\right)}{\mathrm{Cov}\left(Z,D\right)}={\alpha }_D+\frac{\mathrm{Cov}\left(Z,u_2\right)}{\mathrm{Cov}\left(Z,D\right)}.$$

Second term non-zero if instrument is not exogenous. Let ${\sigma }_{u1,u_2}=\mathrm{Cov}\left(u_1,u_2\right)$ and ${\sigma }_{u_1}^2=\mathbb{V}\left[u_2\right]$ [variance of first stage error] and $F$ be F statistic of the first-stage. Then, bias in IV is

$$\mathbb{E}\left[{\hat{\alpha }}_{IV}-\alpha \right]=\frac{{\sigma }_{u_1,u_2}}{{\sigma }_{u_2}^2}\frac{1}{F+1}$$

If first stage is weak, bias approaches $\frac{{\sigma }_{u_1,u_2}}{{\sigma }_{u_2}^2}$. As $F\to \infty$, $B_{IV}\to 0$.

When instruments are weak, AR confidence intervals are preferable to eyeballing F-statistics. Let $\mathbf{M}$ be a $n\times 2$ matrix of $\left(\begin{array}{cc}\mathbf{y}& \mathbf{X}\end{array}\right)$, and let $a_0=({\beta }_0,1)$, $b_0=(1,-{\beta }_0)$ (where ${\beta }_0$ is typically 0), and

$$\hat{\mathbf{\Sigma }}=\frac{{\mathbf{M}}^{\top }{\mathbf{P}}_z\mathbf{M}}{n-L-p}$$

be an estimator for the covariance matrix for the errors.

and let $\hat{\mathbf{s}},\hat{\mathbf{t}}$ be two-dimensional vectors defined as

$$\hat{\mathbf{s}}:=({\mathbf{Z}}^{\top }\mathbf{Z}{)}^{\frac{1}{2}}{\mathbf{Z}}^{\top }\mathbf{M}{b}_0(b_0^{\top }\hat{\mathbf{\Sigma }}{b}_0{)}^{-\frac{1}{2}}$$

and

$$\hat{\mathbf{t}}:=({\mathbf{Z}}^{\top }\mathbf{Z}{)}^{\frac{1}{2}}{\mathbf{Z}}^{\top }\mathbf{M}{\hat{\mathbf{\Sigma }}}^{-1}{a}_0(a_0^{\top }\hat{\mathbf{\Sigma }}{a}_0{)}^{-\frac{1}{2}}$$

Define the scalars ${\hat{Q}}_1={\hat{\mathbf{s}}}^{\top }\hat{\mathbf{s}},{\hat{Q}}_2={\hat{\mathbf{s}}}^{\top }\hat{\mathbf{t}},Q_3={\hat{\mathbf{t}}}^{\top }\hat{\mathbf{t}}$

based on these scalars, two tests that are fully robust to weak instruments for testing $H_0:\beta ={\beta }_0$ - Anderson Rubin test (AR1949) and Conditional Likelihood Test (Moriera 2003)

$$AR\left({\beta }_0\right)=\frac{{\hat{Q}}_1}{L}$$ $$CLR\left({\beta }_0\right)=\frac{1}{2}\left({\hat{Q}}_1-{\hat{Q}}_3\right)+\frac{1}{2}\sqrt{{\left({\hat{Q}}_1+{\hat{Q}}_3\right)}^2-4\left({\hat{Q}}_1{\hat{Q}}_3-{\hat{Q}}_2^2\right)}$$

IV with Heterogeneous Treatment Effects / LATE Theorem

binary instrument $Z_i\in \{0,1\}$

binary treatment $D_z\in \{0,1\}$ is potential treatment status given $Z=z$

potential outcomes: $Y_i(D,Z)=\{Y(1,1),Y(1,0),Y(0,1),Y(0,0)\}$

heterogeneous treatment effects ${\beta }_i=Y_i(1)-Y_i(0)$

Compliers: $D_1>D_0$, $D_0=0$, $D_1=1$

Always takers: $D_0=D_1=1$

Never Takers : $D_0=D_1=0$

Defiers: $D_1<D_0$

A1: Independence of Instrument : $\{Y_0,Y_1,D_0,D_1\}\perp \perp Z$

A2: Exclusion restriction : $Y_i(d,0)=Y_i(d,1)\equiv Y_{di}\text{for }d=0,1$

A3: First Stage: $\mathbb{E}\left[D_{1i}-D_{0i}\right]\ne 0$

A4: Monotonicity / No defiers: $D_{1i}-D_{0i}\ge 0\forall i$ or vice versa

Under A1-A4,

$${\alpha }_{IV}=\frac{\mathbb{E}\left[Y|Z=1\right]-\mathbb{E}\left[Y|Z=0\right]}{\mathbb{E}\left[D|Z=1\right]-\mathbb{E}\left[D|Z=0\right]}=\mathbb{E}\left[Y_{1i}-Y_{0i}|D_{1i}>D_{0i}\right]$$

If A1:A4 are satisfied, the IV estimate is the Local Average Treatment Effect for the compliers.

$$\text{LATE}=\text{ATE}+\frac{\mathrm{Cov}\left({\beta }_{1i},{\pi }_{1i}\right)}{\mathbb{E}\left[{\pi }_{1i}\right]}$$

So, LATE is a weighted average for people with large ${\pi }_{1i}$; i.e. treatment effect for those whose probability of treatment is most influenced by $Z_i$.

IV in Randomized Trials with one-sided noncompliance. Conditional on A1:A4 holding, and $\mathbb{E}\left[D|Z_i=0\right]=\mathbf{Pr}\left(D=1|Z=0\right)=0$. Then,

$$\frac{\mathbb{E}\left[Y|Z=1\right]-\mathbb{E}\left[Y|Z=0\right]}{\mathbf{Pr}\left(D=1|Z=1\right)}=\frac{\text{ITT}}{\text{Compliance}}=\mathbb{E}\left[Y_1-Y_0|D=1\right]=ATT$$

Precision for LATE Estimation

$S{E}_{\widehat{LATE}}\approx \frac{S{E}_{\widehat{ITT}}}{\text{Compliance}}$

Characterising Compliers

PO Model of IV allows for heterogeneous treatment effects but does not formally identify LATE conditional on X.

extends methods by allowing the treatment inducer to be randomized conditionally on the covariates and by allowing the outcome to depend on the covariates besides the treatment intake. The paper also provided semiparametric estimations of the probability of receiving the treatment inducement, which helps to identify the treatment effects in a more robust way.

Need the following assumptions (all conditional on $\mathbf{X})$:

Independence of instrument: $Z\perp \perp (D(z),Y(z^{\prime },d))|\mathbf{X}\forall z,z^{\prime },d\in \left\{0,1\right\}$: SOO w.r.t. instrument. Exclusion restriction: $\mathbf{Pr}\left(Y(1,d)=Y(0,d)=Y(d)|\mathbf{X}\right)=1$ Monotonicity: $\mathbf{Pr}\left(D(1)\ge D(0)|\mathbf{X}\right)=1$ First Stage: $\mathbb{E}\left[D|Z=1,\mathbf{X}\right]-\mathbb{E}\left[D|Z=0,\mathbf{X}\right]\ne 0$ Common Support : $0<\mathbf{Pr}\left(Z=1,\mathbf{X}\right)<1$

Specifically, when the treatment inducer Z is as good as randomized after conditioning on covariates X, Abadie proposed a two-stage procedure to estimate treatment effects.

Estimate the probability of receiving the treatment inducement $P(Z=1|\mathbf{X})$ (preferably using a semiparametric estimator) in order to provide a set of pseudo-weights.

Second, the pseudo-weights are used to estimate the local average response function (LARF) of the outcome conditional on the treatment and covariates.

The estimated coefficient for the treatment intake D reflects the conditional treatment effect.

Given monotonicity, we can identify the proportion of compliers, never-takers, and always-takers respectively.

$${\pi }_{\text{compliers}}=\mathbf{Pr}\left(D_1>D_0|\mathbf{X}\right)=\mathbb{E}\left[D|\mathbf{X},Z=1\right]-\mathbb{E}\left[D|\mathbf{X},Z=0\right]$$ $${\pi }_{\text{always-takers}}=\mathbf{Pr}\left(D_1=D_0=1|\mathbf{X}\right)=\mathbb{E}\left[D|\mathbf{X},Z=0\right]$$ $${\pi }_{\text{never-takers}}=\mathbf{Pr}\left(D_1=D_0=0|\mathbf{X}\right)=1-\mathbb{E}\left[D|\mathbf{X},Z=1\right]$$

If nobody in the treatment group has access to the treatment (i.e. $\mathbb{E}\left[D|Z=0\right]=0$), the $\text{LATE}=\text{ATT}$.

By Bayes rule,

$$\mathbf{Pr}\left(D_1>D_0|D=1\right)=\frac{\mathbf{Pr}\left(D=1|D_1>D_0\right)\mathbf{Pr}\left(D_1>D_0\right)}{\mathbf{Pr}\left(D=1\right)}=\frac{\mathbf{Pr}\left(Z=1\right)[\mathbb{E}\left[D|Z=1\right]-\mathbb{E}\left[D|Z=0\right]]}{\mathbf{Pr}\left(D=1\right)}$$

Suppose assumptions of LATE thm hold conditional on covariates $X$. Let $g(\cdot )$ be any measurable real function of $Y,D,X$ with finite expectation. We can show that the expectation of $g$ is a weighted sum of the expectation in the three groups

$$\mathbb{E}\left[g|\mathbf{X}\right]=\underset{\text{Compliers}}{\underbrace{\mathbb{E}\left[g|\mathbf{X},D_1>D_0\right]\mathbf{Pr}\left(D_1>D_0|\mathbf{X}\right)}}+\underset{\text{Always takers}}{\underbrace{\mathbb{E}\left[g|\mathbf{X},D_1=D_0=1\right]\mathbf{Pr}\left(D_1=D_0=1|\mathbf{X}\right)}}+\underset{\text{Never Takers}}{\underbrace{\mathbb{E}\left[g|\mathbf{X},D_1=D_0=0\right]\mathbf{Pr}\left(D_1=D_0=0\mid \mathbf{X}\right)}}$$

Rearranging terms gives us

Then,

$$\mathbb{E}\left[g(Y,D,\mathbf{X})|D_1>D_0\right]=\frac{\mathbb{E}\left[\kappa \cdot g(Y,D,\mathbf{X})\right]}{\mathbf{Pr}\left(D_1>D_0\right)}=\frac{\mathbb{E}\left[\kappa \cdot g(Y,D,\mathbf{X})\right]}{\mathbb{E}\left[\kappa \right]}$$

where

$${\kappa }_i=1-\frac{D(1-Z)}{1-\mathbf{Pr}\left(Z=1\mid \mathbf{X}\right)}-\frac{(1-D)Z}{\mathbf{Pr}\left(Z=1\mid \mathbf{X}\right)}$$

This result can be applied to any characteristic or outcome and get its mean for compliers by removing the means for never and always takers. [p 181-183] provides overview of estimation. Trick is to construct a weighting scheme with positive weights so that ${\kappa }_i$, which is negative for always-takers and never-takers.

To compute $\kappa$, we need $\mathbf{Pr}\left(Z=1|\mathbf{X}\right)$, which can be computed using a standard logit/probit or a power-series.

Standard example: average covariate value among compliers:

$\mathbb{E}\left[X|D_1>D_0\right]=\frac{\mathbb{E}\left[\kappa X\right]}{\mathbb{E}\left[\kappa \right]}$

is the weighted average of covariate $X$ using Kappa weights.

Likelihood that Complier has a given value of (Bernoulli distributed) characteristic X relative to the rest of the population is given by

$$\frac{\mathbb{E}\left[D|Z=1,X=1\right]-\mathbb{E}\left[D|Z=0,X=1\right]}{\mathbb{E}\left[D|Z=1\right]-\mathbb{E}\left[D|Z=0\right]}=\frac{\text{FS in Subgroup}}{\text{Overall FS}}$$

Assume A1-A4 from LATE. Generalise $D$ to take values in the set$\left\{0,1,\dots ,\check{D}\right\}$; Let $Y_{di}:=f_i(d)$ denote the potential (or latent) outcome for person $i$ for treatment level $d$. Then,

$$\frac{\mathbb{E}\left[Y|Z=1\right]-\mathbb{E}\left[Y|Z=0\right]}{\mathbb{E}\left[D|Z=1\right]-\mathbb{E}\left[D|Z=0\right]}=\sum _{d=1}^{\check{D}}{\omega }_d\mathbb{E}\left[Y_{di}-Y_{d-1,i}|d_{1i}\ge d>d_{0i}\right]$$

where the weights

$${\omega }_d=\frac{\mathbf{Pr}\left(d_{1i}>d>d_{0i}\right)}{{\sum }_{j=1}^{\check{D}}\mathbf{Pr}\left(d_{1i}\ge j\ge d_{0i}\right)}$$

are non-negative and sum to 1.

CEF of $Y|X,D$ for the subpopulation of compliers: $\mathbb{E}\left[Y|X,D,D_1>D_0\right]$

$\mathbb{E}\left[Y|X,D,D_1>D_0\right]=\frac{\mathbb{E}\left[\kappa Y|X,D\right]}{\mathbb{E}\left[\kappa \right]}$

Estimate $\kappa$

Estimate $\mathbb{E}\left[Y|X,D\right]$ in the whole population, weighting by $\kappa$

implemented in LARF::larf in R.

Treatment is $W$. First define two additional quantities

$P_{A,C,i}:=\mathbf{Pr}\left(W_1>W_0\cup W_0=1|{\mathbf{X}}_{\mathbf{i}}=\mathbf{x}\right)=\mathbb{F}({\mathbf{x}}_i^{\prime }{\mathbf{\theta }}_{A,C})$ is the conditional probability that unit $i$ is either a complier *or* an always taker assume that this probability is a function of covariates ${\mathbf{X}}_i$ , with corresponding parameter vector ${\mathbf{\theta }}_{A,C}$ and CDF $\mathbb{F}$ that transforms it to the probability scale [taken to be the normal CDF $\Phi$ henceforth, but can be relaxed] $P_{A|A,C,i}:=\mathbf{Pr}\left(W_0=1|W_1>W_0\cup W_0=1,{\mathbf{X}}_{\mathbf{i}}=\mathbf{x}\right)=\mathbb{F}({\mathbf{x}}_i^{\prime }{\mathbf{\varphi }}_{A|A,C})$ is the conditional probability that unit $i$ is an always taker *conditional* on being either a complier or a never taker assume that this probability is a function of covariates with corresponding covariate vector ${\mathbf{\varphi }}_{A|A,C}$

Next, they note that the probability of treatment for stratum ${\mathbf{X}}_{\mathbf{i}}={\mathbf{x}}_{\mathbf{i}}$ can be written as

$$\mathbf{Pr}\left(W=1\mid {\mathbf{X}}_{\mathbf{i}}={\mathbf{x}}_{\mathbf{i}}\right)=\stackrel{\text{Compliers assigned to treatment}}{\overbrace{\mathbf{Pr}\left(W_1>W_0\mid {\mathbf{X}}_{\mathbf{i}}={\mathbf{x}}_{\mathbf{i}}\right)Z_i}}+\stackrel{\text{Always takers}}{\overbrace{\mathbf{Pr}\left(W_0=1\mid {\mathbf{X}}_{\mathbf{i}}={\mathbf{x}}_{\mathbf{i}}\right)}}$$

Using the two conditional probabilties defined above, this can be written as

$$\mathbf{Pr}\left(W=1\mid {\mathbf{X}}_{\mathbf{i}}={\mathbf{x}}_{\mathbf{i}}\right)=P_{A,C,i}(1-P_{A|A,C,i})Z_i+P_{A,C,i}{P}_{A|A,C,i}$$

which, for binary treatment $W_i$ lets us write a Bernoulli likelihood for an observation

$${\ell }_i(P_{A|A,C,i},P_{A,C,i}|W,Z)=(P_{A,C,i}(1-P_{A|A,C,i})Z_i+P_{A,C,i}{P}_{A|A,C,i}{)}^{W_i}(1-P_{A,C,i}(1-P_{A|A,C,i})Z_i-P_{A,C,i}{P}_{A|A,C,i}{)}^{1-W_i}$$

Plugging in the definitions of $P_{A,C,i}$ and $P_{A|A,C,i}$ gives us the likelihood and its argmax defines the solution for ${\hat{\mathbf{\theta }}}_{A,C}$ and ${\hat{\mathbf{\varphi }}}_{A|A,C}$. This is generically a difficult optimisation problem and improving its computation is a promising avenue for future research.

$$\mathcal{L}(P_{A|A,C,i},P_{A,C,i}|W,Z)=\prod _{i=1}^n{\left(\mathbb{F}({\mathbf{x}}_i^{\prime }{\mathbf{\theta }}_{A,C})(1-\mathbb{F}({\mathbf{x}}_i^{\prime }{\mathbf{\varphi }}_{A|A,C})Z_i+\mathbb{F}({\mathbf{x}}_i^{\prime }{\mathbf{\theta }}_{A,C})\mathbb{F}({\mathbf{x}}_i^{\prime }{\mathbf{\varphi }}_{A|A,C}))\right)}^{W_i}{\left((1-\mathbb{F}({\mathbf{x}}_i^{\prime }{\mathbf{\theta }}_{A,C}))(1-\mathbb{F}({\mathbf{x}}_i^{\prime }{\mathbf{\varphi }}_{A|A,C})Z_i-\mathbb{F}({\mathbf{x}}_i^{\prime }{\mathbf{\theta }}_{A,C})\mathbb{F}({\mathbf{x}}_i^{\prime }{\mathbf{\varphi }}_{A|A,C}))\right)}^{1-W_i}$$

The maximum likelihood estimates of the two parameter vectors can be plugged into $\mathbb{F}$ to compute individual compliance scores

$${\hat{P}}_{C,i}=\mathbf{Pr}\left(W_1>W_0\mid {\mathbf{X}}_{\mathbf{i}}={\mathbf{x}}_{\mathbf{i}}\right)=\mathbb{F}({\mathbf{x}}_i^{\prime }{\hat{\mathbf{\theta }}}_{A,C})\left[1-\mathbb{F}({\mathbf{x}}_i^{\prime }{\hat{\mathbf{\varphi }}}_{A\mid A,C})\right]$$

The inverse compliance score weighted estimator for the ATE with weights ${\omega }_{C,i}:=1/{\hat{P}}_{C,i}$ is then

$${\hat{\tau }}_{\text{ICSW}}^{\text{ATE}}=\frac{\frac{{\sum }_{i=1}^n{\hat{\omega }}_{Ci}{Z}_i{Y}_i}{{\sum }_{i=1}^n{\hat{\omega }}_{Ci}{Z}_i}-\frac{{\sum }_{i=1}^n{\hat{\omega }}_{Ci}\left(1-Z_i\right)Y_i}{{\sum }_{i=1}^n{\hat{\omega }}_{Ci}\left(1-Z_i\right)}}{\frac{{\sum }_{i=1}^n{\hat{\omega }}_{Ci}{Z}_i{W}_i}{{\sum }_{i=1}^n{\hat{\omega }}_{Ci}{Z}_i}-\frac{{\sum }_{i=1}^n{\hat{\omega }}_{Ci}\left(1-Z_i\right)W_i}{{\sum }_{i=1}^n{\hat{\omega }}_{Ci}\left(1-Z_i\right)}}$$

which is a weighted version of the familiar Wald estimator with a Hajek correction that normalises each expectation by the sum of weights in that treatment group.

Shift Share / Bartik Instruments

SSIV setting from [notation and exposition from PGP’s slides]. We want to estimate the causal effect or structural parameter $\tau$ in

$y_l=\tau w_l+{\gamma }^{\top }{\mathbf{x}}_l+{\epsilon }_l$

where $\mathrm{Cov}\left({\epsilon }_l,w_l\right)\ne 0$ because the ‘treatment’ $w_l$ is typically a change in an economic quantity (e.g. employment) that is correlated with unobserved shocks to the outcome $y_l$ (e.g. wages). $l$ indexes locations.

An accounting identity that decomposes the treatment is

$$w_l=\sum _{k=1}^K\stackrel{\text{Location-Industry Shares}}{\overbrace{z_{lk}}}\underset{\text{Location-Industry Shifts}}{\underbrace{g_{lk}}}$$

where $k$ indexes industries. 2nd accounting identity for location-industry shifts is

$$\underset{\text{Location-industry}}{\underbrace{g_{lk}}}=\underset{\text{industry}}{\underbrace{g_k}}+\underset{\text{location-industry shocks: unobserved}}{\underbrace{{\tilde{g}}_{lk}}}$$

As a GMM system

$$y_{lt}={\mathbf{D}}_{lt}^{\top }{\mathbf{\beta }}_0+\tau w_{lt}+{\epsilon }_{lt}$$ $$w_{lt}={\mathbf{D}}_{lt}^{\top }{\mathbf{\gamma }}_0+\psi B_{lt}+{\eta }_{lt}$$ $$g_{lkt}=g_{kt}+{\tilde{g}}_{lkt}$$ $$B_{lt}=\sum _{k=1}^K{z}_{lk0}{g}_{kt}$$

${\mathbf{D}}_{lt}$ denotes exogenous controls and fixed effects.

$${\left\{{\left\{w_{lt},{\mathbf{D}}_{lt},{\epsilon }_{lt}\right\}}_{t=1}^T\right\}}_{l=1}^L\text{ are IID},L\to \infty$$

Under constant $\tau$, need Exogeneity $\mathbb{E}\left[B_{lt}{\epsilon }_{lt}\mid {\mathbf{D}}_{lt}\right]=0$ Relevance $\mathrm{Cov}\left(B_{lt},w_{lt}\mid {\mathbf{D}}_{lt}\right)\ne 0$

$${\hat{\tau }}_{\text{Bartik}}=\frac{{\sum }_{l=1}^L{\sum }_{t=1}^T{\sum }_{k=1}^K{z}_{lkt}{g}_{kt}{y}_{lt}^{\perp }}{{\sum }_{l=1}^L{\sum }_{t=1}^T{\sum }_{k=1}^K{z}_{lkt}{g}_{kt}{w}_{lt}^{\perp }}$$

‘shares’: focus on $z_{lk0}$ : Analogy to DiD: ${\Delta }_{gt}=$ Changes in industry composition $g_{kt}$ ‘shifts’: focus on $g_{kt}$ : requires argument for why shocks are randomly assigned

${\hat{\tau }}_{\text{bartik}}={\sum }_k{\hat{\alpha }}_k{\hat{\tau }}_k$

with Rotemberg weight

${\hat{\alpha }}_k=\frac{g_k{Z}_k^{\top }W}{{\sum }_{k=1}^K{g}_k{Z}_k^{\top }W}$

Marginal Treatment Effects: Treatment effects under self selection

propose the marginal treatment effect (MTE) setup that generalises the IV approach for continuous instruments and nests many estimands (and is a generalisation of the Roy (1951) model). It also has a clearer treatment of self-selection.

Exposition based on . Define potential outcomes

$$Y_{0i}={\mu }_0({\mathbf{x}}_i)+{\upsilon }_{0i}$$ $$Y_{1i}={\mu }_1({\mathbf{x}}_i)+{\upsilon }_{1i}$$

where ${\mu }_j(\cdot )$ is the conditional mean function and ${\upsilon }_{ji}$ captures deviations, with $\mathbb{E}\left[{\upsilon }_{ji}|{\mathbf{x}}_i\right]=0$.

Treatment assignment assumes a weakly separable choice model

$$D_i^{\ast }={\mu }_d({\mathbf{x}}_i,z_i)+v_i$$ $$D_i=1_{D_i^{\ast }\ge 0}$$

where $d_i^{\ast }$ is the latent propensity to take the treatment, and is interpreted as the net gain from treatment since treatment is only taken up if $D_i^{\ast }\ge 0$. $z_i$ is an instrument. $v_i$ enters the selection equation negatively, and thus represents latent resistance to treatment.

The condition $D_i^{\ast }\ge 0$ can be rewritten as ${\mu }_d({\mathbf{x}}_i,z_i)\ge v_i$. Applying the CDF of $v$ ${\mathbb{F}}_v$ to both sides yields

$${\mathbb{F}}_v({\mu }_d({\mathbf{x}}_i,z_i))\ge {\mathbb{F}}_v(v_i)$$

Define $P({\mathbf{x}}_i,z_i)=:{\mathbb{F}}_v({\mu }_d({\mathbf{x}}_i,z_i))$ and ${\upsilon }_{di}=:{\mathbb{F}}_v(v_i)$.

Both RHS and LHS are distributed on $[0,1]$. The treatment decision can now be written as

$D_i=1_{P({\mathbf{x}}_i,z_i)\ge {\upsilon }_{di}}$.

Now, we define treatment effects

$$Y_i=(1-D_i)Y_{0i}+D_i{Y}_{1i}=Y_{0i}+D_i(Y_{1i}-Y_{0i})={\mu }_0({\mathbf{x}}_i)+D_i\left[{\mu }_1({\mathbf{x}}_i)-{\mu }_0({\mathbf{x}}_i)+{\upsilon }_{1i}-{\upsilon }_{0i}\right]+{\upsilon }_{0i}$$

Aggregating over different parts of the covariate distribution yields different estimates.

$$\text{ATE}(\mathbf{x}):=\mathbb{E}\left[{\Delta }_i|{\mathbf{x}}_i=\mathbf{x}\right]={\mu }_1(\mathbf{x})-{\mu }_0(\mathbf{x})$$ $$\text{ATT}(\mathbf{x}):=\mathbb{E}\left[{\Delta }_i|{\mathbf{x}}_i=\mathbf{x},D_i=1\right]={\mu }_1(\mathbf{x})-{\mu }_0(\mathbf{x})+\mathbb{E}\left[{\upsilon }_{1i}-{\upsilon }_{0i}|D_i=1\right]$$ $$\text{ATU}(\mathbf{x}):=\mathbb{E}\left[{\Delta }_i|{\mathbf{x}}_i=\mathbf{x},D_i=0\right]={\mu }_1(\mathbf{x})-{\mu }_0(\mathbf{x})+\mathbb{E}\left[{\upsilon }_{1i}-{\upsilon }_{0i}|D_i=0\right]$$

Integrating these over $\mathbf{x}$ yields the conventional estimators. With self-selection based on $D_i=1_{D_i^{\ast }\ge 0}$, typically ATT > ATE > ATU.

The covariate-specific Wald estimator is

$$\text{Wald}(\mathbf{x})=\frac{\mathbb{E}\left[Y_i|z_i=1,{\mathbf{x}}_i=\mathbf{x}\right]-\mathbb{E}\left[Y_i|z_i=0,{\mathbf{x}}_i=\mathbf{x}\right]}{\mathbb{E}\left[D_i|z_i=1,{\mathbf{x}}_i=\mathbf{x}\right]-\mathbb{E}\left[D_i|z_i=0,{\mathbf{x}}_i=\mathbf{x}\right]}$$

Under the standard A1-A4 from AIR96,

$$\text{LATE}(\mathbf{x}):=\mathbb{E}\left[Y_{1i}-Y_{0i}|D_{1i}>D_{0i},{\mathbf{x}}_i=\mathbf{x}\right]={\mu }_1({\mathbf{x}}_i)-{\mu }_0({\mathbf{x}}_i)+\mathbb{E}\left[{\upsilon }_{1i}-{\upsilon }_{0i}|D_{1i}>D_{0i},{\mathbf{x}}_i=\mathbf{x}\right]$$

These can be aggregated using the ‘saturate and weight’ theorem (Angrist and Imbens)

$\text{IV}={\sum }_{\mathbf{x}\in \mathcal{X}}\omega (\mathbf{x})\text{LATE}(\mathbf{x})$

with weights

$$\omega ({\mathbf{x}}_i)=\frac{p_{\mathbf{x}}\mathbb{V}\left[{\hat{D}}_i\mid {\mathbf{x}}_i,z_i\right]}{\mathbb{V}\left[{\hat{D}}_i\right]}$$

For a continuous instrument, for a pair of instrument values $z,z^{\prime }$, $\text{LATE}(z,z^{\prime },\mathbf{x})=\mathbb{E}\left[Y_{1i}-Y_{0i}|D_{zi}>D_{z^{\prime }i},{\mathbf{x}}_i=\mathbf{x}\right]$.

$$\text{MTE}({\mathbf{x}}_i=\mathbf{x},V_i=v):=\mathbb{E}\left[Y_{1i}-Y_{0i}|{\mathbf{x}}_i=\mathbf{x},V_i=v\right]=\frac{\partial \mathbb{E}\left[Y_i|{\mathbf{x}}_i=\mathbf{x},p(Z,X)=p(z,x)\right]}{\partial p(z,x)}$$

MTE is defined as a continuum of treatment effects along the distribution of ${\upsilon }_D$.

Define two marginal treatment response (MTR) functions

$$m_0(u,\mathbf{x})=\mathbb{E}\left[Y_0|U=u,\mathbf{X}=\mathbf{x}\right]$$ $$m_1(u,\mathbf{x})=\mathbb{E}\left[Y_1|U=u,\mathbf{X}=\mathbf{x}\right]$$

Many useful parameters are identified using the following expression

$${\beta }^{\star }\equiv \mathbb{E}\left[{\int }_0^1{m}_0(u,X){\omega }_0^{\star }(u,X,Z)\mathrm{d}u\right]+\mathbb{E}\left[{\int }_0^1{m}_1(u,X){\omega }_1^{\star }(u,X,Z)\mathrm{d}u\right]$$

with weights specified in the figure below.

Parametric Model: Assuming joint normality for $U_0,U_1,V$,

$$\mathbb{E}\left[U_{0i}\mid D_i=0,X_i,Z_i\right]=\mathbb{E}\left[U_{0i}\mid V_i\ge \left(X_i,Z_i\right){\beta }_d,X_i,Z_i\right]={\rho }_0\left(\frac{\varphi \left(\left(X_i,Z_i\right){\beta }_d\right)}{1-\Phi \left(\left(X_i,Z_i\right){\beta }_d\right)}\right)$$ $$\mathbb{E}\left[U_{1i}\mid D_i=1,X_i,Z_i\right]=\mathbb{E}\left[U_{1i}\mid V_i<\left(X_i,Z_i\right){\beta }_d,X_i,Z_i\right]={\rho }_1\left(\frac{-\varphi \left(\left(X_i,Z_i\right){\beta }_d\right)}{\Phi \left(\left(X_i,Z_i\right){\beta }_d\right)}\right)$$

where ${\rho }_0$ is the correlation $\mathrm{Corr}\left(U_{0i},V_i\right)$, and ${\rho }_1=\mathrm{Corr}\left(U_{1i},V_i\right)$.

yields MTE estimator

$$\mathrm{MTE}\left(x,u_D\right)=\mathbb{E}\left[Y_{1i}-Y_{0i}\mid X_i=x,U_{Di}=u_D\right]=x\left({\beta }_1-{\beta }_0\right)+\left({\rho }_1-{\rho }_0\right){\Phi }^{-1}\left(u_D\right)$$

Let ${\tilde{\mathbf{x}}}_i={\mathbf{x}}_i-\overline{\mathbf{x}}$. Write

$$Y_i={\mathbf{x}}_i^{\top }\mathbf{\alpha }+D_i{\tilde{\mathbf{x}}}_i^{\prime }\mathbf{\theta }+D_i{\delta }_i+{\epsilon }_i$$

where ${\delta }_i$ is a random effect that captures treatment effect heterogeneity . We can rewrite this and by demeaning ${\delta }_i=\overline{\delta }-\widetilde{{\delta }_i}$.

$$Y_i={\mathbf{x}}_i^{\top }\mathbf{\alpha }+D_i{\tilde{\mathbf{x}}}_i^{\top }\mathbf{\theta }+D_i\overline{\delta }+{\upsilon }_{0i}$$

where $\overline{\delta }$ captures the ATE at means of $X$, which is the unconditional ATE under the linear specification.

Write the selection equation

$$D_i={\mathbf{x}}_i^{\top }{\pi }_1+z_i{\pi }_2+{\nu }_i$$

with $\mathbb{E}\left[{\nu }_i\mid {\mathbf{x}}_i,z_i\right]=0$.

Assumptions

$\mathbb{E}\left[{\epsilon }_i|{\nu }_i\right]=\eta {\nu }_i$: Conventional selection bias.

$\mathbb{E}\left[\widetilde{{\delta }_i}|{\nu }_i\right]=\psi {\nu }_i$: unobservable part of treatment effect ${\tilde{\delta }}_i$ depends linearly on the unobservables that affect treatment selection.

Including ${\hat{\nu }}_i$ and ${\hat{\nu }}_i{D}_i$ in the control-function outcome equation yields a consistent estimate of the ATE: $\overline{\delta }$.

High Dimensional IV selection

setup:

$$y_i=\tau d_i+{\mathbf{x}}_i^{\prime }\mathbf{\beta }+{\epsilon }_i$$ $$d_i={\mathbf{x}}_i^{\prime }{\mathbf{\gamma }}_0+{\mathbf{z}}_i^{\prime }{\mathbf{\delta }}_0+{\upsilon }_i$$

where

${\mathbf{x}}_i$ is a vector of $p_n^x$ exogenous controls, including a constant.

${\mathbf{z}}_i$ is a vector of $p_n^z$ instruments

$d_i$ is an endogenous variable

$p_n^x>>n$ and $p_n^z>>n$

Run (post)LASSO of $d_i$ on ${\mathbf{x}}_i,{\mathbf{z}}_i$ to obtain $\hat{\gamma },\hat{\delta }$

Run (post)LASSO of $y_i$ on ${\mathbf{x}}_i$ to get $\hat{\theta }$.

Run (post)LASSO of ${\hat{d}}_i={\mathbf{x}}_i^{\prime }\hat{\gamma }+{\mathbf{z}}_i^{\prime }\hat{\delta }$ on ${\mathbf{x}}_i$ to get $\hat{\vartheta }$.

Construct ${\hat{\rho }}_i^y:=y_i-{\mathbf{x}}_i\hat{\theta }$, ${\hat{\rho }}_i^d:=d_i-{\mathbf{x}}_i^{\prime }\hat{\vartheta }$ and ${\hat{v}}_i:={\mathbf{x}}_i\hat{\gamma }+{\mathbf{z}}_i^{\prime }\hat{\delta }+{\mathbf{x}}_i^{\prime }\hat{\vartheta }$.

Estimate $\hat{\tau }$ by using standard IV regression of ${\hat{\rho }}_i^y$ on ${\hat{\rho }}_i^d$ with ${\hat{v}}_i$ as instrument. Perform inference using score statistics or conventional heteroskedasticity-robust SEs.

implemented in hdm::rlassoIV(., select.X = T, select.Z = T).

Discussion in .

Principal Stratification

Treatment comparisons often need to be adjusted for post-treatment variables.

Binary treatment $Z_i\in \left\{0,1\right\}$. post-treatment Intermediate variable $S_i(z_i)\in \left\{0,1\right\}$, Outcome $Y_i\in \left\{0,1\right\}$. For each individual, the treatment assumes a single value, so only one of the two potential intermediate values are observed. Based on joint potential outcomes of the intermediate variable $(S_i,(0),S_i(1))$, we have 4 strata

$$00=\left\{i:S_i(0)=0,S_i(1)=0\right\}$$

Never takers.

$$10=\left\{i:S_i(0)=1,S_i(1)=0\right\}$$

Defiers.

$$01=\left\{i:S_i(0)=0,S_i(1)=1\right\}$$

Compliers.

$$11=\left\{i:S_i(0)=1,S_i(1)=1\right\}$$

Always takers.

The basic principal stratification $P_0$ w.r.t post treatment variable $S$ is the partition of units $i=1,\dots ,n$ such that, forall units in any set of $P_0$, all units have the same vector of $(S_i(0),S_i(1))$. The principal stratum $G_i\in \left\{00,10,01,11\right\}$ to which unit $i$ belongs is not affected by treatment assignment for any principal stratification, so can be considered pre-treatment.

Treatment Ignorability implies $(Y_i(0),Y_i(1))\perp \perp Z_i|S_i(0),S_i(1),\mathbf{X}$ (i.e. treatment and control units can be compared conditional on stratum)

Principal Causal Effect (PCE) ${\tau }_{s_0,s_1}:=\mathbb{E}\left[Y_i(1)-Y_i(0)|S_i(0)=s_0,S_i(1)=s_1\right]$

A common example is the

Complier Average Causal Effect (CACE) = Causal Effect on Principal Stratum of Compliers (AIR96)

$\text{CACE}=\mathbb{E}\left[Y_i(1)-Y_i(0)|S_i(0)=0,S_i(1)=1\right]$

Recall that $G_i=(S_0,S_1)$ concatenated. So, AIR96 in PS terms: Monotonicity: $S_1\ge S_0⟹\left\{G_i=10\right\}$ must be empty: no defiers. Exclusion: ${\tau }_{11}={\tau }_{00}$

Estimation under principal ignorability

Treatment ignorability $Z\perp \perp (S_0,S_1,Y_0,Y_1)\mid X$ monotonicity: $S_1\ge S_0⟹G_i=10$ is not allowed principal ignorability

$$\mathbb{E}\left[Y_1\mid G=11,X\right]=\mathbb{E}\left[Y_1\mid G=01,X\right]$$ $$\mathbb{E}\left[Y_0\mid G=00,X\right]=\mathbb{E}\left[Y_1\mid G=01,X\right]$$

	S = 0	S = 1
Z = 0	G = 00 or 01	G = 11
Z = 1	G = 00	G = 11 or 01

Disentangle mixture distribution within strata by assuming same conditional expectation across mixture components (complier, never taker, always taker).

Define nuisance functions:

Treatment probability: $\pi (X)=\mathbf{Pr}\left(Z=1\mid X\right)$ Principal Score: $e_g(X)=\mathbf{Pr}\left(G=g\mid X\right)$ identified by

$$e_{01}(X)=p_1(X)-p_0(X)$$ $$e_{00}(X)=1-p_1(X)$$ $$e_{11}(X)=p_0(X)$$

where $p_z(X)=\mathbf{Pr}\left(S=1\mid Z=z,X\right)$. Outcome mean: ${\mu }_{zs}(X)=\mathbb{E}\left[Y\mid Z=z,S=s,X\right]$.

Treatment Probability and Principal Score

$${\tau }_{01}=\mathbb{E}\left\{\frac{e_{01}(X)}{p_1-p_0}\frac{S}{p_1(X)}\frac{Z}{\pi (X)}Y\right\}-\mathbb{E}\left\{\frac{e_{01}(X)}{p_1-p_0}\frac{1-S}{1-p_0(X)}\frac{1-Z}{1-\pi (X)}Y\right\}$$ $${\tau }_{00}=\mathbb{E}\left\{\frac{1-S}{1-p_1}\frac{Z}{\pi (X)}Y\right\}-\mathbb{E}\left\{\frac{e_{00}(X)}{1-p_1}\frac{1-S}{1-p_0(X)}\frac{1-Z}{1-\pi (X)}Y\right\}$$ $${\tau }_{11}=\mathbb{E}\left\{\frac{e_{11}(X)}{p_0}\frac{S}{p_1(X)}\frac{Z}{\pi (X)}Y\right\}-\mathbb{E}\left\{\frac{S}{p_0}\frac{1-Z}{1-\pi (X)}Y\right\}$$

Treatment Probability and Outcome Mean

$${\tau }_{01}=\mathbb{E}\left[\frac{SZ/\pi (X)-S(1-Z)/\{1-\pi (X)\}}{p_1-p_0}\{{\mu }_{11}(X)-{\mu }_{00}(X)\}\right]$$ $${\tau }_{00}=\mathbb{E}\left[\frac{1-SZ/\pi (X)}{1-p_1}\{{\mu }_{10}(X)-{\mu }_{00}(X)\}\right]$$ $${\tau }_{11}=\mathbb{E}\left[\frac{S(1-Z)/\{1-\pi (X)\}}{p_0}\{{\mu }_{11}(X)-{\mu }_{01}(X)\}\right]$$

Principal Score and Outcome Mean

$${\tau }_{01}=\mathbb{E}\left[\frac{p_1(X)-p_0(x)}{p_1-p_0}\left\{{\mu }_{11}(X)-{\mu }_{00}(X)\right\}\right]$$ $${\tau }_{00}=\mathbb{E}\left[\frac{1-p_1(X)}{1-p_1}\left\{{\mu }_{10}(X)-{\mu }_{00}(X)\right\}\right]$$ $${\tau }_{11}=\mathbb{E}\left[\frac{p_0(X)}{p_0}\left\{{\mu }_{11}(X)-{\mu }_{01}(X)\right\}\right]$$

Direct and Indirect Effects via Principal Stratification

Direct effect of $Z$ conditional on $S$ exists if there is a causal effect of $Z$ on $Y$ for observations for whom the treatment does not affect selection $S$, i.e. principal strata $00,11$. This is a zero-first-stage sample in IV-terms.

The Indirect Effect is mediated through $S$.

Attrition as Selection Bias

Let $S$ denote a binary selection indicator for when $Y$ is observed. Let $S(1),S(0)$ denote potential selection states under treatment and nontreatment.

$S(1)=0,S(0)=0$ : never-selected $S(1)=1,S(0)=1$ : always selected $S(0)=0,S(1)=1$ : selection compliers $S(0)=1,S(1)=0$ : selection defiers (ruled out by Lee bounds)

Dominance assumption: $\mathbb{E}\left[Y(1)|S(1)=1,S(0)=1\right]\ge \mathbb{E}\left[Y(1)|S(1)=1,S(0)=0\right]$ and $\mathbb{E}\left[Y(0)|S(1)=1,S(0)=1\right]\ge \mathbb{E}\left[Y(0)|S(1)=1,S(0)=0\right]$. The average potential outcome of the always selected dominates that of compliers under either treatment state.

Then, Zhang and Rubin (2003) bounds are

$${\Delta }^{UB}=\mathbb{E}\left[Y|D=1,S=1,Y\ge y^{\ast }\right]-\mathbb{E}\left[Y|D=0,S=1\right]$$ $${\Delta }^{LB}=\mathbb{E}\left[Y|D=1,S=1\right]-\mathbb{E}\left[Y|D=0,S=1\right]$$

where $y^{\ast }$ is chosen such that the lowest outcomes among those with $D=1,S=1$ correspond to the share of compliers among those with $D=1,S=1$ are smaller than this value.

Assuming

randomisation: $\left\{Y(1),Y(0),S(0),S(1),\mathbf{X}\right\}\perp \perp D$ monotonicity: $S(1)\ge S(0)a.s.$

Lee (2009) focuses on the ATE among the always observed

$\mathbb{E}\left[Y(1)-Y(0)|S(0)=S(1)=1\right]$

The second quantity: $\mathbb{E}\left[Y(0)|S(1)=1,S(0)=1\right]$ is point identified. In contrast, the outcome in the treatment group can be either an always-selected’s outcome or a selection complier’s outcome.

Always selected share among the treated is

$$p_0=\mathbf{Pr}\left(S(1)=1,S(0)=1|S(1)=1\right)=\mathbf{Pr}\left(S(0)=1|S(1)=1\right)=\frac{\mathbf{Pr}\left(S=1|D=0\right)}{\mathbf{Pr}\left(S=1|D=1\right)}$$

In the best case, the always-selected comprise the top $p_0$ quantile of the treatment outcomes. Then the largest possible value of $\beta$ is

$${\beta }_U=\mathbb{E}\left[Y|Y\ge Q_{y|S=1,D=1}(1-p_0),D=1,S=1\right]-\mathbb{E}\left[Y|S=1,D=0\right]$$

The smallest possible one is

$${\beta }_L=\mathbb{E}\left[Y|Y\le Q_{y|S=1,D=1}{p}_0,D=0,S=1\right]-\mathbb{E}\left[Y|S=1,D=0\right]$$

this can be implemented conditional on covariates by constructing $p_0(x)$ within each $x$ stratum.

Regression Discontinuity Design

Setup

Treatment ($D$) changes discontinuously at some particular value $x_0$ in $x$ [and nothing else does], so

$$D_i=1_{x_i\ge x_0}$$

Standard identification assumptions violated by definition because although unconfoundedness holds trivially since we have $D_i=1_{x_i\ge c}$, this also means overlap is always violated. Need to invoke continuity to do causal inference.

Identified at $x=c$, i.e. ${\tau }_c={\mu }_{(1)}(c)-{\mu }_{(0)}(c)$ via

$${\tau }_c:=\mathbb{E}\left[Y_1-Y_0|X=c\right]=\underset{x\downarrow c}{\lim }\mathbb{E}\left[Y|X=c\right]-\underset{x\uparrow c}{\lim }\mathbb{E}\left[Y|X=c\right]$$ $$\underset{x\downarrow c}{\lim }\mathbb{E}\left[y|X\right]-\underset{x\uparrow c}{\lim }\mathbb{E}\left[y|X\right]={\tau }_{SRDD}+\left(\underset{x\downarrow c}{\lim }\mathbb{E}\left[u|X\right]-\underset{x\uparrow c}{\lim }\mathbb{E}\left[u|X\right]\right)\approx {\tau }_{SRDD}$$

• Conditional mean function $\mathbb{E}\left[u|X\right]$ is continuous at $c$

• Mean Treatment effect function $\mathbb{E}\left[{\tau }_i|X\right]$ is right continuous at $c$

Estimators

Normalise running variable $c:=x_0$. Then, the linear regression implementation is the following:

$$Y={\alpha }_l+\tau D+{\beta }_{l}f(X-c)+({\beta }_r-{\beta }_l)\times D\times g(X-c)+ϵ$$

where $f$ and $g$ are local or global polynomials. Since the design relies on identification at infinity (i.e. at the cutoff), choice of polynomial / functional form matters a lot.

Calonico, Cattaneo, Titiunik (2014) recommend local-linear regressions. Older literature relies on global higher-order polynomials, which often yields strange estimates.

$${\hat{\tau }}_c=\arg \underset{a,\tau ,{\beta }_{(0)},{\beta }_{(1)}}{\min }\sum _{i=1}^{n}K\left(\frac{\left|X_i-c\right|}{h_n}\right){\left(Y_i-a-\tau D_i-{\beta }_{(0)}(X_i-c{)}_{-}-{\beta }_{(1)}(X_i-c{)}_{+}\right)}^2$$

Where $K(\cdot )$ is a kernel function. Common choices are the window function $K(x)=1_{|x|\le 1}$ or the triangular kernel $K(x)=(1-\left|x\right|{)}_{+}$

Assumptions for Local Linear Estimator

Loosely, we need CEFs ${\mu }_{(w)}$ to be smooth. More precisely, we need ${\mu }_{(w)}(x)$ to be twice-differentiable with uniformly bounded second derivative.

$$\left|\frac{d^2}{d{x}^2}{\mu }_{(w)}(x)\right|\le B,\forall x\in \mathbb{R},w\in \left\{0,1\right\}$$

Taking a taylor expansion around $c$, we can write the CEFs as

$${\mu }_{(w)}(x)=a_{(w)}+{\beta }_{(w)}(x-c)+\frac{1}{2}{\rho }_{(w)}(x-c{)}^2,\left|{\rho }_{(w)}(x)\right|\le B$$

with ${\tau }_c=a_{(1)}-a_{(0)}$. The local linear regression with a window kernel can be solved in closed form

$${\hat{a}}_{(1)}=\sum _{c\le X_i\le c+h_n}{\gamma }_i{Y}_i$$ $${\gamma }_i=\frac{{\hat{\mathbb{E}}}_{(1)}[(X_i-c{)}^2]-{\hat{\mathbb{E}}}_{(1)}[X_i-c]\cdot (X_i-c)}{{\hat{\mathbb{E}}}_{(1)}[(X_i-c{)}^2]-{\hat{\mathbb{E}}}_{(1)}[X_i-c{]}^2}$$

where $\hat{\mathbb{E}}(\cdot )$ denote sample averages over the regression window. Then, the error term can be written as

$${\hat{a}}_{(1)}=a_{(1)}+\sum _{c\le X_i\le c+h_n}{\gamma }_i{\rho }_{(1)}(X_i-c)+\sum _{c\le X_i\le c+h_n}{\gamma }_i(Y_i-{\mu }_{(1)}(X_i))$$

Curvature bias is bounded by $B{h}_n^2$.

$${\hat{\tau }}_c={\tau }_c+\mathsf{O}\left(n^{-2/5}\right),\text{with }{h}_n\sim n^{-1/5}$$

This rate is a consequence of working with the 2nd derivative. In general, if we assume ${\mu }_{(w)}(\cdot )$ has a bounded $k-$th derivative, we can achieve an $n^{-k/(2k+1)}$ rate using local polynomial regression of order $k-1$ with a bandwidth scaling as $h_n\sim n^{-1/(2k+1)}$.

The local linear regression estimator for ${\tau }_c$

$${\hat{\tau }}_c=\arg \underset{a,\tau ,{\beta }_{(0)},{\beta }_{(1)}}{\min }\sum _{i=1}^{n}K\left(\frac{\left|Z_i-c\right|}{h_n}\right){\left(Y_i-a-\tau W_i-{\beta }_{(0)}(Z_i-c{)}_{-}-{\beta }_{(1)}(Z_i-c{)}_{+}\right)}^2$$

which can be written as a local linear estimator ${\hat{\tau }}_c={\sum }_{i=1}^n{\gamma }_i{Y}_i$ where weights ${\gamma }_i$ only depend on the running variable $Z$. show that local linear regression is not the best estimator in this class.

Under an assumption that $\left|{\mu }_{(w)}^{\prime \prime }(z)\right|\le B|\left\{Z_1,\dots ,Z_n\right\}$, the minimax linear estimator is the one that minimises the MSE $\text{MSE}({\hat{\tau }}_c(\gamma )|\left\{Z_1,\dots ,Z_n\right\})\le \sigma {\Vert \gamma \Vert }_2^2+I_B^2(\gamma )$ and is given by

$${\hat{\tau }}_c({\gamma }^B)=\sum _{i=1}^n{\gamma }_i^B{Y}_i,{\gamma }^B=\arg \underset{\gamma }{\min }\left\{\sigma {\Vert \gamma \Vert }_2^2+I_B^2(\gamma )\right\}$$

These weights can be solved for using quadratic programming.

Fuzzy RD

Discontinuity doesn’t deterministically change treatment, but affects probability of treatment. Analogue of IV with one-sided non-compliance.

$$P[D_i=1|x_i]=g_0(x_i)1_{x_i<x_0}+g_1(x_i)1_{x_i\ge x_0}$$

$g_0(x_i)\ne g_1(x_i)$. Assuming $g_1(x_0)>g_0(x_0)$, the probability of treatment relates to $x_i$ via:

$E[D_i|x_i]=P[D_i=1|x_i]=g_0(x_i)+[g_1(x_i)-g_0(x_i)]T_i$ where $T_i=1_{x_i\ge x_0}$ := point of discontinuity

Regression Kink Design

First-derivative version of the fuzzy RD. Continuous treatment, where the treatments are a function of the running variable $X$ with kink at $x_0$. This implies that the first derivative $\frac{\partial D}{\partial X}$ of continuous treatment D is discontinuous at the threshold.

The marginal treatment effect at the threshold is defined as

$${\Delta }_{X=x_0}(d_0)=\frac{\partial \mathbb{E}\left[Y(d_0)|X=x_0\right]}{\partial D}=\frac{{\lim }_{\epsilon \to 0}\frac{\partial \mathbb{E}\left[Y|X=x_0+\epsilon \right]}{\partial X}-{\lim }_{\epsilon \to 0}\frac{\partial \mathbb{E}\left[Y|X=x_0-\epsilon \right]}{\partial X}}{{\lim }_{\epsilon \to 0}\frac{\partial \mathbb{E}\left[D|X=x_0+\epsilon \right]}{\partial X}-{\lim }_{\epsilon \to 0}\frac{\partial \mathbb{E}\left[D|X=x_0-\epsilon \right]}{\partial X}}$$

Differences-in-Differences

DiD with 2 periods

Binary treatment $d\in \left\{0,1\right\}$, 2 time periods $t\in \left\{0,1\right\}$.

Potential outcomes denoted $Y_t^d$.

ATT in the $2$nd period.

$${\tau }_{ATT}:=\mathbb{E}\left[Y_1^1-Y_1^0|D=1\right]$$

$\mathbb{E}\left[Y_1^0|D=1\right]$ not observed, so must be imputed.

Naive Estimation Strategies

Before-After Comparison: $\tau = \Exp{Y_1^1 | D = 1}

• \hlblu{\Exp{Y_0^0 | D = 1}}$

assumes $\mathbb{E}\left[Y_1^0|D=1\right]=\mathbb{E}\left[Y_0^{0}D=1\right]$ (No trending)

Post Treatment-Control Comparison: $\tau =\mathbb{E}\left[Y_1^1|D=1\right]-\mathbb{E}\left[Y_1^0|D=0\right]$

Assumes $\mathbb{E}\left[Y_1^0|D=1\right]=\mathbb{E}\left[Y_1^0|D=0\right]$ (Random Assignment in the 2nd period)

Both typically untenable in practice, so we need parallel trends.

Sample analogue of

Impute $\mathbb{E}\left[Y_1^0|D=1\right]$ with

$$\mathbb{E}\left[Y_1^0|D=1\right]\approx \mathbb{E}\left[Y_0^0|D=1\right]+\left(\mathbb{E}\left[Y_1^0|D=0\right]-\mathbb{E}\left[Y_0^0|D=0\right]\right)$$ $${\Delta }_{D=1}:=\left(\mathbb{E}\left[Y_1^1|D=1\right]-\mathbb{E}\left[Y_0^0|D=1\right]\right)-\left(\mathbb{E}\left[Y_1^0|D=0\right]-\mathbb{E}\left[Y_0^0|D=0\right]\right)$$ $$\mathbb{E}\left[Y_1^0-Y_0^0|D=1\right]=\mathbb{E}\left[Y_1^0-Y_0^0|D=0\right]$$

Often justified using a figure [with transformed $y$ if necessary], or control for time trends [which relies on a strong functional form assumption], or a clear falsification test [on a placebo group].

If $\mathbb{E}\left[Y_0^0|D=1\right]=\mathbb{E}\left[Y_0^0|D=0\right]$, this collapses to a selection-on-observables assumption in period 2.

$$\mathbb{E}\left[Y_1^0|D=1\right]=\mathbb{E}\left[Y_1^0|D=0\right]$$

For a two-period difference, we can also write the standard OLS exogeneity condition in differences form

$$\mathbb{E}\left[\Delta {\mathbf{x}}^{\prime }\Delta \mathbf{ϵ}\right]=0$$ $$\mathbb{E}\left[{\mathbf{x}}_2^{\prime }{\mathbf{ϵ}}_2\right]+\mathbb{E}\left[{\mathbf{x}}_1^{\prime }{\mathbf{ϵ}}_1\right]-\left(\mathbb{E}\left[{\mathbf{x}}_1^{\prime }{\mathbf{ϵ}}_2\right]-\mathbb{E}\left[{\mathbf{x}}_2^{\prime }{\mathbf{ϵ}}_1\right]\right)=0$$

Which makes a direct link with the strong exogeneity assumption in panel data models that asserts that ${ϵ}_t\perp \perp {\mathbf{x}}_1,\dots {\mathbf{x}}_t$.

Regression Estimator

We typically prefer the following regression estimator (for automatic standard errors etc).

$$Y_{it}=\alpha +\gamma {\text{Treat}}_i+\lambda {\text{Post}}_t+\tau ({\text{Treat}}_i\times {\text{Post}}_t)+{\epsilon }_{it}$$

Triple Differences (DDD) Estimator

Regular Diff-in-Diff estimate - Diff-in-diff estimate for placebo group.

Nonparametric Identification Assumptions with Covariates

Estimand:

$${\tau }_{ATT}:=\mathbb{E}\left[Y_t^1-Y_t^0|D=1\right]=\mathbb{E}\left[{\theta }_t(\mathbf{x})|D=1\right]={\mathbb{E}}_{\mathbf{X}|D=1}\left[{\theta }_t(\mathbf{x})\right]$$

Identification Assumptions:

SUTVA $Y_t=D{Y}_t^1+(1-D)Y_t^0,t\in \left\{0,1\right\}$

Covariate exogeneity

${\mathbf{X}}^1={\mathbf{X}}^0=\mathbf{X}x\in \mathcal{X}$

No effect before treatment

${\theta }_0(\mathbf{x})=0;\forall x\in \mathcal{X}$

Common Trend

(parallel trends within $\mathbf{x}$ strata)

$$\mathbb{E}\left[Y_1^0|\mathbf{X}=\mathbf{x},D=1\right]-\mathbb{E}\left[Y_0^0|\mathbf{X}=\mathbf{x},D=1\right]=\mathbb{E}\left[Y_1^0|\mathbf{X}=\mathbf{x},D=0\right]-\mathbb{E}\left[Y_0^0|\mathbf{X}=\mathbf{x},D=0\right]=\mathbb{E}\left[Y_1^0|\mathbf{X}=\mathbf{x}\right]-\mathbb{E}\left[Y_0^0|\mathbf{X}=\mathbf{x}\right]$$

Common support

$$\mathbf{Pr}\left(T=1,D=1|\mathbf{X}=\mathbf{x},(T,D)\in \left\{(t,d),(1,1)\right\}\right)<1,\forall (t,d)\in \left\{(0,1),(0,0),(1,0)\right\},x\in \mathcal{X}$$

This allows us to estimate the conditional ATT as the standard DiD within each $X$ stratum.

$$\mathbb{E}\left[Y_1|D=1,X\right]-\mathbb{E}\left[Y_0|D=1,X\right]-\mathbb{E}\left[Y_1\mid D=0,X\right]+\mathbb{E}\left[Y_0\mid D=0,X\right]$$

Averaging these over $\mathsf{d}X$ gives us the ATT

$${\tau }_1^{\text{ATT}}=\mathbb{E}\left[\left\{{\mu }_1(1,X)-{\mu }_1(0,X)\right\}-\left\{{\mu }_0(1,X)-{\mu }_0(0,X)\right\}|D=1,T=1\right]$$

where regression functions ${\mu }_d(t,x)$ denote conditional expectations for treatment $d$ at time $t$ given covariates $x$.

Denote potential outcomes under treatment and control for unit $i$ as $Y_{it}^1$ and $Y_{it}^0$. For some observed covariates $X_i$, we are interested in the CATT

${\tau }_0(X_i):=\mathbb{E}\left[Y_{i1}^1-Y_{i1}^0|X_i,D_i=1\right]$

For identification, we need Conditional parallel trends: $\mathbb{E}\left[Y_{i1}^0-Y_{i0}^0|D_i=1,X_i\right]=\mathbb{E}\left[Y_{i1}^0-Y_{i0}^0|D_i=0,X_i\right]$ Overlap: $\exists c>0$ such that $\mathbb{E}\left[D_i=1|X_i\right]>c$ and $\mathbb{E}\left[D_i|X_i\right]<1-c$

The Abadie estimand can be defined as

$$\mathbb{E}\left[Y_{i1}^1-Y_{i1}^0|X_i,D_i\right]=\mathbb{E}\left[\frac{D_i-\mathbb{E}\left[D_i=1|X_i\right]}{\mathbb{E}\left[D_i=1|X_i\right](1-\mathbb{E}\left[D_i=1|X_i\right])}(Y_{i1}-Y_{i0})|X_i\right]$$

Defining $\Delta Y_i:=Y_{i1}-Y_{i0}$, we then have

$$\mathbb{E}\left[Y_{i1}^1-Y_{i1}^0|X_i,D_i\right]=\mathbb{E}\left[\frac{D_i-\mathbb{E}\left[D_i=1|X_i\right]}{\mathbb{E}\left[D_i=1|X_i\right]\mathbb{E}\left[D_i=0|X_i\right]}\Delta Y_i|X_i\right]$$ $$\mathbb{E}\left[Y_{i1}^1-Y_{i1}^0|X_i,D_i\right]=\mathbb{E}\left[\frac{D_i\Delta Y_i}{\mathbb{E}\left[D_i=1|X_i\right]}|X_i\right]-\mathbb{E}\left[\frac{(1-D_i)\Delta Y_i}{(1-\mathbb{E}\left[D_i=1|X_i\right])}|X_i\right]$$

This is an IPW Estimator.

Integrating this over $\mathsf{d}P(X|D=1)$ gives us the ATT

$$\mathbb{E}\left[Y_1^1-Y_0^0|D=1\right]=\mathbb{E}\left[\frac{Y_1-Y_0}{\mathbf{Pr}\left(D=1\right)}\cdot \frac{D-\mathbb{E}\left[D=1|X_i\right]}{1-\mathbb{E}\left[D=1|X_i\right]}\right]$$

THe full IPW estimator can be written

$${\Delta }_{D=1,T=1}=\mathbb{E}\left[Y\cdot \left\{\frac{DT}{\Pi }-\frac{D(1-T){\rho }_{1,1}(X)}{{\rho }_{1,0}(X)\Pi }-\left(\frac{(1-D)T{\rho }_{1,1}(X)}{{\rho }_{0,1}(X)\Pi }-\frac{(1-D)(1-T){\rho }_{1,1}(X)}{{\rho }_{0,0}(X)\Pi }\right)\right\}\right]$$

where $\Pi =\mathbf{Pr}\left(D=1,T=1\right)$ is the unconditional probability of being treated in the post-treatment period, and ${\rho }_{d,t}(X)=\mathbf{Pr}\left(D=dT=t|X\right)$ are conditional probabilities of specific treatment-group combinations.

Double-robust version - Zimmert (2020)

$${\Delta }_{D=1,T=1}=\mathbb{E}\left[\left\{\frac{DT}{\Pi }-\frac{D(1-T){\rho }_{1,1}(X)}{{\rho }_{1,0}(X)\Pi }-\left(\frac{(1-D)T{\rho }_{1,1}(X)}{{\rho }_{0,1}(X)\Pi }-\frac{(1-D)(1-T){\rho }_{1,1}(X)}{{\rho }_{0,0}(X)\Pi }\right)\right\}\times (Y-{\mu }_d(T,X))+\frac{DT}{\Pi }\left(({\mu }_1(1,X)-{\mu }_1(0,X))-({\mu }_0(1,X)-{\mu }_0(0,X))\right)\right]$$

Panel Data

Setup: We observe a sample of $i=1,\dots ,N$ cross-sectional units for $t=1,\dots ,T$ time periods $⟹\text{Data : }\{(y_{it},{\mathbf{x}}_{it}^{\prime }):t=1,\dots ,T{\}}_{t=1}^T$

One-way fixed effects and Random effects both use the form

$$y_{it}={\mathbf{x}}_{it}^{\prime }\beta +\underset{e_{it}}{\underbrace{{\theta }_i+{ϵ}_{it}}}$$

although they make different assumptions about the error.

Error assumptions for panel regressions

(1) FE: $\mathbb{E}\left[{ϵ}_{it}|{\mathbf{x}}_i,{\theta }_i\right]=0⟺{\theta }_i/\perp \perp {\mathbf{x}}_i$.

(2) RE: (1) and $\mathbb{E}\left[e_{it}|{\mathbf{x}}_i\right]=0$ [Absorb unobserved unit effect into error term, impose orthogonality it] $⟹{\theta }_i\perp \perp {\mathbf{x}}_i$. Equivalent to Pooled OLS with FGLS.

Fixed Effects Regression

Identification Assumption []

• Strict Exogeneity - errors are uncorrelated with lags and leads of x $\mathbb{E}\left[{ϵ}_{it}|x_i\right]=\mathbb{E}\left[{ϵ}_{it}|x_{i1},\cdots x_{iT}\right]=0⟺\mathbb{E}\left[{\mathbf{x}}_{is}^{\prime }{\mathbf{ϵ}}_{it}\right]=0\forall s,t=1,\dots T$

  Equivalent statement for $y_{it}$ is $\mathbb{E}\left[y_{it}|{\mathbf{x}}_{i1},\dots ,{\mathbf{x}}_{iT}\right]=\mathbb{E}\left[y_{it}|{\mathbf{x}}_{it}\right]={\mathbf{x}}_{it}^{\prime }\mathbf{\beta }$

  • Rules out feedback loops i.e. $x_{it}$ correlated with ${ϵ}_{i,t-1}$ because $X$s are set in response to prior error, e.g. Policing and crime.

• regressors vary over time for at least some $i$.

Setup[]

Define an individual fixed effect for individual $i$

$$A_i=1_{\text{the observation involves unit }i}$$

and define the same for each time period for panel data.

If $D_{it}$ is as good as randomly assigned conditional on $A_i$:

$E[Y_{0it}|A_i,X_{it},t,D_{it}]=E[Y_{0it}|A_i,X_{it},t]$

Then, assuming $A_i$ enter linearly,

$$E[Y_{0it}|A_i,X_{it},t,D_{it}]=\alpha +{\lambda }_t+A_i^{\prime }\gamma +X_{it}^{\prime }\beta$$

Assuming the causal effect of the treatment is additive and constant,

$E[Y_{1it}|A_i,X_{it},t]=E[Y_{0it}|A_i,X_{it},t]+\rho$

where $\rho$ is the causal effect of interest.

Then, we can write:

$$Y_{it}={\alpha }_i+{\lambda }_t+\rho D_{it}+X_{it}^{\prime }\beta +{ϵ}_{it}$$ $${ϵ}_{it}=Y_{0it}-E[Y_{0it}|A_i,X_{it},t]$$ $${\alpha }_i=\alpha +A_i\gamma$$

Restrictions

• Linear

• Additive functional form

• Variation in $D_{it}$, over time, for $i$, must be as good as random

Estimate the specification

${\ddot{y}}_i={\ddot{\mathbf{x}}}_i^{\prime }\beta +{ϵ}_i$

where ${\ddot{k}}_i=M_i{k}_i$ individual demeaned values from pre-multiplying by the Individual specific demeaning operator ${\mathbf{M}}_i:={\mathbf{I}}_i-1_i{\left(1_i^{\prime }{1}_i\right)}^{-1}{1}_i^{\prime }$ with every component in the general model above, which removes the fixed effect ${\theta }_i$.

Lag the general model one period and subtracting gives

$\Delta y_{it}=\Delta {\mathbf{x}}_{it}^{\prime }\mathbf{\beta }+\Delta {ϵ}_{it}$

where $\Delta y_{it}=y_{it}-y_{i,t-1}$ and so on. This naturally eliminates the time-invariant fixed effect ${\theta }_i$. The pooled OLS estimation of $\beta$ in the above regression is called the first differences (FD) estimator ${\hat{\mathbf{\beta }}}_{FD}$.

FE estimator is more efficient under the assumption that ${ϵ}_{it}$ are serially uncorrelated [$\mathbb{E}\left[{\mathbf{e}}_i{\mathbf{e}}_i^{\prime }|{\mathbf{x}}_i,{\theta }_i\right]={\sigma }_e^2{\mathbf{I}}_T$]

FD more efficient when ${ϵ}_{it}$ follows a random-walk.

For Individual Fixed Effects/Within estimation, using the regression anatomy formula, write:

$${\hat{\rho }}_{\text{FE}}=\frac{\mathrm{Cov}\left((\right){\ddot{Y}}_{it},{\ddot{D}}_{it})}{\mathrm{Var}\left((\right){\ddot{D}}_{it})}$$

Since $t=2$, ${\overline{Y}}_i=Y_{it}+\frac{\Delta Y_{it}}{2}$ and ${\overline{D}}_i=D_{it}+\frac{\Delta D_{it}}{2}$

$${\hat{\rho }}_{\text{FE}}=\frac{\mathrm{Cov}\left((\right){\ddot{Y}}_{it},{\ddot{D}}_{it})}{\mathrm{Var}\left((\right){\ddot{D}}_{it})}=\frac{\mathrm{Cov}\left((\right)\Delta Y_{it},\Delta D_{it})}{\mathrm{Var}\left((\right)\Delta D_{it})}={\hat{\rho }}_{\text{FD}}$$

Random Effects

Identification Assumption []

Assume ${\theta }_i\perp \perp X_i⟺\mathbb{E}\left[{\theta }_i|{\mathbf{x}}_i\right]=\mathbb{E}\left[{\theta }_i\right]=0$ - strong assumption

In other words, entire error term $e_{it}={\nu }_{it}+{\theta }_i$ is independent of $X$. This assumes OLS is consistent but inefficient, which is why it is of limited use in observational settings.

When there is autocorrelation in time series (i.e. ${ϵ}_t$ s are correlated over time ), GLS estimates can be obtained by estimating OLS on quasi-differenced data. This allows us to estimate the effects of time-invariant characteristics (assuming the independence condition is met). $y_{it}-\lambda \overline{y_i}=(x_{it}-\lambda \overline{x_i})\beta +(1-\lambda ){\theta }_i+{\nu }_{it}-\lambda {\bar{\nu }}_i$

where $\lambda =1-{\left[\frac{{\sigma }_{\nu }^2}{{\sigma }_{\nu }^2+T{\sigma }_{\theta }^2}\right]}^{\frac{1}{2}}$

Idiosyncratic errors ${\nu }_{it}$ have constant finite variance: $\mathbb{E}\left[{\nu }_{it}^2\right]={\sigma }_{\nu }^2$

Idiosyncratic errors ${\nu }_{it}$ are serially uncorrelated: $\mathbb{E}\left[{\nu }_{it}{\nu }_{is}\right]=0\forall t\ne s$.

$\mathbb{E}\left[{\theta }_i^2|{\mathbf{x}}_i\right]={\sigma }_{\theta }^2$

Under these assumptions, the FGLS matrix $\mathbf{\Omega }$ takes a special form

$\mathbf{\Omega }={\sigma }_{\nu }^2{\mathbf{I}}_T+{\sigma }_{\theta }^2{\mathbf{j}}_T{\mathbf{j}}_T^{\prime }$

where ${\mathbf{j}}_T{\mathbf{j}}_T^{\prime }$ is a $T\times T$ matrix of $1$s. Estimators for the variance components are in [c 10, pp 260-61]. A robust estimator of $\hat{\mathbf{\Omega }}$ is constructed using pooled OLS residuals $\widehat{{\mathbf{v}}_i}$

$\hat{\mathbf{\Omega }}=\frac{1}{n}{\sum }_{i=1}^n{\hat{\mathbf{v}}}_i{\hat{\mathbf{v}}}_i^{\prime }$

With this, we can apply the FGLS estimator

$${\hat{\mathbf{\beta }}}_{RE}={\left({\mathbf{X}}^{\prime }{\hat{\mathbf{\Omega }}}^{-1}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\prime }{\hat{\mathbf{\Omega }}}^{-1}\mathbf{y}$$

Hausman Test: Choosing between FE and RE

${\beta }_{FE}$ is assumed to be consistent. Oft-abused test as a result.

• H0: ${\beta }_{FE}-{\beta }_{RE}=0$

• H0: ${\beta }_{FE}-{\beta }_{RE}\ne 0$

$$H={\left({\hat{\mathbf{\beta }}}_{FE}-{\hat{\mathbf{\beta }}}_{RE}\right)}^{\prime }{\left[\widehat{\mathrm{Var}}\left({\hat{\mathbf{\beta }}}_{FE}\right)-\widehat{\mathrm{Var}}\left({\hat{\mathbf{\beta }}}_{RE}\right)\right]}^{-1}\left({\hat{\mathbf{\beta }}}_{FE}-{\hat{\mathbf{\beta }}}_{RE}\right)\stackrel{d}{\to }{\chi }_k^2$$

If the error component $\theta$ is correlated with $x$, RE estimates are not consistent. Perform Hausman test for random vs fixed effects (where under the null, $\mathrm{Cov}\left({\theta }_i,x_{it}\right)=0$)

• When the idiosyncratic error variance ${\hat{\sigma }}_{\nu }^2$ is large relative to $T_i{\hat{\sigma }}_{\theta }^2$, $\lambda \to 0$ and ${\hat{\beta }}_{RE}\approx {\hat{\beta }}_{pool}$. In words, the individual effect is relatively small, so Pooled OLS is suitable.

• When the idiosyncratic error variance ${\hat{\sigma }}_{\nu }^2$ is small relative to $T_i{\hat{\sigma }}_{\theta }^2$, $\lambda \to 1$ and ${\hat{\beta }}_{RE}\approx {\hat{\beta }}_{FE}$. Individual effects are relatively large, so FE is suitable.

Time Trends

Linear Time Trend[]

$$y_{it}={\mathbf{x}}_{it}\mathbf{\beta }+c_i+t+{\epsilon }_{it},t=1,2,\dots ,T$$

Time Fixed Effects (a.k.a. Two-way Fixed Effects)[]

$$y_{it}={\mathbf{x}}_{it}\mathbf{\beta }+c_i+t_t+{\epsilon }_{it},t=1,2,\dots ,T$$

Unit Specific Time Trends[]

$$y_{it}={\mathbf{x}}_{it}\mathbf{\beta }+c_i+g_i\cdot t+t_t+{\epsilon }_{it},t=1,2,\dots ,T$$

Distributed Lag

Define switching indicator $D_{it}$ as 1 if $i$ switched from control to treatment between $t-1$ and $t$.

$$Y_{ist}={\gamma }_s+{\lambda }_t+\sum _{\tau =0}^m{\delta }_{-\tau }{D}_{s,t-\tau }+\sum _{\tau =1}^q{\delta }_{+\tau }{D}_{s,t+\tau }+X_{ist}^{\prime }\beta +{ϵ}_{ist}$$

where the sums on the RHS allow for m lags / post-treatment effects, and q leads / pre-treatment effects. Leads should be close to 0.

Staggered Adoption

Let $T$ denote multiple time periods such that $t\in \left\{0,1,\dots ,\mathcal{T}\right\}$, with nobody treated at $t=0$ and staggered adoption. Let $G_t$ be a dummy that is equal to one if a subject experiences treatment introduction in period $t$ (e.g. $G_2=1$ implies the treatment is introduced in period $2$ in said group).

Under parallel trends for the untreated potential outcomes, $Y_{g,t}(0)$, the treatment effect ${\hat{\beta }}_{\text{FE}}$ in the vanilla two-way fixed effects regression

$$Y_{g,t}={\hat{\beta }}_{\text{FE}}{D}_{g,t}+{\gamma }_g+{\psi }_t+{\epsilon }_{gt}$$

can be decomposed as

$$\mathbb{E}\left[{\hat{\beta }}_{\text{FE}}\right]=\mathbb{E}\left[\sum _{(g,t):D_{g,t}\ne 0}{W}_{g,t}{\Delta }_{g,t}\right]$$

where ${\Delta }_{g,t}=Y_{g,t}(1)-Y_{g,t}(0)$.

The weights $W_{g,t}$ sum to one and are proportional to and the same sign as

$$N_{g,t}\left(D_{g,t}-D_{g,\cdot }-D_{\cdot ,t}+D_{\cdot ,\cdot }\right)$$

where $D_{g,\cdot }$ is the average treatment of group $g$ across periods (share of periods treated), $D_{\cdot ,t}$ is the average treatment at period $t$ across groups, and $D_{\cdot ,\cdot }$ is the grand mean of the treatment indicator. These weights can be negative.

This means that ${\hat{\beta }}_{\text{FE}}$ is biased for the ATT because $W_{g,t}$ is in not (only) proportional to $N_{g,t}$. $\hat{\beta }$ is only unbiased when

the treatment is binary AND the treatment is staggered and absorbing (i.e. groups get treated once and stay treated) AND there is no variation in treatment timing

Under these conditions, the pesky weight is constant across treated units, so the weights are proportional to $N_{g,t}$.

OR, ${\hat{\beta }}_{\text{FE}}$ is also unbiased if $(D_{g,t}-D_{g,\cdot }-D_{\cdot ,t}+D_{\cdot ,\cdot })$ is uncorrelated with the treatment effects ${\Delta }_{g,t}$. This is only plausible when treatment has been randomly staggered, otherwise, it is entirely plausible that groups with larger treatment effects selected into treatment early, and so on.

Consider a dataset comprising $K$ timing groups ordered by the time at which they first receive treatment and a maximum of one never-treated group $U$. The OLS estimate from a two-way fixed effects regression is

$${\hat{\beta }}_{DD}=\sum _{k\ne U}{s}_{kU}{\hat{\beta }}_{kU}^{DD}+\sum _{k\ne U}\sum _{j>k}\left(s_{kj}{\hat{\beta }}_{kj}^{DD}+s_{jk}{\hat{\beta }}_{jk}^{DD}\right)$$

where weights depend on sample size and variance of treatment within each DD. This maximises the weights of groups treated in the middle of the panel. The Late vs Early comparison is particularly problematic (and is typically incorrect when treatment effects are heterogeneous in time).

Visually, this involves decomposing the setup in the first figure below into its constituent two-way parts in the second figure.

Estimand: Group-time average treatment effect

$$ATT(g,t)=\mathbb{E}\left[Y_t(g)-Y_t(\infty )|G_g=1\right],\forall t\ge g$$

where $Y_t(g)$ is the potential outcome for group treated at $g$.

Separate (1) identification, (2) estimation and inference, and (3) aggregation.

A1: No anticipation $\forall i,t\text{ and }t<g,g^{\prime }$, $Y_{it}(g)=Y_{i,t}(g^{\prime })$

A2: Parallel trends based on ‘never treated’ group: $\forall t\in \left\{2,\dots ,\mathcal{T}\right\}$, $g\in \mathcal{G}$ s.t. $t\ge g$, $\underset{\text{Trend in group treated at 1}}{\underbrace{\mathbb{E}\left[Y_t(0)-Y_{t-1}(0)|G_g=1\right]}}=\underset{\text{Trend in never treated}}{\underbrace{\mathbb{E}\left[Y_t(0)-Y_{t-1}(0)|C=1\right]}}$

Estimators for Group-time ATEs

$$AT{T}_{\text{unc}}^{\text{never}}(g,t)=\mathbb{E}\left[Y_t-Y_{g-1}|G_g=1\right]-\mathbb{E}\left[Y_t-Y_{g-1}|C=1\right]$$ $$AT{T}_{\text{unc}}^{\text{notyet}}(g,t)=\mathbb{E}\left[Y_t-Y_{g-1}|G_g=1\right]-\mathbb{E}\left[Y_t-Y_{g-1}|D_t=0,C=1\right]$$

Aggregation: event-study type estimand.

$${\theta }_D(e)=\sum _{g=2}^{\mathcal{T}}{1}_{g+e\le \mathcal{T}}ATT(g,g+e)\mathbf{Pr}\left(G=g|G+e\le \mathcal{T},C\ne 1\right)$$

Implemented in did and DRDID.

The negative weighting problem with 2WFE under staggered adoption can be remedied easily by using the following procedure, which is termed Imputation by . This nests the procedures in etc.

Fit a model for $Y_0^{([)}][it]$ using only untreated observations for all units (i.e. untreated periods for units that eventually got treated) Impute $\widehat{Y_0^{([)}]}$ for treated units and treated time periods compute ${\hat{\tau }}_{it}=Y_{it}-{\hat{Y}}_{it}^{(0)}\mid \forall i,t\text{ where }{W}_{it}=1$ Average for (equal weighting) ATT or average over time for event study

This works well when the outcome model for $Y_0^{([)}][it]$ is good, i.e. when the fixed effects or latent factors are well estimated. This will not work well for short panels.

Changes-in-Changes

Given a continuous outcome $Y$ and a monotonicity in unobserved heterogeneity, CiC allows us to identify both the ATT and Quantile effect on the treated (QTT).

Assume the following about untreated potential outcomes

$Y_T^0=\mathcal{H}(U,T)U\perp \perp T|D$

where $U$ is a scalar unobservable or an index of unobservables. $\mathcal{H}(u,t)$ is a general function assumed to be strictly monotonically increasing in values of $u$ for periods $t\in \left\{0,1\right\}$. The conditional independence assumption requires that the unobserved heterogeneity is constant over time within treatment groups.

Denote ${\mathbb{F}}_{Y(d)|dt}(y)=\mathbf{Pr}Y(d)\le y|D=d,T=t$ the conditional CDF of potential outcome $Y(d)$, and ${\mathbb{F}}_{dt}(y)=\mathbf{Pr}D=d,T=t$ corresponding CDF for observed outcome. Conditional outcome distributions ${\mathbb{F}}_{01},{\mathbb{F}}_{00},{\mathbb{F}}_{10}$ are observed. The inverse of the latter is ${\mathbb{F}}_{dt}^{-1}(y)$, the conditional quantile function. The unobserved CDF is identified as

${\mathbb{F}}_{Y(0)|11}(y)={\mathbb{F}}_{10}\left({\mathbb{F}}_{00}^{-1}({\mathbb{F}}_{01}(y))\right)$

The QTT at quantile $\tau$ is then identified as

$${\Delta }_{D=1}(\tau )={\mathbb{F}}_{11}^{-1}(\tau )-{\mathbb{F}}_{01}^{-1}\left({\mathbb{F}}_{00}({\mathbb{F}}_{10}^{-1}(\tau ))\right)$$

and the ATT is identified as

$${\Delta }_{D=1}=\mathbb{E}\left[Y|D=1,T=1\right]-\mathbb{E}\left[{\mathbb{F}}_{01}^{-1}({\mathbb{F}}_{00}(Y_{10}))\right]$$

Implemented in qte::CiC.

Synthetic Control

Original setup.

Observe $n_0{+}_1$ units in periods $t=1,\dots ,T$. Unit 1 is treated starting from period $T_0+1$, while $2,\dots ,n_0+1$ are never treated, and are therefore called the donor pool.

$$Y_{it}^{obs}=Y_{it}(D_{it})=(1-D_{it})Y_{it}(0)+D_{it}{Y}_{it}(1)$$

Since there is only 1 treated unit, the effect of interest ${\tau }_t:=Y_{1i}(t)-Y_{0i}(t),t=T_0{+}_1,\dots ,T$

Observed data matrix ()

$${\mathbf{Y}}^{obs}:=(Y_{it}^{obs}{)}_{t=1,\dots ,T,i=1,\dots ,n_0+1}$$

FPCI applies; potential outcome matrices are:

$$\mathbf{Y}(0):=(Y_{it}(0){)}_{t=1,\dots ,T,i=1,\dots ,n_0+1}$$ $$\mathbf{Y}(1):=(Y_{it}(1){)}_{t=1,\dots ,T,i=1,\dots ,n_0+1}$$

Let ${\mathbf{X}}_{treat}$ be a $p-$vector of a pre-intervention characteristics, and ${\mathbf{X}}_c$ is a $p\times n_0$ matrix containing the same values for control units. This typically includes pre-treatment outcomes, in which case $p=T_0$, but predictors (even time invariant ones, $Z_i$) are usually available.

$${\mathbf{X}}_{treat}:=(Y_{1,1}^{obs},Y_{1,2}^{obs},\dots ,Y_{1,T_0}^{obs},{\mathbf{Z}}_i{)}^{\top }$$

For some $p\times p$ PSD matrix $\mathbf{V}$, define ${\left|\left|\mathbf{X}\right|\right|}_{\mathbf{V}}=\sqrt{{\mathbf{X}}^{\prime }\mathbf{VX}}$, where $\mathbf{V}$ is typically diagonal. Consider weights $\mathbf{\omega }=({\omega }_2,\dots ,{\omega }_{n_0+1})$ satisfying

$${\omega }_i\ge 0,i=2,\dots ,n_0+1$$ $$\sum _{i\ge 2}{\omega }_i=1$$

This forces interpolation, i.e. the counterfactual cannot take a value greater than the maximal value or smaller than the minimal value of for a control unit. The synthetic control solution ${\mathbf{\omega }}^{\ast }$ solves

$$\underset{\mathbf{\omega }}{\min }{\left|\left|{\mathbf{X}}_{treat}-{\mathbf{X}}_c\mathbf{\omega }\right|\right|}_{\mathbf{V}}^2\text{s.t. non-negativity and sum-to-1 constraints}$$

The Synthetic Control Estimator is then

${\hat{\tau }}_t:=Y_{1,t}^{obs}-{\sum }_{i=2}^{n_0+1}{\omega }_i^{\ast }{Y}_{it}^{obs}$

In contrast, a simple difference-in-differences estimator gives

$${\hat{\tau }}_t^{DID}:=Y_{1,t}^{obs}-\left(Y_{1,T_0}^{obs}-\frac{1}{n_0}\sum _{i=2}^{n_0+1}\left(Y_{it}^{obs}-Y_{i,T_0}^{obs}\right)\right)$$

choose $\mathbf{V}=\mathrm{diag}{v}_1,\dots ,v_p$ using a nested-minimisation of the Mean Square Prediction Error (MSPE) over the pre-treatment period

$$\text{MSPE}(\mathbf{V}):=\sum _{t=1}^{T_0}{\left(Y_{1,t}^{obs}-\sum _{i=2}^{n_0+1}{\omega }_i(\mathbf{V})Y_{it}^{obs}\right)}^2$$

Setup:

$${\mathbf{Y}}^{\text{obs}}=\left[\begin{array}{cc}{\mathbf{Y}}_{\text{t, post}}^{\text{obs}}& {\mathbf{Y}}_{\text{c, post}}^{\text{obs}}\\ {\mathbf{Y}}_{\text{t, pre}}^{\text{obs}}& {\mathbf{Y}}_{\text{c, pre}}^{\text{obs}}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{Y}}_{\text{t, post}}(1)& {\mathbf{Y}}_{\text{c, post}}(0)\\ {\mathbf{Y}}_{\text{t, pre}}(0)& {\mathbf{Y}}_{\text{c, pre}}(0)\end{array}\right]$$ $$\mathbf{Y}(0)=\left[\begin{array}{cc}?& {\mathbf{Y}}_{\text{c, post}}(0)\\ {\mathbf{Y}}_{\text{t, pre}}(0)& {\mathbf{Y}}_{\text{c, pre}}(0)\end{array}\right]$$

relative magnitudes of $T$ and $N$ might dictate whether we impute the missing potential outcome $?$ using or comparison

Many Units and Multiple Periods: $N>>T_0$, $\mathbf{Y}(0)$ is ‘fat’, and comparison becomes challenging relative to . So matching methods are attractive.

$T_0>>N$, $\mathbf{Y}(0)$ is ‘tall’, and matching becomes infeasible. So it might be easier to estimate dependence structure.

Finally, if $T_0\approx N$, regularization strategy for limiting the number of control units that enter into the estimation of $Y_{0,T_0+1}(0)$ may be important

Focus on last period for now: ${\tau }_{0,T}=Y_{0,T}(1)-Y_{0,T}(0)=Y_{0,T}^{\text{obs}}-Y_{0,T}(0)$

Many estimators impute $Y_{0,T}(0)$ with the linear structure ${\hat{Y}}_{0,T}(0)=\mu +{\sum }_{i=1}^n{\omega }_i\cdot Y_{i,T}^{\text{obs}}$

Methods differ in how $\mu$ and $\omega$ are chosen as a function of ${\mathbf{Y}}_{\text{c, post}}^{\text{obs}},{\mathbf{Y}}_{\text{t, pre}}^{\text{obs}},{\mathbf{Y}}_{\text{c, pre}}^{\text{obs}}$

Impose four constraints

No Intercept: $\mu =0$. Stronger than Parallel trends in DiD.

Adding up : ${\sum }_{i=1}^n{\omega }_i=1$. Common to DiD, SC.

Non-negativity: ${\omega }_i\ge 0\forall i$. Ensures uniqueness via ‘coarse’ regularisation + precision control. Negative weights may improve out-of-sample prediction.

Constant Weights: ${\omega }_i=\overline{\omega }\forall i$

DiD imposes 2-4.

ADH(2010, 2014) impose 1-3

1 + 2 imply ‘No Extrapolation’.

Relaxing these assumptions:

Negative weights

If treated units are outliers on important covariates, negative weights might improve fit

Bias reduction - negative weights increase bias-reduction rate

When $N>>T_0$, (1-3) alone might not result in a unique solution. Choose by

Matching on pre-treatment outcomes : one good control unit is better than synthetic one comprised of disparate units

Constant weights - implicit in DiD

Given many pairs of $(\mu ,\omega )$

prefer values s.t. synthetic control unit is similar to treated units in terms of lagged outcomes

low dispersion of weights

few control units with non-zero weights

Optimisation Problem

Ingredients of objective function Balance: difference between pre-treatment outcomes for treated and linear-combination of pre-treatment outcomes for control ${\Vert {\mathbf{Y}}_{\text{t, pre}}-\mu -{\omega }^{\top }{\mathbf{Y}}_{\text{c, pre}}\Vert }_2^2=({\mathbf{Y}}_{\text{t, pre}}-\mu -{\omega }^{\top }{\mathbf{Y}}_{\text{c, pre}}{)}^{\top }({\mathbf{Y}}_{\text{t, pre}}-\mu -{\omega }^{\top }{\mathbf{Y}}_{\text{c, pre}})$ Sparse and small weights: sparsity : ${\Vert \omega \Vert }_1$ magnitude: ${\Vert \omega \Vert }_2$

$$({\hat{\mu }}^{en}(\lambda ,\alpha ),{\hat{\omega }}^{en}(\lambda ,\alpha ))=\underset{\mu ,\omega }{\mathrm{argmin}}Q(\mu ,\omega \mid {\mathbf{Y}}_{\text{t, pre}},{\mathbf{Y}}_{\text{c, pre}};\lambda ,\alpha ).$$ $$Q(\mu ,\omega \mid {\mathbf{Y}}_{\text{t, pre}},{\mathbf{Y}}_{\text{c, pre}};\lambda ,\alpha )={\Vert {\mathbf{Y}}_{\text{t, pre}}-\mu -{\omega }^{\top }{\mathbf{Y}}_{\text{c, pre}}\Vert }_2^2+\lambda \left(\frac{1-\alpha }{2}{\Vert \omega \Vert }_2^2+\alpha {\Vert \omega \Vert }_1\right).$$

Tailored Regularisation

don’t want to scale covariates ${\mathbf{Y}}_{\text{c, pre}}$ to preserve interpretability of weights. Instead, treat each control unit as a ‘pseudo-treated’ unit and compute

$${\hat{Y}}_{j,T}(0)={\hat{\mu }}^{\text{en}}(j;\alpha ,\lambda )+\sum _{i\ne j}{\hat{\omega }}_i(j;\alpha ,\lambda )\cdot Y_{i,T}^{\text{obs}}.$$

where

$$({\hat{\mu }}^{en}(j;\lambda ,\alpha ),{\hat{\omega }}^{en}(j;\lambda ,\alpha ))=\underset{\mu ,\omega }{\mathrm{argmin}}\sum _{t=1}^{T_0}{\left(Y_{j,t}-\mu -\sum _{i\ne 0,j}{\omega }_i{Y}_{i,t}\right)}^2+\lambda \left(\frac{1-\alpha }{2}{\Vert \omega \Vert }_2^2+\alpha {\Vert \omega \Vert }_1\right).$$

pick the value of the tuning parameters $({\alpha }_{opt}^{en},{\lambda }_{opt}^{en})$ that minimises

$$C{V}^{en}(\alpha ,\lambda )=\frac{1}{N}\sum _{j=1}^N\left(Y_{j,T}-{\hat{\mu }}^{en}(j;\alpha ,\lambda )-\sum _{i\ne 0,j}{\hat{\omega }}_i^{en}(j;\alpha ,\lambda )\cdot Y_{i,T}\right).$$

Difference in Differences

assume (2-4) No unique $\mu ,\omega$ solution for $T=2$, so fix $\omega =\frac{1}{N}$

$${\omega }_i^{\text{did}}=\frac{1}{N}\forall i\in \left\{1,\dots N\right\}.$$ $${\hat{\mu }}^{\text{did}}=\frac{1}{T_0}\sum _{s=1}^{T_0}{Y}_{0,s}-\frac{1}{N{T}_0}\sum _{s=1}^{T_0}\sum _{i=1}^N{Y}_{i,s}.$$

Best Subset; One-to-one Matching

$({\hat{\mu }}^S,{\hat{\omega }}^S)={\mathrm{argmin}}_{\mu ,\omega }Q(\cdot ;\lambda =0,\alpha )$ with ${\sum }_{i=1}^N{1}_{{\omega }_i\ne 0}\le k$ (=1 for OtO)

Synthetic Control

assume (1-3) (i.e. $\mu =0$) For $M\times M$ PSD diagonal matrix $\mathbf{V}$

$$(\hat{\omega }(\mathbf{V}),\hat{\mu }(\mathbf{V}))=\underset{\omega ,\mu }{\mathrm{argmin}}\left(({\mathbf{X}}_t-\mu -{\omega }^{\top }\mathbf{X}{)}^{\top }\mathbf{V}({\mathbf{X}}_t-\mu -{\omega }^{\top }\mathbf{X})\right).$$ $$\hat{\mathbf{V}}=\underset{\mathbf{V}=\mathrm{diag}(v_1,\dots ,v_M)}{\mathrm{argmin}}\left(({\mathbf{Y}}_{\text{t, pre}}-\hat{\omega }(\mathbf{V}{)}^{\top }{\mathbf{Y}}_{\text{c, pre}}{)}^{\top }({\mathbf{Y}}_{\text{t, pre}}-\hat{\omega }(\mathbf{V}{)}^{\top }{\mathbf{Y}}_{\text{c, pre}})\right).$$

Constrained regression: When $X_i=Y_{i,t};1\le t\le T_0$ (Lagged Outcomes only) $\mathbf{V}={\mathbf{I}}_N$ and $\lambda =0$

Consider a balanced panel with $N$ units and $T$ time periods, where the first $N_{co}$ units are never treated, while $N_{tr}=N-N_{co}$ treated units are exposed after time $T_{pre}$. We seek to solve for sdid weights ${\hat{\omega }}^{\text{sdid}}$ that align pre-exposure trends in outcomes of unexposed units with those for exposed units

$$\sum _{i=1}^{N_{co}}{\hat{\omega }}^{sdid}{Y}_{it}\approx N_{tr}^{-1}\sum _{i=N_{co}+1}^N{Y}_{it}$$

we also look for time weights ${\hat{\lambda }}_t^{sdid}$ that balance pre-exposure time periods with post-exposure time periods for unexposed units.

Weights are solved using the following optimisation problems

$$({\hat{\omega }}_0,{\hat{\omega }}^{\text{sdid}})=\underset{{\omega }_0\in \mathbb{R},\omega \in \Omega }{\mathrm{argmin}}{\ell }_{\text{unit}}({\omega }_0,\omega ).$$ $${\ell }_{\text{unit}}({\omega }_0,\omega )=\sum _{t=1}^{T_{\text{pre}}}{\left({\omega }_0+\sum _{i=1}^{N_{co}}{\omega }_i{Y}_{it}-\frac{1}{N_{\text{tr}}}\sum _{i=N_{\text{co}}+1}^N{Y}_{it}\right)}^2+{\zeta }^2{T}_{\text{pre}}{\Vert \omega \Vert }_2^2.$$ $$\Omega =\left\{\omega \in {\mathbb{R}}_{+}^N:\sum _{i=1}^{N_{\text{co}}}{\omega }_i=1,{\omega }_i=N_{\text{tr}}^{-1}\text{ for all }i=N_{\text{co}}+1,\dots ,N\right\}.$$

where ${\mathbb{R}}_{+}$ denotes the positive real line. We set the regularization parameter $\zeta$ as

$$\zeta ={\left(N_{\mathrm{tr}}{T}_{\mathrm{post}}\right)}^{1/4}\hat{\sigma },{\hat{\sigma }}^2=\frac{1}{N_{\mathrm{co}}(T_{\mathrm{pre}}-1)}\sum _{i=1}^{N_{\mathrm{co}}}\sum _{t=1}^{T_{\mathrm{pre}}-1}{\left({\Delta }_{it}-\bar{\Delta }\right)}^2.$$ $${\Delta }_{it}=Y_{i(t+1)}-Y_{it},\bar{\Delta }=\frac{1}{N_{\mathrm{co}}(T_{\mathrm{pre}}-1)}\sum _{i=1}^{N_{\mathrm{co}}}\sum _{t=1}^{T_{\mathrm{pre}}-1}{\Delta }_{it}.$$

We implement this for the time weights ${\hat{\lambda }}^{\text{sdid }}$ by solving

$$({\hat{\lambda }}_0,{\hat{\lambda }}^{\text{sdid}})=\underset{{\lambda }_0\in \mathbb{R},\lambda \in \Lambda }{\mathrm{argmin}}{\ell }_{\text{time}}({\lambda }_0,\lambda ).$$ $${\ell }_{\text{time}}({\lambda }_0,\lambda )=\sum _{i=1}^{N_{\mathrm{co}}}{\left({\lambda }_0+\sum _{t=1}^{T_{\mathrm{pre}}}{\lambda }_t{Y}_{it}-\frac{1}{T_{\mathrm{post}}}\sum _{t=T_{\mathrm{pre}}+1}^T{Y}_{it}\right)}^2.$$ $$\Lambda =\left\{\lambda \in {\mathbb{R}}_{+}^T:\sum _{t=1}^{T_{\mathrm{pre}}}{\lambda }_t=1,{\lambda }_t=T_{\text{post}}^{-1}\text{ for all }t=T_{\mathrm{pre}}+1,\dots ,T\right\}.$$

Compute regularisation parameter $\zeta$. Compute unit weights ${\hat{\omega }}^{\text{sdid}}$. Compute time weights ${\hat{\lambda }}^{\text{sdid}}$. Compute the SDID estimator using the following weighted DID regression.

$$({\hat{\tau }}^{\text{sdid}},\hat{\mu },\hat{\alpha },\hat{\beta })=\underset{\tau ,\mu ,\alpha ,\beta }{\mathrm{argmin}}\left\{\sum _{i=1}^N\sum _{t=1}^T{\left(Y_{it}-\mu -{\alpha }_i-{\beta }_t-D_{it}\tau \right)}^2{\hat{\omega }}_i^{sdid}{\hat{\lambda }}_t^{sdid}\right\}.$$

implemented in synthdid::synthdid_estimate

$$Y_{it}={\delta }_{it}{D}_{it}+{\mathbf{x}}_{it}^{\prime }\mathbf{\beta }+{\mathbf{\lambda }}_i^{\prime }{\mathbf{f}}_t+{\epsilon }_{it}.$$

Where $D$ is the treatment, ${\delta }_{it}$ is the heterogeneous treatment effect for unit $i$ at time $t$, ${\mathbf{x}}_{it}$ is a $p-$ vector of time-varying controls. ${\mathbf{f}}_t=[f_{1t},\dots ,f_{rt}{]}^{\prime }$ is a $k\times 1$ vector of unknown common factors, ${\mathbf{\lambda }}_i=[{\lambda }_{i1},\dots ,{\lambda }_{ir}{]}^{\prime }$ is a $r\times 1$ vector of unknown factor loadings. This factor component nests standard functional forms

$$\underset{\text{Confounders}}{\underbrace{U_{it}}}=\underset{\text{Loadings}}{\underbrace{{\lambda }_i}}\times \underset{\text{factors}}{\underbrace{f_t}}.$$

$f_t=1⟹{\lambda }_i\times 1={\lambda }_i$ unit FEs

${\lambda }_i=1⟹1\times f_t=f_t$ time FEs

$f_{1t}=1,f_{2t}={\xi }_t,{\lambda }_{i1}={\alpha }_i,{\lambda }_{i2}=1⟹f_t\times {\lambda }_i={\alpha }_i+{\xi }_t$ two-way FEs.

$f_t=t⟹{\lambda }_i\times f_t={\lambda }_i\times t$ Unit-specific linear time trends

${\lambda }_i=y_{i0},f_t={\alpha }_t⟹{\lambda }_i\times f_t=\alpha y_{i,t-1}-{\nu }_{it}$ Lagged dependent variable

Steps

Get initial value of $\hat{\beta }$ using within estimator

Estimate ${\hat{\lambda }}_i,{\hat{f}}_t$ using $\hat{\beta }$

Re-estimate $\hat{\beta }$ using ${\widehat{{\lambda }_i}}^{\prime }{\hat{f}}_t$

Iterate

Drawback - constant effect

With, $N_{CO}$ control units and $N_{TR}$ treated units, Write DGP for individual unit as

$${\mathbf{Y}}_i={\mathbf{D}}_i\circ {\mathbf{\delta }}_i+{\mathbf{X}}_i^{\prime }\mathbf{\beta }+\mathbf{F}{\mathbf{\lambda }}_i+{\mathbf{\epsilon }}_i,i\in \left\{1,2,\dots ,N\right\}.$$

Where ${\mathbf{Y}}_i=[y_{i1},y_{i2},\dots ,y_{iT}{]}^{\prime }$, ${\mathbf{D}}_i=[D_{i1},\dots ,D_{iT}{]}^{\prime }$, ${\mathbf{X}}_i=[{\mathbf{x}}_{i1},\dots ,{\mathbf{x}}_{iT}{]}^{\prime }$ is $T\times k$, $\mathbf{F}=[f_1,\dots ,f_T{]}^{\prime }$ is $T\times r$.

Stack controls together gives

$$\underset{T\times N_{CO}}{\underbrace{{\mathbf{Y}}_{CO}}}=\underset{T\times N_{CO}\times p}{\underbrace{{\mathbf{X}}_{CO}}}\underset{p\times 1}{\underbrace{\mathbf{\beta }}}+\underset{T\times N_{CO}}{\underbrace{\mathbf{F}}}\underset{\prime }{\underbrace{{\mathbf{\Lambda }}_{CO}}}+{\mathbf{\epsilon }}_{CO}.$$

GSC for treatment effects is an out-of-sample prediction method: the treatment effect for unit $i$ at time $t$ is the difference between the actual outcome and its estimated counterfactual ${\hat{\delta }}_{it}=Y_{it}(1)-{\hat{Y}}_{it}(0)$, where ${\hat{Y}}_{it}(0)$ is imputed in three steps.

Estimate an IFE model using only the control group data and estimate $\hat{\mathbf{\beta }},\hat{\mathbf{F}},{\hat{\mathbf{\Lambda }}}_{CO}$

$$(\hat{\mathbf{\beta }},\hat{\mathbf{F}},{\hat{\mathbf{\Lambda }}}_{CO})=\underset{\tilde{\mathbf{\beta }},\tilde{\mathbf{F}},\tilde{\mathbf{\Lambda }}}{\mathrm{argmin}}\sum _{i\in \mathcal{C}}{\left({\mathbf{Y}}_i-{\mathbf{X}}_i\tilde{\mathbf{\beta }}-\tilde{\mathbf{F}}{\tilde{\mathbf{\Lambda }}}_i\right)}^{\prime }\left({\mathbf{Y}}_i-{\mathbf{X}}_i\tilde{\mathbf{\beta }}-\tilde{\mathbf{F}}{\tilde{\mathbf{\Lambda }}}_i\right).$$ $${\tilde{\mathbf{F}}}^{\prime }\tilde{\mathbf{F}}={\mathbf{I}}_r,{\tilde{\mathbf{\Lambda }}}_{CO}^{\prime }{\tilde{\mathbf{\Lambda }}}_{CO}=\text{Diagonal}.$$

Estimate Factor loadings for each treated unit by minimising mean-squared error of the predicted treated outcome in pretreatment periods

$${\hat{\mathbf{\Lambda }}}_i=\underset{{\tilde{\mathbf{\Lambda }}}_i}{\mathrm{argmin}}{\left({\mathbf{Y}}_i^0-{\mathbf{X}}_i^0\hat{\mathbf{\beta }}-{\hat{\mathbf{F}}}^0{\tilde{\mathbf{\Lambda }}}_i\right)}^{\prime }\left({\mathbf{Y}}_i^0-{\mathbf{X}}_i^0\hat{\mathbf{\beta }}-{\hat{\mathbf{F}}}^0{\tilde{\mathbf{\Lambda }}}_i\right).$$ $${\hat{\mathbf{\Lambda }}}_i={\left({\hat{\mathbf{F}}}^{0^{\prime }}{\hat{\mathbf{F}}}^0\right)}^{-1}{\hat{\mathbf{F}}}^{0^{\prime }}\left({\mathbf{Y}}_i^0-{\mathbf{X}}_i^0\hat{\mathbf{\beta }}\right),i\in \mathcal{T}.$$

where $0$ superscripts denote the pretreatment periods.

Calculate Treated Counterfactuals based on $\hat{\mathbf{\beta }},\hat{\mathbf{F}},{\hat{\mathbf{\Lambda }}}_i$

$${\hat{Y}}_{it}(0)={\mathbf{x}}_{it}^{\prime }\hat{\mathbf{\beta }}+{\hat{\mathbf{\lambda }}}_i^{\prime }{\hat{\mathbf{f}}}_t,i\in \mathcal{T},t>T_0.$$

Choose the number of factors $r$ by cross-validation. Implemented in gsynth.

Dynamic Treatment Effects

We may want to estimate the effects of treatment sequences (‘time-varying exposures’), as in medical settings (Robins 1986, ).

2 period example

Consider a setting with $t=1,2$ and corresponding outcomes $Y_t$ and treatments $D_t$, where the treatment takes on values $d_1,d_2\in \left\{0,1,\dots ,J\right\}$, and baseline covariates ${\mathbf{X}}_0$ and covariates at the end of the first period ${\mathbf{X}}_1$.

Let ${\mathbf{d}}_2:=(d_1,d_2)\in \left\{0,1,\dots ,J\right\}\times \left\{0,1,\dots ,J\right\}$. Accordingly, $Y_2({\mathbf{d}}_2)$ is the potential outcome realised when treatment is set to sequence ${\mathbf{d}}_2$. The ATE (contrast) two distinct treatment sequences ${\mathbf{d}}_2$ vs ${\mathbf{d}}_2^{\prime }$ is

$$\Delta ({\mathbf{d}}_2,{\mathbf{d}}_2^{\prime }):=\mathbb{E}\left[Y_2({\mathbf{d}}_2)-Y_2({\mathbf{d}}_2^{\prime })\right].$$

Estimating this quantity requires a sequential selection on observables assumption

$$Y_2({\mathbf{d}}_2)\perp \perp D_1\mid {\mathbf{X}}_0,Y_2({\mathbf{d}}_2)\perp \perp D_2\mid D_1,{\mathbf{X}}_0,{\mathbf{X}}_1,d_1,d_2\in \left\{0,1,\dots ,J\right\}.$$ $$\mathbf{Pr}\left(D_1=d_1\mid {\mathbf{X}}_0\right)>0,\mathbf{Pr}\left(D_2=d_2\mid D_1,{\mathbf{X}}_0,{\mathbf{X}}_1\right)>0.$$

Under these assumptions, dynamic treatment effects can be estimated based on nested conditional means regressions

$${\Delta }^{\text{snmm}}({\mathbf{d}}_2,{\mathbf{d}}_2^{\prime }):=\mathbb{E}\left[\mathbb{E}\left[\mathbb{E}\left[Y_2\mid {\mathbf{D}}_2={\mathbf{d}}_2,{\mathbf{X}}_0,{\mathbf{X}}_1\right]\mid D_1=d_1,{\mathbf{X}}_0\right]-\mathbb{E}\left[\mathbb{E}\left[Y_2\mid {\mathbf{D}}_2={\mathbf{d}}_2^{\prime },{\mathbf{X}}_0,{\mathbf{X}}_1\right]\mid D_1=d_1^{\prime },{\mathbf{X}}_0\right]\right].$$

where ${\mathbf{d}}_2=(d_1,d_2)$ and ${\mathbf{d}}_2^{\prime }=(d_1^{\prime },d_2^{\prime })$ denote distinct treatment sequences.

or an IPW estimator

$${\Delta }^{\text{ipw}}({\mathbf{d}}_2,{\mathbf{d}}_2^{\prime }):=\mathbb{E}\left[\frac{Y\cdot 1_{D_1=d_1}{1}_{D_2=d_2}}{p^{d_1}({\mathbf{X}}_0)p^{d_2}(D_1,{\mathbf{X}}_0,{\mathbf{X}}_1)}-\frac{Y\cdot 1_{D_1=d_1^{\prime }}{1}_{D_2=d_2^{\prime }}}{p^{d_1^{\prime }}({\mathbf{X}}_0)p^{d_2^{\prime }}(D_1,{\mathbf{X}}_0,{\mathbf{X}}_1)}\right].$$

where $p^{d_1}({\mathbf{X}}_0)$ and $p^{d_2}$ are propensity scores in the two periods.

Finally, a double robust estimator is

$${\Delta }^{\text{dr}}({\mathbf{d}}_2,{\mathbf{d}}_2^{\prime })=\mathbb{E}\left[{\psi }^{{\mathbf{d}}_2}-{\psi }^{{\mathbf{d}}_2^{\prime }}\right].$$ $${\psi }^{{\mathbf{d}}_2}=\frac{1_{D_1=d_1}\cdot 1_{D_2=d_2}\cdot \left(Y_2-{\mu }^{Y_2}({\mathbf{d}}_2,{\mathbf{X}}_1)\right)}{p^{d_1}({\mathbf{X}}_0)p^{d_2}(D_1,{\mathbf{X}}_0,{\mathbf{X}}_1)}+\frac{1_{D_1=d_1}\cdot \left({\mu }^{Y_2}({\mathbf{d}}_2,{\mathbf{X}}_1)-{\nu }^{Y_2}({\mathbf{d}}_2,{\mathbf{X}}_0)\right)}{p^{d_1}({\mathbf{X}}_0)}+{\nu }^{Y_2}({\mathbf{d}}_2,{\mathbf{X}}_0).$$

where

$${\mu }^{Y_2}({\mathbf{d}}_2,{\mathbf{X}}_0,{\mathbf{X}}_1)=\mathbb{E}\left[Y_2\mid {\mathbf{D}}_2={\mathbf{d}}_2,{\mathbf{X}}_0,{\mathbf{X}}_1\right].$$ $${\nu }^{Y_2}({\mathbf{d}}_2,{\mathbf{X}}_0)=\mathbb{E}\left[\mathbb{E}\left[Y_2\mid {\mathbf{D}}_2={\mathbf{d}}_2,{\mathbf{X}}_0,{\mathbf{X}}_1\right]\mid D_1=d_1,{\mathbf{X}}_0\right].$$

are (nested) conditional mean outcomes.

If we assume that $D_2$ is conditionally independent of potential outcomes given pre-treatment covariates ${\mathbf{X}}_0$ and $D_1$ (implying that post-treatment ${\mathbf{X}}_1$ aren’t required to control for confounders jointly affecting the second treatment and the outcome). In this case, the second part of the first SOO assumption can be strengthened to $Y({\mathbf{d}}_2)\perp \perp D_2|D_1,{\mathbf{X}}_1$. This simplifies

$${\psi }^{{\mathbf{d}}_2}=\frac{1_{D_1=d_1}\cdot 1_{D_2=d_2}\cdot \left(Y_2-{\mu }^{Y_2}({\mathbf{d}}_2,{\mathbf{X}}_0)\right)}{p^{d_1}({\mathbf{X}}_0)p^{d_2}(d_1,{\mathbf{X}}_0)}+{\mu }^{Y_2}({\mathbf{d}}_2,{\mathbf{X}}_0).$$

implemented in causalweight::dyntreatDML.

Generalisation to arbitrary panels

Let $D_{it}$ denote treatment status at time $t$, and collect them into a $t-$ vector for each unit to form a Treatment History ${\mathbf{D}}_i:=(D_{i1},D_{i2},\dots ,D_{iT})$. A partial treatment history up to time $t$ is denoted ${\mathbf{D}}_{i,1:t}$. Time varying covariates are arranged analogously $X_{it},{\mathbf{X}}_{it},{\mathbf{X}}_{i,1:t}$.

Potential outcomes are defined on treatment histories and rely on the standard consistency assumption / SUTVA, which assumes that the potential outcome for the same observed history $Y_{it}:=Y_{it}({\mathbf{d}}_{1:t})$ when ${\mathbf{D}}_{i,1:t}={\mathbf{d}}_{1:t}$. This generates $2^t$ potential outcomes for the outcome in period $t$, which permits many hypothetical comparisons.

The estimand typically of interest the average causal effect of a treatment history

$$\tau ({\mathbf{d}}_{1:t},{\mathbf{d}}_{1:t}^{\prime }):=\mathbb{E}\left[Y_{it}({\mathbf{d}}_{1:t})-Y_{it}({\mathbf{d}}_{1:t}^{\prime })\right].$$

Define potential outcomes just intervening on the last $j$ periods as $Y_{it}({\mathbf{D}}_{i,1:t-j-1},{\mathbf{d}}_{t-j:t})$, which is the ‘marginal’ potential outcome if the treatment history runs its natural course up to $t-j-1$ and sets the last $j$ lags to ${\mathbf{d}}_{t-j:t}$.

This allows us to define a contemporaneous treatment effect (CET)

$${\tau }_c(t)=\mathbb{E}\left[Y_{it}({\mathbf{D}}_{i,1:t-1},1)-Y_{it}({\mathbf{D}}_{i,1:t-1},0)\right]=\mathbb{E}\left[Y_{it}(1)-Y_{it}(0)\right].$$

The $j-$step lagged effect is defined analogously

$${\tau }_l(t,j):=\mathbb{E}\left[Y_{it}({\mathbf{D}}_{i,1:t-j-1},1,0_j)-Y_{it}({\mathbf{D}}_{i,1:t-j-1},0,0_j)\right].$$

and the step response function (SRF) describes how this effect varies by time period and distance between the shift and the outcome

$${\tau }_s(t,j)=\mathbb{E}\left[Y_{it}(1_j)-Y_{it}(0_j)\right].$$

These effects are (clunkily) parametrised in an autoregressive distributed-lag (ADL) models of the form

$$Y_{it}={\beta }_0+\alpha Y_{i,t-1}+{\beta }_1{D}_{it}+{\beta }_2{D}_{i,t-1}+{\epsilon }_{it}.$$

with assumption ${\epsilon }_{it}\perp \perp {\mathbf{D}}_{i,s}\forall t,s$. This implies the following form for potential outcomes

$$Y_{it}({\mathbf{d}}_{1:t})={\beta }_0+\alpha Y_{i,t-1}({\mathbf{d}}_{1:t-1})+{\beta }_1{d}_{it}+{\beta }_2{d}_{i,t-1}+{\epsilon }_{it}.$$

hence, changes in $d_{t-1}$ can have both a direct and indirect effect on $Y_{it}$.

$$\left\{Y_{it}({\mathbf{d}}_{1:t}):t=1,\dots ,T\right\}\perp \perp {\mathbf{D}}_{i,1:t}\mid X_{i,0}.$$

This relates to linear panel models of the form

$$Y_{it}={\beta }_0+{\beta }_1{D}_{it}+{\beta }_2{D}_{i,t-1}+{\eta }_{it}.$$

where strict exogeneity $\mathbb{E}\left[{\eta }_{it}|{\mathbf{D}}_{i,1:T}\right]=\mathbb{E}\left[{\eta }_{it}\right]=0$ is assumed.

For every treatment history ${\mathbf{d}}_{1:T}$ and period $t$,

$$\left\{Y_{is}({\mathbf{d}}_{1:s}):s=1,\dots ,T\right\}\perp \perp {\mathbf{D}}_{i,1:t}\mid {\mathbf{V}}_{it}.$$

where ${\mathbf{V}}_{it}$ is a set of covariates such as $\left\{Y_{i,t-1},D_{i,t-1},{\mathbf{X}}_{it}\right\}$.

This relates to sequential exogeneity in panel models

$$\mathbb{E}\left[{\epsilon }_{it}\mid {\mathbf{D}}_{i,1:t},{\mathbf{X}}_{i,1:t},Y_{i,1:t-1}\right]=\mathbb{E}\left[{\epsilon }_{it}\mid D_{it},{\mathbf{V}}_{it}\right]=0.$$

Under sequential ignorability, an ADL approach would be to write the outcome regression with time-varying covariates

$$Y_{it}={\beta }_0+\alpha Y_{i,t-1}+{\beta }_1{D}_{it}+{\beta }_2{D}_{i,t-1}+{\mathbf{X}}_{it}^{\prime }\delta +{\epsilon }_{it}.$$

This generates post-treatment bias because ${\mathbf{X}}_{it}$ may be affected by ${\mathbf{D}}_{i,1:t-1}$.

Define the impulse response functions (‘blip-down’ functions) as

$$b_t({\mathbf{d}}_{1:t},j):=\mathbb{E}\left[Y_{it}({\mathbf{d}}_{1:t-j},0_j)-Y_{it}({\mathbf{d}}_{1:t-j-1},0_{j+1})\mid {\mathbf{D}}_{1:t-j}={\mathbf{d}}_{1:t-j}\right].$$

which is the effect of a change from $0$ to $d_{t-j}$ in terms of the treatment on the outcome at time $t$, conditional on treatment history up to time $t-j$.

These functions are parametrised as a function of lag length

$$b_t({\mathbf{d}}_{1:t},j;\mathbf{\gamma })={\gamma }_{1j}{d}_{t-j}+{\gamma }_{2j}{d}_{t-j}{d}_{t-j-1}+\dots$$

This then allows us to construct blipped-down / demediated outcomes

$${\tilde{Y}}_{it}^j=Y_{it}-\sum _{s=1}^{j-1}{\gamma }_s{D}_{i,t-s}.$$

Intuitively, this transformation subtracts off the effects of $j$ lags of treatment, creating an estimate of the counterfactual level of the outcome at time $t$ if the treatment had been set to $0$ for $j$ periods before $t$. Under sequential ignorability, the transformed outcome ${\tilde{Y}}_{it}^j$ has the same expectation as the counterfactual $Y_{it}({\mathbf{d}}_{1:t-j},0_j)$, and can be used to construct ${\tilde{Y}}_{it}^{j+1}$ by modelling the relationship between ${\tilde{Y}}_{it}^j$ and $D_{i,t-j}$ to estimate the lagged effect for $j+1$. This is recursive, hence the ‘nested’.

Sequential g-estimation can be used to estimate effects. Suppose we’re interested in the contemporaneous effect and the first-lagged effect and we adopt an impulse response function $b_t({\mathbf{d}}_{1:t,j;\gamma })={\gamma }_j{d}_{t-j}$ for both these effects. We assume sequential ignorability conditional on ${\mathbf{V}}_{it}:=\left\{D_{i,t-1},Y_{i,t-1},{\mathbf{X}}_{it}\right\}$. Sequential g-estimation proceeds as follows

For $j=0$ regress the un-transformed outcome on $\left\{D_{it},D_{i,t-1},Y_{i,t-1},{\mathbf{X}}_{it}\right\}$ as in an ADL model. If this is correctly specified, we estimate the blip-down parameter ${\gamma }_0$ (contemporaneous effect) correctly. We use ${\hat{\gamma }}_0$ to construct the one-lag blipped-down outcome ${\hat{Y}}_{it}^1=Y_{it}-{\hat{\gamma }}_0{D}_{it}$ This blipped-down outcome would be regressed on $\left\{D_{i,t-1},D_{i,t-2},Y_{i,t-2},{\mathbf{X}}_{i,t-1}\right\}$ to estimate the next blip-down parameter ${\gamma }_1$ (the first lagged effect) (repeat for further lags, standard error estimated via block-bootstrap)

To specify a marginal structural model, we choose a potential outcome lag length and write a model for the marginal model of those potential outcomes in terms of treatment history

$$\mathbb{E}\left[Y_{it}({\mathbf{d}}_{1:t})\right]=g({\mathbf{d}}_{1:t};\beta ).$$

for example, for a contemporaneous and two lagged effects, we write $\mathbb{E}\left[Y_{it}({\mathbf{d}}_{t-2:t})\right]=g({\mathbf{d}}_{t-2:t};\beta )$, marginalising over further lags and covariates.

The average causal effect is then

$${\tau }^{\text{msm}}:=g({\mathbf{d}}_{1:t};\beta )-g({\mathbf{d}}_{1:t}^{\prime };\beta ).$$

This motivates an IPW approach where weights are constructed as

$${\widehat{\text{SW}}}_{it}:=\prod _{s=1}^t\frac{\hat{\mathbb{P}}(D_{is}\mid D_{i,s-1},\hat{\gamma })}{\hat{\mathbb{P}}(D_{is}\mid {\mathbf{X}}_{is},Y_{i,s-1},D_{i,s-1},\hat{\alpha })}.$$

where the denominator of each term is the product of the predicted probability of observing unit $i$‘s observed treatment status conditional on covariates that satisfy conditional ignorability. Multiplying this over time produces the probability of seeing this unit’s treatment history conditional on the past.

These weights can be used in a regression of the form

$$g({\mathbf{d}}_{t-2:t};\beta )={\beta }_0+{\beta }_1{d}_t+{\beta }_2{d}_{t-1}+{\beta }_3{d}_{t-2}.$$

Decomposition Methods

Basic idea of decomposition

$${\mathbb{F}}_M(y)-{\mathbb{F}}_F(y)=\int {\mathbb{F}}_M(y\mid x)f_M(x)dx-\int {\mathbb{F}}_F(y\mid x)f_F(x)dx.$$ $${\mathbb{F}}_M(y)-{\mathbb{F}}_F(y)=\int [{\mathbb{F}}_M(y\mid x)-{\mathbb{F}}_F(y\mid x)]f_M(x)dx+\int {\mathbb{F}}_M(y\mid x)[f_M(x)-f_F(x)]dx.$$

Oaxaca-Blinder Decomposition

$${\overline{y}}_A-{\overline{y}}_B={\overline{x}}_A^{\prime }{\beta }_A-{\overline{x}}_B^{\prime }{\beta }_B={\overline{x}}_B^{\prime }({\beta }_A-{\beta }_B)+({\overline{x}}_A-{\overline{x}}_B{)}^{\prime }{\beta }_A.$$

We consider two groups, $A$ and $B$, and an outcome $Y$, and a vector of predictors $\mathbf{x}$. Main question for decomposition is how much of the mean outcome difference [or another summary statistic / quantile of CDF] is accounted for by group differences in the predictors $\mathbf{x}$. The Oaxaca-Blinder decomposition refers to the following decompositions:

$$R:=\mathbb{E}\left[Y_A\right]-\mathbb{E}\left[Y_B\right]=\mathbb{E}{\left[{\mathbf{x}}_A\right]}^{\prime }{\mathbf{\beta }}_A-\mathbb{E}{\left[{\mathbf{x}}_B\right]}^{\prime }{\mathbf{\beta }}_B.$$ $$R=\underset{\text{Explained}}{\underbrace{(\mathbb{E}\left[{\mathbf{x}}_A\right]-\mathbb{E}\left[{\mathbf{x}}_B\right]{)}^{\prime }{\mathbf{\beta }}_A}}+\underset{\text{Unexplained}}{\underbrace{\mathbb{E}{\left[{\mathbf{x}}_B\right]}^{\prime }({\mathbf{\beta }}_A-{\mathbf{\beta }}_B)}}.$$ $$R=\underset{\text{endowments}}{\underbrace{(\mathbb{E}\left[{\mathbf{x}}_A\right]-\mathbb{E}\left[{\mathbf{x}}_B\right]{)}^{\prime }{\mathbf{\beta }}_B}}+\underset{\text{discrimination}}{\underbrace{\mathbb{E}{\left[{\mathbf{x}}_B\right]}^{\prime }({\mathbf{\beta }}_A-{\mathbf{\beta }}_B)}}+\underset{\text{interaction}}{\underbrace{(\mathbb{E}\left[{\mathbf{x}}_A\right]-\mathbb{E}\left[{\mathbf{x}}_B\right]{)}^{\prime }({\mathbf{\beta }}_A-{\mathbf{\beta }}_B)}}.$$ $$R=\underset{\text{Explained}}{\underbrace{(\mathbb{E}\left[{\mathbf{x}}_A\right]-\mathbb{E}\left[{\mathbf{x}}_B\right]{)}^{\prime }{\mathbf{\beta }}^{\ast }}}+\underset{\text{Unexplained}}{\underbrace{\mathbb{E}{\left[{\mathbf{x}}_A\right]}^{\prime }({\mathbf{\beta }}_A-{\mathbf{\beta }}^{\ast })+\mathbb{E}{\left[{\mathbf{x}}_B\right]}^{\prime }({\mathbf{\beta }}^{\ast }-{\mathbf{\beta }}_B)}}.$$ $${\mathbf{\delta }}_A:={\mathbf{\beta }}_A-{\mathbf{\beta }}^{\ast },{\mathbf{\delta }}_B:={\mathbf{\beta }}_B-{\mathbf{\beta }}^{\ast }.$$

Detailed Decomposition

To examine the ‘contribution’ of each variable to the observed gap, estimate

$$y_i=\sum _{j=1}^k{x}_{ji}{\beta }_j+\sum _{j=1}^k{d}_i{x}_{ji}{\delta }_j+{\epsilon }_i,d_i:=\{\begin{array}{ll}1,& i\in B,\\ 0,& \text{otherwise}.\end{array}$$

so, ${\beta }_j$ is the coefficient for group $A$, and ${\beta }_j+{\delta }_j$ is the coefficient for group $B$. A t-test for ${\delta }_j$ is used to establish whether a variable is a source of the observed gap. The contribution of each variable to the explained part is

$$c_k^{\ast }=\frac{({\overline{x}}_k^A-{\overline{x}}_k^B){\hat{\beta }}_k^A}{({\overline{\mathbf{x}}}^A-{\overline{\mathbf{x}}}^B){\hat{\mathbf{\beta }}}^A}.$$

Let outcome models be linear $Y_i=X_i{\mathbf{\beta }}_1+{\nu }_{1i}\text{ if}{W}_i=1$ and $Y_i=X_i{\mathbf{\beta }}_0+{\nu }_{0i}\text{ if}{W}_i=0$ where $\mathbb{E}\left[{\nu }_{1i}\right]=\mathbb{E}{\left[\nu \right]}_{0i}=0$.

The difference in means decomposition is

$$\mathbb{E}\left[Y\mid W=1\right]-\mathbb{E}\left[Y\mid W=0\right]=\mathbb{E}\left[\mathbf{X}\mid W=1\right]({\mathbf{\beta }}_1-{\mathbf{\beta }}_0)+\left(\mathbb{E}\left[\mathbf{X}\mid W=1\right]-\mathbb{E}\left[\mathbf{X}\mid W=0\right]\right){\mathbf{\beta }}_0.$$ $$\mathbb{E}\left[Y\mid W=1\right]-\mathbb{E}\left[Y\mid W=0\right]=\mathbb{E}\left[Y^1-Y^0\mid W=1\right]+\left(\mathbb{E}\left[Y^0\mid W=1\right]-\mathbb{E}\left[Y^0\mid W=0\right]\right)={\tau }_{PATT}+\left(\mathbb{E}\left[Y^0\mid W=1\right]-\mathbb{E}\left[Y^0\mid W=0\right]\right).$$

Sloczynski: SATT can be estimated by running the following regression:

$Y_i=\alpha +\tau W_i+{\mathbf{X}}_i^{\prime }\mathbf{\beta }+\psi W_i({\mathbf{X}}_i-{\overline{\mathbf{X}}}_1)+{\epsilon }_i$

Kline (2011) shows that this is ‘doubly robust’ and equivalent to a reweighting estimator based on the weights

$$w(\mathbf{x}):=\frac{d{\mathbb{F}}_{\mathbf{X}\mid W=1}(\mathbf{x})}{d{\mathbb{F}}_{\mathbf{X}\mid W=0}(\mathbf{x})}=\frac{1-\rho }{\rho }\frac{e(\mathbf{x})}{1-e(\mathbf{x})}.$$

where $\rho :=\mathbf{Pr}\left(W_i=0\right)$ is the treated share.

Distributional Regression

Section based on counterfactual distribution decomposition methods.

Let $F_{X_k}$ denote the distribution of job-relevant characteristics (education, experience, etc.) for men when $k=m$ and for women when $k=w$. Let $F_{Y_j|X_j}$ denote the conditional distribution of wages given job-relevant characteristics for group $j\in \{w,m\}$, which describes the stochastic wage schedule that a given group faces. Using these distributions, we can construct $F_{⟨j\mid k⟩}$, the distribution of wages for group $k$ facing group $j$‘s wage schedule as

$F_{⟨j\mid k⟩}(y)=\int F_{Y_j\mid X_j}(y\mid x)d{F}_{X_k}(x),y\in {\mathcal{Y}}_j$

For example, $F_{⟨0\mid 0⟩}$ is the distribution of wages for men who face men’s wage schedule, and $F_{⟨1\mid 1⟩}$ is the distribution of wages for women who face women’s wage schedule, which are both observed distributions. We can also study $F_{⟨0\mid 1⟩}$, the counterfactual distribution of wages for women if they faced the men’s wage schedule $F_{Y_0\mid X_0}$.

$$F_{Y⟨0\mid 1⟩}(y)\equiv {\int }_{{\mathcal{X}}_1}{F}_{Y_0\mid X_0}(y\mid x)d{F}_{X_1}(x).$$

is the counterfactual distribution constructed by integrating the conditional distribution of wages for men with respect to the distribution of characteristics for women.

We can Interpret $F_{Y⟨0|1⟩}$ as the distribution of wages for women in the absence of gender discrimination, although it is predictive and cannot be interpreted as causal without further (strong) assumptions.

$$F_{Y⟨1\mid 1⟩}^{\leftarrow }-F_{Y⟨0\mid 0⟩}^{\leftarrow }=\underset{\text{structure}}{\underbrace{\left[F_{Y⟨1\mid 1⟩}^{\leftarrow }-F_{Y⟨0\mid 1⟩}^{\leftarrow }\right]}}+\underset{\text{composition}}{\underbrace{\left[F_{Y⟨0\mid 1⟩}^{\leftarrow }-F_{Y⟨0\mid 0⟩}^{\leftarrow }\right]}}.$$

Assumptions for Causal Interpretation

Under conditional exogeneity / selection on observables, CE can be interpreted as causal effects. Sec 2.3 in ECTA 2013 paper spells this out in detail. Let $(Y_j^{\ast }:j\in \mathcal{J})$ be the vector of potential outcomes for various values of a policy $j\in \mathcal{J}$, and $X$ be a vector of covariates. Let $J$ denote the random variable that describes the realised policy and let $Y:=Y_J^{\ast }$ denote the realised outcome variable. When $J$ is not randomly assigned, the distribution of $Y|J=j$ may differ from the distribution of $Y_j^{\ast }$. However, under conditional exogeneity, the distribution of $Y|X,J=j$ and $Y_j^{\ast }|X$ agree, and the observed conditional distributions have a causal interpretation, and so do counterfactual distributions generated from these conditionals by integrating out $X$.

Let $F_{Y_j^{\ast }|J}(y|k)$ denote the distribution of the potential outcome $Y_j^{\ast }$ in the population with $J=k\in \mathcal{J}$. The causal effect of exogenously changing the policy from $l$ to $j$ on the distribution of the potential outcome in the population with the realised policy $J=k$ is $F_{Y_j^{\ast }|J}(y|k)-F_{Y_l^{\ast }|J}(y|k)$. Under conditional exogeneity, for any $j,k\in \mathcal{J}$, the counterfactual distribution $F_{Y⟨j|k⟩}(y)$ exactly corresponds to $F_{Y_j^{\ast }|J}(y|k)$, and hence the causal effect of exogenously changing the policy from $l$ to $j$ in the population with $J=k$ corresponds to the CE of changing the conditional distribution from $l$ to $j$, that is

$F_{Y_j^{\ast }|J}(y|k)-F_{Y_l^{\ast }|J}(y|k)=F_{Y⟨j|k⟩}(y)-F_{Y⟨l|k⟩}(y)$

Conditional exogeneity assumption for this section:

$$(Y_j^{\ast }:j\in \mathcal{J})\perp \perp J\mid X$$

$\mathcal{K}$ groups that partition the sample. For each population $k,\exists X_k\in {\mathbb{R}}^d$ and outcome $Y_k$. Covariate vector is observable in all populations, but the outcome is only observable in populations $j\in \mathcal{J}\subset \mathcal{K}$. Let $F_{X_k}$ denote the covariate distribution in the population $k\in \mathcal{K}$, and $F_{Y_j|X_j}$ and $Q_{Y_j|X_j}$ denote the conditional distribution and quantile functions in population $j\in \mathcal{J}$. We denote the support of $X_k$ by ${\mathcal{X}}_k\subset {\mathbb{R}}^{d_x}$ and the region of interest $Y_j$ by ${\mathcal{Y}}_j\subseteq \mathbb{R}$. We refer to $j$ as the reference population and $k$ as the counterfactual population.

The reference and counterfactual populations in the wage example correspond to different groups. We can also generate counterfactual populations by artificially transforming a reference population. We can think of $X_k$ as being created through a known transformation of $X_j$: $X_k=g_k(X_j),\text{where}{g}_k:{\mathcal{X}}_j\to {\mathcal{X}}_k$

Counterfactual distribution and quantile functions are formed by combining the conditional distribution in population $j$ with the covariate distribution in population $k$, namely:

$$F_{Y⟨j\mid k⟩}(y)\equiv {\int }_{{\mathcal{X}}_k}{F}_{Y_j\mid X_j}(y\mid x)d{F}_{X_k}(x),y\in {\mathcal{Y}}_j.$$ $$Q_{Y⟨j\mid k⟩}(\tau )\equiv F_{Y⟨j\mid k⟩}^{\leftarrow }(\tau ),\tau \in (0,1).$$

where $(j,k)\in \mathcal{J}\mathcal{K}$ and $F_{Y⟨j|k⟩}^{\leftarrow }(\tau )=\inf \{y\in {\mathcal{Y}}_j:F_{Y⟨j|k⟩}(y)\ge \tau \}$ is the left-inverse function of $F_{Y⟨j|k⟩}$.

The main interest lies in the quantile effect (QE) function, defined as the difference of the two counterfactual quantile functions over a set of quantile indexes $\mathcal{T}\subset (0,1)$

$$\Delta (\tau )=Q_{Y⟨j\mid k⟩}(\tau )-Q_{Y⟨j\mid j⟩}(\tau ),\tau \in \mathcal{T}.$$

Estimation of Conditional distribution

$$F_{Y_j\mid X_j}(y\mid x)\equiv {\int }_{(0,1)}{1}_{Q_{Y_j\mid X_j}(u\mid x)\le y}du.$$

• method = "qr" default implements

$${\hat{F}}_{Y_j\mid X_j}(y\mid x)=\epsilon +{\int }_{(ϵ,1-ϵ)}{1}_{x^{\prime }{\hat{\beta }}_j(u)\le y}du.$$

where $\epsilon$ is a small constant that avoids estimation of tail quantiles, and $\hat{\beta }(u)$ is the quantile regression estimator

$${\hat{\beta }}_j(u)=\underset{b\in {\mathbb{R}}^{d_x}}{\mathrm{argmin}}\sum _{i=1}^{n_j}\left(u-1_{Y_{ji}\le X_{ji}^{\prime }b}\right)\left(Y_{ji}-X_{ji}^{\prime }b\right).$$

• method = "logit" implements the distribution regression estimator of the conditional distribution with the logistic link function

${\hat{F}}_{Y_j|X_j}(y|x)=\Lambda (x^{\prime }\hat{\beta }(y))$

where $\Lambda$ is the standard logistic CDF and $\hat{\beta }(y)$ is the distribution regression estimator

$$\hat{\beta }(y)\equiv \underset{b\in {\mathbb{R}}^{d_x}}{\mathrm{argmax}}\sum _{i=1}^{n_j}\left(1_{Y_{ji}\le y}\log \Lambda (X_{ij}^{\prime }b)+1_{Y_{ij}>y}\log \Lambda (-X_{ji}^{\prime }b)\right).$$

Causal Directed Acyclic Graphs

based on ,Pearl (2009), Morgan and Winship (2014), Cunningham (2020).

For an undirected graph between $X,Y,\text{ and }Z$, there are four possible directed graphs:

$X\to Y\to Z$ (a chain)

$X\leftarrow Y\leftarrow Z$ (another chain)

$X\leftarrow Y\to Z$ (a fork on Y)

$X\to Y\leftarrow Z$ (collision on Y)

With the fork or either chain, we have $X\perp \perp Z|Y$. However, With a collider, $X/\perp \perp Z|Y$.

Causal effect of $X$ on $Y$ is written $\mathbf{Pr}(Y|\mathsf{do}(X=x))$. Basic idea is condition on adequate controls (i.e. not every observed control). Here, controlling for $U$ is unnecessary and would bias the estimate of $\mathbf{Pr}(Y|\mathsf{do}(X=x))$.

Basics / Terminology

A backdoor path is a non-causal path from $A$ to $Y$. They are ‘backdoor’ because they flow backwards out of $A$: all of these paths point into $A$.

Here, $A\leftarrow U\to Y$, where $U$ is a common cause for treatment and the outcome. So, $U$ is a confounder.

A worse problem arises with the following DAG, where dotted lines indicate that $U$ is unobserved. Because $U$ is unobserved, this backdoor path is open.

Colliders, when left alone, always close a backdoor path. Conditioning on them, however, opens a backdoor path, and yields biased estimates of the causal effect of $A$ on $Y$.

Common colliders are post-treatment controls $A\to C\leftarrow Y$

Another insidious type of collider is of the form $A\leftarrow \cdots \to C\leftarrow \cdots \to Y$, where $C$ is typically a lagged outcome.

Vector of measured controls $\mathbf{S}$ satisfies the backdoor criterion if (i) $S$ blocks every path from $A$ to $Y$ that has an arrow into $A$ (i.e. blocks the back door) and (ii) no node in $S$ is a descendant of $A$. Then,

$\mathbf{Pr}\left(Y|do(A=a)\right)={\sum }_s\mathbf{Pr}\left(Y|A=a,S=s\right)\mathbf{Pr}\left(S=s\right)$

Which is the same as the subclassification estimator. The conditional Expectation $\mathbb{E}\left[Y|A=a,S=s\right]$ can be computed using a nonparametric regression / ML algorithm of choice.

$\mathbf{M}$ satisfies the frontdoor criterion if (i) $M$ blocks all directed paths from $A$ to $Y$, (ii) there are no unblocked back-door paths from $A$ to $M$, and (iii) $A$ blocks all backdoor paths from $M$ to $Y$.

Then,

$\mathbf{Pr}\left(Y|do(A)\right)=\underset{\text{Pr(M∣do(A))}}{\underbrace{\sum _M\mathbf{Pr}\left(M=m|A=a\right)}}\underset{\text{Pr(Y∣M,do(A))}}{\underbrace{\sum _{a^{\prime }}\mathbf{Pr}\left(Y|A=a^{\prime },M=M\right)Pr(A=a^{\prime })}}$

The above DAG in words

The only way $A$ influences $Y$ is through $M$, so there is no arrow bypassing $M$ between $X$ and $Y$. In other words, $M$ intercepts all directed paths from $A$ to $Y$.

Relationship between $A$ and $M$ is not confounded by unobservables - i.e. no back-door paths between A and M.

Conditional on $A$, the relationship between $M$ and $Y$ is not confounded, i.e. every backdoor path between $M$ and $Y$ has to be blocked by $A$.

With a single mediator $M$ that is not caused by $U$, the ATE can be estimated by multiplying estimates $\hat{\gamma }\times \hat{\delta }$ .

The FDC estimates the ATE because it decomposes a reduced-form relationship that is not causally identified into two causally identified relationships.

Implementation through linear regressions:

$$M_i=\kappa +\gamma A_i+{\omega }_i.$$ $$Y_i=\alpha +\delta M_i+\psi A_i+{\nu }_i.$$

Since $\mathbb{E}\left[M|A\right]=\gamma$ is identified, $\mathrm{Cov}\left({\omega }_i{A}_i\right)=0$ in the first-stage equation and $\mathrm{Cov}\left(M,\nu \right)=0$ in the second-stage equation. Assume $\psi =0$. Then, write

${\tau }_{\text{FDC}}=\mathbb{E}\left[Y|do(A)\right]=\hat{\delta }\times \hat{\gamma }$

Mediation Analysis

Pearl (2001), Robins(2003)

Consider SRS where we observe $(D_i,M_i,{\mathbf{X}}_i,Y_i)$, where $D_i$ is a treatment indicator, $M_i$ is a mediator, ${\mathbf{X}}_i$ is a vector of pre-treatment controls, and $Y_i$ is the outcome. The supports are $\mathcal{M},\mathcal{X},\mathcal{Y}$ respectively. $\mathbf{X}$s are partialled out.

Let $M_i(d)$ denote potential value for the mediator under treatment status $D_i=d$. The outcome $Y_i(d,m)$ is the potential outcome for unit $i$ when $D_i=d,M_i=m$. The observed variables can be written as $M_i=M_i(D_i),Y_i=Y_i(D_i,M_i(D_i))$.

d, a used interchangeably for treatment.

$$\left\{Y_i(d,m),M_i(d)\right\}\perp \perp D_i\mid X_i=x.$$ $$Y_i(d,m)\perp \perp M_i(d)\mid D_i,X_i=x.$$

$\forall d^{\prime },d\in \left\{0,1\right\}\text{ and }(m,x)\in \mathcal{M}\times \mathcal{X}$

This requires the treatment to be conditionally independent of the potential mediator states and outcomes given X, ruling out unobserved confounders jointly affecting the treatment on the one hand and the mediator and/or the outcome on the other hand conditional on the covariates. (5) postulates independence between the counterfactual outcome and mediator values ‘across-worlds’.

Effectively, Need $M$ to be randomly assigned (approx).

${\text{NIE}}_i(d)\equiv {\delta }_i(d):=Y_i(d,M_i(1))-Y_i(d,M_i(0))$

Difference in $Y$ holding treatment status constant, and varying the mediator. Sample Average: Average Causal Mediation Effect (ACME)

$$\overline{\delta }(d):=\mathbb{E}\left[{\delta }_i(d)\right]=\mathbb{E}\left[Y_i(d,M_i(1))-Y_i(d,M_i(0))\right].$$ $${\text{NDE}}_i(d)\equiv {\theta }_i(d):=Y_i(1,M_i(d))-Y_i(0,M_i(d)).$$

Difference in $Y$ holding mediator constant, and varying the treatment.

$${\tau }_i=Y_i(1,M_i(1))-Y_i(0,M_i(0)).$$ $${\tau }_i=\left(Y_i(1,M_i(1))-Y_i(0,M_i(1))\right)+\left(Y_i(0,M_i(1))-Y_i(0,M_i(0))\right)={\theta }_i(1)+{\delta }_i(0).$$ $${\tau }_i=\left(Y_i(1,M_i(0))-Y_i(0,M_i(0))\right)+\left(Y_i(1,M_i(1))-Y_i(1,M_i(0))\right)={\theta }_i(0)+{\delta }_i(1).$$ $${\tau }_i={\delta }_i(d)+{\theta }_i(1-d).$$

NDE conditions on potential mediator effects.For CDE, we set mediator at a prescribed value $m$.

${\text{CDE}}_i(d,d^{\prime },m)=Y_i(d,m)-Y_i(d^{\prime },m)m\in \mathcal{M}$

Difference between NDE and CDE is what value mediator is fixed at. Restated:

${\psi }_i(d,d^{\prime },m)=Y_i(d,m)-Y_i(d^{\prime },m)$

Effect of changing the treatment while fixing the value of the mediator at some level $m$.

$\overline{\psi }(d,d^{\prime },m)=\mathbb{E}\left[Y_i(d,m)-Y_i(d^{\prime },m)\right]$

Decomposing total effect with binary mediator

$$\tau (d,d^{\prime })=\text{ACDE}(d,d^{\prime },0)+\text{ANIE}(d,d^{\prime })+\mathbb{E}\left[M(a^{\prime })\left[\text{CDE}(d,d^{\prime },1)-\text{CDE}(d,d^{\prime },0)\right]\right].$$

Assume linear models for mediator $M=\psi T+U_m$ and $Y=\beta T+\gamma M+U_Y$.

Then fit the following regressions

$$Y_i={\alpha }_1+\tau D_i+{\epsilon }_{i1}.$$ $$M_i={\alpha }_2+\psi D_i+{\epsilon }_{i2}.$$ $$Y_i={\alpha }_3+\beta D_i+\gamma M_i+{\epsilon }_{i3}.$$

Baron and Kenny (1986) suggest testing $\tau =\psi =\beta =0$. If all nulls rejected, Mediation effect $\overline{\delta }=\psi \gamma$. Equivalently, mediation effect is $\tau -\beta =\psi \times \gamma$. Estimate variance using bootstrap / delta method.

Assume selection on observables w.r.t. D, M.

Huber(2014)

Average direct effect identified by

$$\theta (d)=\mathbb{E}\left[\left(\frac{Y\cdot D}{\mathbf{Pr}\left(D=1\mid M,X\right)}-\frac{Y\cdot (1-D)}{1-\mathbf{Pr}\left(D=1\mid M,X\right)}\right)\cdot \frac{\mathbf{Pr}\left(D=d\mid M,X\right)}{\mathbf{Pr}\left(D=d\mid X\right)}\right].$$

Average Indirect Effect identified by

$$\delta (d)=\mathbb{E}\left[\frac{Y\cdot 1_{D=d}}{\mathbf{Pr}\left(D=d\mid M,X\right)}\left(\frac{\mathbf{Pr}\left(D=1\mid M,X\right)}{\mathbf{Pr}\left(D=1\mid X\right)}-\frac{1-\mathbf{Pr}\left(D=1\mid M,X\right)}{1-\mathbf{Pr}\left(D=1\mid X\right)}\right)\right].$$

implemented in causalweight::medweight.