15. Selection on Observables and Matching

Non-parametric estimation aims to estimate an unknown quantity while making as few assumptions about the data generating process (DGP) as possible. We will think about the following two problems: - density estimation - regression: local constant, local linear

Both problems can be solved by non-parametric methods, which are basically local averaging methods. And the key question is: how to define local.

Sample analog of density. For a random variable \(X\), the density \(f(x)\) is defined as $$f(x)=\lim_{h\to0}\frac{F(x+h)-F(x-h)}{2h}=\lim_{h\to0}\frac{\mathbb P(x-h

Kernel. Generalize the histogram estimator to a family of estimators of density: $$\hat f(x)=\frac1{Nh}\sum_{i=1}^N K\Big(\frac{X_i-x}h\Big)$$ - the weighting function \(K(\cdot)\) is called a kernel function; - \(h\) is a smoothing parameter called the bandwidth.

For \(\tilde f\to f\) we require \(Nh\to\infty\) and \(h\to0\). We usually assume \(K(\cdot)\) satisfies: 1. Symmetry: \(K(z)=K(-z)\). 2. Integrates to 1: \(\int K(z)\,dz=1\). 3. Zero "mean" property: \(\int z\cdot K(z)\,dz=0\) (15.2). 4. Finite second moment property: \(\int z^2\cdot K(z)\,dz=\kappa_2<\infty\) (15.3).

Desired properties: \(\int\hat f(x)\,dx=1\); \(\int x\hat f(x)\,dx=\frac1N\sum_{i=1}^N X_i\).

Estimation bias. Consider the bias of the estimator \(\hat f(x)\): $$\mathbb E[\hat f(x)]=\frac1h\mathbb E[K(z_i)]=\int K(z_i)f(u_i)\,d\tfrac1h u_i=\int K(z_i)f(z_ih+x)\,dz_i\tag{15.4}$$ which typically cannot be solved analytically, so we do the 2nd-order Taylor approximation of the last line of (15.4): $$\begin{aligned}\int K(z_i)f(z_ih+x)\,dz_i&\approx\int K(z_i)\big[f(x)+f'(x)z_ih+f''(x)z_i^2h^2\big]\,dz_i\\&\approx f(x)\underbrace{\int K(z_i)\,dz_i}_{=1}+f'(x)h\underbrace{\int K(z_i)z_i\,dz_i}_{=0}+f''(x)h^2\underbrace{\int K(z_i)z_i^2\,dz_i}_{=\kappa_2}\\&\approx f(x)+f''(x)h^2\kappa_2\end{aligned}$$ i.e. $$\mathbb E[\hat f(x)]-f(x)\approx\underbrace{f''(x)h^2\kappa_2}_{\text{bias}}\tag{15.5}$$ So \(\hat f(x)\) is a biased estimator, with bias approximately \(f''(x)h^2\kappa_2\) when \(h\) is small.

Estimation variance. Consider the variance of the estimator \(\hat f(x)\) (assume \(X_i\) i.i.d.): $$\operatorname{Var}(\hat f(x))=\frac1{Nh^2}\mathbb E\Big[K\Big(\tfrac{X_i-x}h\Big)^2\Big]-\frac1N\Big(\frac1h\mathbb E\Big[K\Big(\tfrac{X_i-x}h\Big)\Big]\Big)^2\tag{15.6}$$ Part 2 goes to 0 as \(N\to\infty\) (since \(\frac1h\mathbb E[K(\frac{X_i-x}h)]=\mathbb E[\hat f(x)]\approx f(x)\)); Part 1 does not necessarily go to 0 (since we only assumed \(Nh\to\infty\), not \(\frac1{Nh^2}\to0\)). Rearranging, $$\operatorname{Var}(\hat f(x))\approx\frac1{Nh}f(x)R(K)\tag{15.7}$$ where \(R(K)\equiv\int K(z_i)^2\,dz_i\) is called the roughness of the kernel.

Bias variance trade-off. From (15.5) and (15.7), at each amount \(N\) of data, lower \(h\) leads to smaller bias but higher variance. However, since bias does not depend on sample size but variance decreases in sample size, we can decrease the bandwidth \(h\) as the sample grows larger but at a slower rate than the growth of sample size. In this way, we can have both smaller bias and smaller variance.

We would always like to estimate the relationship between an outcome \(Y\) and an explanatory variable \(X\), generally characterized by $$g(x)\equiv\mathbb E[Y\mid X=x]$$ using the same technique as in non-parametric density estimation. From $$\mathbb E[Y\mid X=x]=\int y f(y\mid X=x)\,dy=\int y\frac{f(y,x)}{f(x)}\,dy\tag{15.8}$$ the sample analog is \(\hat{\mathbb E}[Y\mid X=x]=\int y\frac{\hat f(y,x)}{\hat f(x)}\,dy\). We already know how to estimate \(\hat f(x)\), so it only remains to estimate \(\hat f(y,x)\) — use the bivariate kernel \(K(u,v)=K_1(u)K_2(v)\): $$\hat f(y,x)=\frac1{Nh^2}\sum_{i=1}^N K\Big(\frac{X_i-x}h,\frac{Y_i-y}h\Big)=\frac1{Nh^2}\sum_{i=1}^N K_1\Big(\frac{X_i-x}h\Big)K_2\Big(\frac{Y_i-y}h\Big)$$ Finally we obtain the non-parametric estimator of \(g(x)\): $$\begin{aligned}\hat g(x)\equiv\hat{\mathbb E}[Y\mid X=x]&=\frac{\int y\hat f(y,x)\,dy}{\hat f(x)}=\frac{\frac1{Nh}\sum_{i=1}^N K_1(\frac{X_i-x}h)Y_i}{\frac1{Nh}\sum_{i=1}^N K_1(\frac{X_i-x}h)}\\&=\frac{\sum_{i=1}^N K_1(\frac{X_i-x}h)Y_i}{\sum_{i=1}^N K_1(\frac{X_i-x}h)}=\sum_{i=1}^N\omega_h(X_i,x)Y_i\end{aligned}$$ with weight \(\omega_h(X_i,x)=\frac{K_1(\frac{X_i-x}h)}{\sum_{i=1}^N K_1(\frac{X_i-x}h)}\). When \(K_1(z)=\frac12\mathbf 1\{x-h

For a specific kernel, the bandwidth \(h\) could be determined in a sense of optimality, e.g. leave-one-out cross-validation $$h^\star=\arg\max_h\sum_{i=1}^N\big(\hat g_{h,(-i)}(X_i)-Y_i\big)^2$$ where \(\hat g_{h,(-i)}(x)\) is estimated without using \((X_i,Y_i)\).

Theoretically \(g(x)=\mathbb E[Y\mid X=x]\) is point identified. However, when doing estimation, we have to allow for a small band (local neighborhood) around \(X=x\) and estimate the estimand (e.g. \(\text{ATE}(x)\)) using the data in that local neighborhood of \(X=x\) instead of only at \(X=x\). Then it becomes important how we believe the \(g(\tilde x)\) changes for \(\tilde x\in(x-h,x+h)\) slightly off \(X=x\). There come two kinds of beliefs: - Local constant regression: assume the estimand within each small enough local neighborhood is constant, which is only an approximation valid for small enough neighborhood. - Local linear regression: assume the estimand within each small enough local neighborhood is varying linearly with \(X\), which is still just an approximation valid for small enough neighborhood, but relaxes the constant assumption a little bit so should perform better in approximation than local constant regression.

Local constant regression. Kernel regression is equivalent to the local constant regression below: $$g(x)=\alpha,\qquad\hat\alpha=\arg\min_a\sum_{i=1}^N K_1\Big(\frac{X_i-x}h\Big)(Y_i-a)^2$$

Local linear regression. In each small band around \(X=x\), we relax \(g(x)=\alpha\) for local constant regression by allowing the estimand \(g(x)\) to vary linearly with \(X\) in each small neighborhood: $$g(x)\equiv\mathbb E[Y\mid X=x]=\alpha+\beta(X-x)\quad\text{for }X\in(x-h,x+h)$$ The solution \((\hat\alpha,\hat\beta)\): $$(\hat\alpha,\hat\beta)=\arg\min_{a,b}\sum_{i=1}^N K_1\Big(\frac{X_i-x}h\Big)(Y_i-a-b(X_i-x))^2$$ When the random variable \(X\) takes the exact value \(x\), the local linear regression gives \(\hat g(x)=\hat\alpha+\hat\beta(x-x)=\hat\alpha\).

Selection into treatment actually assumes that agents have knowledge about their potential outcomes \((Y_0,Y_1)\), and make a decision on whether to take treatment (\(D=0\) or \(1\)) depending on their knowledge of \((Y_0,Y_1)\), i.e. \((Y_0,Y_1)\not\perp D\) (dependence). In contrast, random assignment of treatment state (RCT) assumes that agents choose \(D\) without taking into consideration the knowledge of \(\{Y_d\}_{d\in\mathcal D}\), i.e. \((Y_0,Y_1)\perp D\) (independence).

We always use OLS regression to estimate \(\beta^{\text{OLS}}\). As discussed before, when \((Y_0,Y_1)\perp D\) (independence) is true, \(\beta^{\text{OLS}}\) is the ATE of the whole group: $$\beta^{\text{OLS}}=\mathbb E[Y\mid D=1]-\mathbb E[Y\mid D=0]=\mathbb E[Y_1\mid D=1]-\mathbb E[Y_0\mid D=0]\overset{\perp}{=}\mathbb E[Y_1]-\mathbb E[Y_0]=\mathbb E[Y_1-Y_0]$$ However, when \((Y_0,Y_1)\not\perp D\) (dependence) is true, we can only do the following decomposition: $$\begin{aligned}\mathbb E[Y\mid D=1]-\mathbb E[Y\mid D=0]&=\mathbb E[Y_1\mid D=1]-\mathbb E[Y_0\mid D=0]\\&=\underbrace{\mathbb E[Y_1\mid D=1]-\mathbb E[Y_0\mid D=1]}_{\text{ATT}}+\underbrace{\mathbb E[Y_0\mid D=1]-\mathbb E[Y_0\mid D=0]}_{\text{selection bias}}\\&=\underbrace{\mathbb E[Y_1\mid D=0]-\mathbb E[Y_0\mid D=0]}_{\text{ATU}}+\underbrace{\mathbb E[Y_1\mid D=1]-\mathbb E[Y_1\mid D=0]}_{\text{selection bias}}\end{aligned}\tag{15.10}$$ Selection bias is the part that comes from selection into treatment. For example in (15.10), it is possible that the treated select to be treated because they know they have higher potential outcome \(Y_1\) than the potential \(Y_1\) of those untreated, and then $$\underbrace{\mathbb E[Y_1\mid D=1]-\mathbb E[Y_1\mid D=0]}_{\text{selection bias}}>0\Rightarrow\beta^{\text{OLS}}>\text{ATU}$$ Selection bias is zero under random assignment (RCT), which implies \(\text{ATT}=\text{ATU}=\text{ATE}\) under RCT.

We could solve the problem of selection into treatment by conducting random assignment (RCT) in experiments, but, as we have discussed, experiment could have low external validity. So we can instead choose to assume selection on observables, which is a relaxation of random assignment (RCT).

Selection on observables means that agents select into treatment completely based on variables in the vector \(X\) that is observable to the researcher. So, conditional on observables, treatment is as good as randomly assigned.

Common support (overlapping condition) and balancing. We will also need the overlapping assumption \(\mathbb P(D=d\mid X)<1\) for all \(d\in\mathcal D\), which means that treated and untreated have common support in each matched cell \(X=x\) (both treated and untreated agents exist in each matched cell). Common support is crucial since we are first estimating the average treatment effect \(\text{ATE}(x)\) in each matched cell, and then averaging them up. If common support is broken, then the \(\text{ATE}(x)\) in that cell is not well-defined, and it ruins our overall effect estimation. We also need the balancing condition, which means that the treated and untreated should be comparable in size within each cell to reduce the variance of the estimation. Otherwise, even though the common support is satisfied, it could still be problematic if we have thousands of treated and only one or two untreated in some cells.

Identify counterfactual potential outcomes in expectation with conditional independence. Now we have the conditional version of random assignment: $$F_d(y\mid x)\equiv\mathbb P(Y_d\le y\mid X=x)=\mathbb P(Y_d\le y\mid D=d,X=x)=\mathbb P(Y\le y\mid D=d,X=x)$$ where the second line requires the overlapping condition. Then we can average over \(X\) to point identify the unconditional CDF of \(Y_d\): $$\begin{aligned}F_d(y)\equiv\mathbb P(Y_d\le y)&=\mathbb E_X[\mathbb P(Y_d\le y\mid D=d,X)]=\mathbb E_X[\mathbb P(Y\le y\mid D=d,X)]\\&\overset{\text{discrete}}{=}\sum_x\mathbb P(Y\le y\mid D=d,X=x)\mathbb P(X=x)\end{aligned}$$ which could thus identify \(\mathbb E[Y_d]=\int y\,dF_d(y)\).

The next question is when might conditional independence hold? There are two typical cases: 1. We have detailed information about the treatment assignment mechanism, which informs us that only a certain set of observable variables influence both potential outcomes and treatment selection. After taking out those variables (within groups matched based on those variables), the independence holds between potential outcomes and treatment selection. 2. We have extremely rich data. So, we can control for whatever we think is relevant in selection decision making.

Possibility of selection on unobservables: example from Neale and Johnson (1996, JPE). The question of interest is what is the effect of pure racial factor in the Black-White wage gap. To look at this, they regress earnings on race (\(D=1\): white; \(D=0\): black) and some covariates, and without test scores \(T\) (\(T=0\): high score; \(T=0\): low score). They find that the test scores play a large component in earnings. Conditional on test scores there doesn't seem to be a large difference in black and white wages, i.e. $$\mathbb E[Y\mid D=1,X\backslash T]-\mathbb E[Y\mid D=0,X\backslash T]\gg\mathbb E[Y\mid D=1,X\backslash T,T]-\mathbb E[Y\mid D=0,X\backslash T,T]$$ Does it mean that test score is the factor differentiating black from white in wages? Many found that this explanation was implausible. Suppose we regress \(Y\) on \(D\), \(X\backslash T\) and \(T\), in the regression $$Y=\alpha+\beta D+\gamma T+\eta X\backslash T+\varepsilon$$ we control for test score. However, there might exist some unobservables that also determine wage, such as effort \(E\). Maybe the black, as a disadvantaged group, need to put much more effort to achieve the same test score, and thus \(E\) and \(D\) are correlated. Also, clearly \(E\) and \((Y_0,Y_1)\perp D\mid X\) is ruined by the existence of unobservable \(E\), which is triggered by including \(T\) in \(X\). The lesson is that we should always be careful about the choices of control variables, and think twice about the story behind the conditional independence and possible unobservables in the story.

Bias and the inclusion of additional controls. Is it always true that the bias becomes smaller when you include more variables?

Suppose that \((Y_0,Y_1)\perp D\mid X\). Let \(X_1\subset X_2\subset X\). When it is possible to use \(X\), we can always reduce bias by adding more controls to \(X_1\) or \(X_2\) and make it become \(X\). However, when we cannot observe \(X\), it is not sure which one of \(X_1\) and \(X_2\) will give smaller bias. Including more is not always better: as we discussed before, if we throw a lot of variables in \(X\), it is very likely that we may be including outcome variables, which for sure increases the bias. Eliminating bias is just a limiting case.

Suppose that \(D\in\{0,1\}\) is binary. We then have that $$\begin{aligned}\text{ATE}&\equiv\mathbb E[Y_1-Y_0]=\mathbb E[\mathbb E[Y_1\mid X]-\mathbb E[Y_0\mid X]]=\mathbb E\big[\underbrace{\mathbb E[Y\mid D=1,X]}_{\text{data}}-\underbrace{\mathbb E[Y\mid D=0,X]}_{\text{data}}\big]\\&\overset{\text{discrete}}{=}\sum_x\big(\mathbb E[Y\mid D=1,X=x]-\mathbb E[Y\mid D=0,X=x]\big)\underbrace{\mathbb P(X=x)}_{\text{weights}}\\&\overset{\text{continuous}}{=}\int\big(\mathbb E[Y\mid D=1,X=x]-\mathbb E[Y\mid D=0,X=x]\big)\underbrace{f_X(x)}_{\text{weights}}\,dx\end{aligned}$$ (which follows by conditional independence and the overlapping support condition.) Likewise $$\text{ATT}\equiv\mathbb E[Y_1-Y_0\mid D=1]=\int\big(\mathbb E[Y\mid D=1,X=x]-\mathbb E[Y\mid D=0,X=x]\big)\underbrace{f_{X\mid D=1}(x)}_{\text{weights}}\,dx$$ $$\text{ATU}\equiv\mathbb E[Y_1-Y_0\mid D=0]=\int\big(\mathbb E[Y\mid D=1,X=x]-\mathbb E[Y\mid D=0,X=x]\big)\underbrace{f_{X\mid D=0}(x)}_{\text{weights}}\,dx$$ where $$f_{X\mid D=1}(x)=\frac{\mathbb P(D=1,X=x)}{\mathbb P(D=1)}=f(x)\frac{\mathbb P(D=1\mid X=x)}{\mathbb P(D=1)},\quad f_{X\mid D=0}(x)=f(x)\frac{\mathbb P(D=0\mid X=x)}{\mathbb P(D=0)}$$ Every item inside these expressions can be estimated by non-parametric methods (see §15.1).

Using \(p(X)\), we can rewrite the ATE in each matched cell \(X=x\), i.e. \(\text{ATE}(x)\), and then rewrite ATE, ATT and ATU as a weighted average of \(\text{ATE}(x)\). First look at \(\text{ATE}(x)\): since $$\mathbb E[YD\mid X=x]=\mathbb E[Y_1\mid D=1,X=x]\mathbb P(D=1\mid X=x)\overset{\text{CIA}}{=}\mathbb E[Y_1\mid X=x]p(x)$$ $$\Rightarrow\mathbb E[Y_1\mid X=x]=\frac{\mathbb E[YD\mid X=x]}{p(x)}=\mathbb E\Big[\frac{YD}{p(x)}\mid X=x\Big]$$ similarly \(\mathbb E[Y_0\mid X=x]=\mathbb E\big[\frac{Y(1-D)}{1-p(x)}\mid X=x\big]\). Then $$\begin{aligned}\text{ATE}(x)&=\mathbb E[Y_1-Y_0\mid X=x]=\mathbb E\Big[\frac{YD}{p(x)}\mid X=x\Big]-\mathbb E\Big[\frac{Y(1-D)}{1-p(x)}\mid X=x\Big]\\&=\mathbb E\Big[\frac{YD(1-p(x))-Y(1-D)p(x)}{p(x)(1-p(x))}\mid X=x\Big]=\mathbb E\Big[\frac{(D-p(x))Y}{p(x)(1-p(x))}\mid X=x\Big]\end{aligned}$$ Then in the continuous case we can rewrite $$\text{ATE}\equiv\mathbb E[\text{ATE}(x)]=\int\text{ATE}(x)\underbrace{f_X(x)}_{\text{weights}}\,dx=\int\frac{(D-p(x))Y}{p(x)(1-p(x))}\underbrace{f_X(x)}_{\text{weights}}\,dx$$ $$\text{ATT}\equiv\mathbb E[\text{ATE}(x)\mid D=1]=\int\frac{(D-p(x))Y}{p(x)(1-p(x))}\underbrace{\frac{\mathbb P(D=1\mid X=x)}{\mathbb P(D=1)}f(x)}_{\text{weights}}\,dx$$ $$\text{ATU}\equiv\mathbb E[\text{ATE}(x)\mid D=0]=\int\frac{(D-p(x))Y}{p(x)(1-p(x))}\underbrace{\frac{\mathbb P(D=0\mid X=x)}{\mathbb P(D=0)}f(x)}_{\text{weights}}\,dx$$ Every item inside these expressions can be estimated by non-parametric methods (§15.1).

Up to now, we are matching on \(X\), which is generally a vector with high dimensions. Theorem 15.1 implies that we can also match on the propensity score \(p(X)\) (a scalar). Different methods imply different identification and thus lead to different construction of estimators. Then, those estimators may have different efficiency, rates of convergence, and finite sample performance etc.

Since we will be using non-parametric methods in the estimation, it is very important that we put as many data points as possible within each cell of match to reduce both bias and variance. However, given the limited size of data, high dimensional control variables \(X\) need exponentially increasing amount of data points. For example, suppose that you have \(K=30\) binary covariates as control variables \(X\); then this would imply that \(\min N=2^{30+1}=2147483648\). This is the curse of dimensionality.

To overcome the curse we can do the following: 1. Regression analysis. (a) OLS regression cannot give us ATE, ATT, or ATU. Instead, it gives us a differently weighted average of the ATE in each matched cell; (b) see below for details. 2. Match on propensity score: solve the dimensionality problem by partly using parametric models. (a) This is similar to what we have done for matching on \(X\); (b) however, although we can avoid the curse of dimensionality for matching, we are still faced with this problem when estimating propensity score; (c) in practice, we could use probit or logit regression model instead of non-parametric methods to estimate propensity score and avoid the curse of dimensionality; (d) so, we could end up using non-parametric methods for averages for matched cells but parametric models for estimation of propensity score. 3. Inexact matching when doing non-parametric estimates. (a) First, define a metric on \(X\); (b) then, allow for a positive bandwidth of matching; (c) but we need to think about the optimal choice of bandwidth, taking into consideration of the trade-off between bias and variance. See more discussion on this in §15.1.

In conclusion, if we insist on using completely non-parametric method in estimation to make as less assumptions as possible, then the cost to pay is the dimensionality problem. But if we combine parametric and non-parametric methods (e.g. matching on propensity score), then the dimensionality problem can be solved and it is again a trade-off and balance between different objectives.

OLS Regression. Suppose we have found correct control variables \(X\) and add it into the OLS regression; then the coefficient estimand \(\beta^{\text{OLS}}\) for treatment \(D\) is a weighted average of \(\text{ATE}(x_i)\), but the weights are different from what we used for ATE, ATT and ATU, so \(\beta^{\text{OLS}}\) is an object different from ATE, ATT and ATU with special weights. To be more specific, consider the following saturated-in-\(X\) regression model: there are \(K\) possible discrete states for \(X_i\), \(\forall i\). By imposing the overlapping condition, we have common support for treated and untreated, i.e. \(X_i\) partitioned the agents into \(K\) subgroups, within which there are both treated and untreated agents. The model is: $$Y_i=\sum_{k=1}^K d_{ik}\alpha_k+\beta^{\text{OLS}}D_i+U_i,\qquad d_{ik}=\mathbf 1\{X_i=x_k\}\tag{15.14}$$ where: \(\alpha_k\) is the coefficient for \(X_k=x_k\); \(\beta^{\text{OLS}}\) is the regression estimand for treatment; CIA \((Y_{0i},Y_{1i})\perp D_i\mid X_i\) holds.

Clearly, in the OLS estimator, weight \(\omega_k\) is higher for the cell \(x_k\) in which variance of treatment (\(D_i\)) is higher. This makes sense because the cell with higher variance has more variations in \(D_i\) thus more accurate estimator and less residuals. Since OLS minimizes the sum of squared distance, this result is natural. Also note that \(\delta_{x_k}=\text{ATE}(x_k)\) can also be estimated by non-parametric methods.

There are also some other alternative methods including matching on different objects. But every method faces the same trade-off: more assumptions on the data generating process (parametric methods) or more severe dimensionality problem (non-parametric methods). Inside the dimensionality problem, we have the trade-off between bias and variance.