Distribution Regression and Counterfactual Analysis

14.382 L7. DISTRIBUTION REGRESSION AND COUNTERFACTUAL ANALYSIS Abstract. We use binary regressions to model conditional distributions of real outcomes given covariates and develop models for counterfactual analysis. An application to (predictive) gender discrimination effects in the labor market finds that these effects explain most of the observed difference between the distribution of wages for men and women. An application to racial differences in mental ability of young children finds that two thirds of the gap between black and white children is explained by differences in background characteristics, but the remaining unexplained gap is still significant in both economic and statistical terms. 1. Distribution Regression We use the binary regression framework to model and estimate the conditional distribution of an outcome given covariates. The outcome is a real-valued variable and can be either of the following: • • • • continuous (log wages), count (number of patents), nonnegative (durations, capital levels), discrete (rounded wages) or binary as in the previous lecture. The key, simple observation is that the conditional distribution of the outcome Y given the set of covariates X can be expressed as FY |X (y | x) = E[1{Y ≤ y} | X = x]. Accordingly, even when the outcome is not binary, we can always construct a collection of binary response variables, which record the events that the outcome falls bellow a set of thresholds: 1{Y ≤ y}, y ∈ T ⊂ R, where T is countable subset of R. For estimation and other practical purposes, we take T to be a finite collection of grid points. We can then use binary regressions for a collection of binary response variables to model the conditional distribution of Y given X, FY |X (y | x) = P(Y ≤ y | X = x) = Fy (B(x)0 β(y)), 1 2 where Fy is a link function, which may vary with the threshold y, B(X) is a dictionary of transformations of X including a 1 as the first entry, and β(y) is the parameter vector that may vary with the threshold y. This is the distribution regression model. The distribution regression model is quite flexible and nests a variety of classical models for conditional distribution functions. Example 1 (Classical Normal Regression Model). In the classical normal regression model, Y | X ∼ N (B(X)0 γ, σ 2 ), the conditional distribution of Y given X is FY |X (y | x) = Φ((y − B(x)0 γ)/σ), where Φ is the standard normal distribution. This conditional distribution is a special case of the distribution regression model with Fy equal to the probit link Φ, and β(y) = 0 )0 , for γ = (γ , γ 0 )0 . Note that the slopes here don’t vary with y. (y − γ1 , −γ−  1 −1 1 Example 2 (Cox Proportional Hazard Model). The Cox duration regression model is: FY |X (y | x) = 1 − exp(− exp(t(y) − B(x)0 γ)), where t(·) is an unknown monotonic transformation, is a common approach to model conditional distributions in duration and survival analysis. It has also been used to model non-negative outcomes, such as capital in (S, s) models and wages. It corresponds to the following location-shift representation: t(Y ) = B(X)0 γ + V, where V has an extreme value distribution and is independent of X. The model is called proportional hazard model, because the hazard rate, ∂t(y) ∂ d ln(1 − FY |X (y | x)) = − exp(t(y)) exp(B(x)0 γ), ∂y ∂y depends proportionally on exp(B(x)0 γ). The conditional distribution is a special case of the distribution regression model with Fy equal to the complementary log-log link, 0 )0 , for γ = (γ , γ 0 )0 . Note that Fy (u) = 1 − exp(− exp(u)), and β(y) = (t(y) − γ1 , −γ− 1 −1 1 the slopes here don’t vary with y, while distributional regression allows for slopes to be varying with y.  h(y | x) = Example 3 (Poissson Regression Model). The Poisson distribution is frequently used to model count variables taking values in 0, 1, 2, .... In the conditional (regression) version of the Poisson model, the conditional distribution of the count variable Y given X takes the form: y X

exp (B(x)0 γ)k exp (− exp (B(x)0 β)) = Q y, exp B(x)0 γ FY |X (y | x) = k! k=0 where Q called the incomplete Gamma function. Thus, the distribution regression model with link function Fy (u) = Q (y, exp(u)) and β(y) = γ nests the Poisson regression model. The Poisson regression model assumes that the same index governs the whole distribution. This assumption has been often challenged in applications. The zero-inflated Poisson 3 regression model makes one step in the direction of a more flexible model by allowing the coefficients to be different at 0. The distribution regression model does not restrict the heterogeneity of the coefficients at any level.  Given the conditional distribution we can also look at the conditional quantiles: FY←|X (u | x), u ∈ [0, 1], where ← denotes the left-inverse of the map y 7→ FY |X (y | x) on T . Here we define the left-inverse of a function G : T → [0, 1] on T as G← (u) := inf{t ∈ T : G(t) ≥ u} ∧ sup{t ∈ T }. (1.1) Given a graph of a distribution function t 7→ G(t) we can obtain the graph of the quantile function u 7→ G← (u) by simply flipping the axes and mirroring the resulting image. There are many ways to estimate the distribution regression models. Here we focus on the following method. We can estimate the conditional distribution by: ˆ FˆY |X (y|x) = Fy (B(x)0 β(y)), y ∈ T, ˆ where for each y ∈ T the estimator β(y) is the maximum likelihood estimator, β̂(y) ∈ arg max En [1(Yi ≤ y) ln Fy (B(Xi )0 b(y)) + 1(Yi > y) ln(1 − Fy (B(Xi )0 b(y)))], b(y)∈B ˆ where B is the parameter space for β(y). For example, β(y) is the probit estimator if we use the normal link Fy = Φ or the logit estimator if we use the logistic link Fy = Λ. We can ˆ handle inference on y 7→ Fy (B(x)0 β(y)) for y ∈ T by using the delta method in conjunction with the GMM formulation of the problem. Indeed, we can view estimation of θ0 = vec(β(y) : y ∈ T ) as a GMM problem with the score: g(Zi , θ) = vec (gy (Zi , b(y)) : y ∈ T ) , ∂  1(Yi ≤ y) ln Fy (B(Xi )0 b(y)) + 1(Yi > y) ln(1 − Fy (B(Xi )0 b(y))) , ∂b(y) which simply stacks the scores of many binary regressions. The joint parameter vector is θ = vec(b(y) : y ∈ T ), which simply stacks the parameters of many binary regressions. The ˆ ˆ map y 7→ Fy (B(x)0 β(y)), y ∈ T , is a smooth transformation of the estimators β(y), y ∈ T, so the delta method delivers the large sample properties of the estimators Fy (B(x)0 β̂(y)), y ∈ T . This also means that we can use the bootstrap for inference. gy (Zi , b(y)) = ˆ In practice, the map y 7→ Fy (B(x)0 β(y)) may be non-monotone, in which case we can rearrange it into monotone function by simply sorting the values of function in a nondecreasing order. This typically improves the finite-sample properties of the estimator. We discuss the rearrangement procedure in Section 3.2. 4 2. Unconditional and Counterfactual Distributions and Quantiles Once we get the conditional distribution, we can construct an unconditional distribution of the form: Z F (y) = FY |X (y|x)dM (x), y ∈ T, which is simply the conditional distribution integrated against a marginal distribution M for the covariate values. When M is the distribution of X, namely M = FX , we obtain the marginal distribution of Y : F (y) = FY (y), y ∈ T. We can also look at the marginal quantile functions Q = F ← which is the left-inverse of the distribution function y 7→ F (y) on T , as defined in (1.1). When M is not the distribution of X, then F is a counterfactual distribution which corresponds to a sampling experiment where covariates X 0 are sampled from M , but outcomes are sampled from the conditional distribution FY |X (·|X 0 ). This construction is very useful for counterfactual analysis. For example, let FXk denote the distribution of job-relevant characteristics (education, experience, etc.) for men when k = m, and women when k = w. Let FYj |Xj denote the conditional distributions of wages given job-relevant characteristics for group j ∈ {w, m}. This conditional distribution describes the stochastic wage schedule that a given group faces. Using these distributions we can construct Fhj |ki , the distribution of wages for group k facing group j wage schedule as Z Fhj |ki (y) = FYj |Xj (y |x)dFXk (x), y ∈ T. For example, Fhm|mi is the distribution of wages for men who face men’s wage schedule, and Fhw|wi is the distribution of wages for women who face women’s wage schedule. These are observed distributions. We can also look at Fhm|wi , the counterfactual distribution of wage for women if they would face the men’s wage schedule. We can interpret Fhm|wi as the distribution of wages for women in the absence of “gender discrimination”, although we do have to be careful and say that this is just a predictive distribution, which does not have a causal interpretation without further (strong) assumptions. We can use the counterfactual distributions to decompose the differences in the observed wage distributions: Fhm|mi − Fhw|wi = (Fhm|mi − Fhm|wi ) + (Fhm|wi − Fhw|wi ), composition effect discrimination effect where the first term on the right is the composition effect, which results from the populations of men and women having different distributions of job market characteristics, and 5 the second is the discrimination effect or price effect, which results from women facing different wage schedules than men. Analogous decomposition could be made for the observed quantile functions: Qhm|mi − Qhw|wi = (Qhm|mi − Qhm|wi ) + (Qhm|wi − Qhw|wi ) . composition effect (2.1) discrimination effect What we see above are the distribution and quantile versions of the Oaxaca-Blinder decomposition. In the empirical example based on U.S. data, which we describe below, we find that the distributions of wages for men and women are different, with the distribution of wages for men being shifted to the right by 20% or 30% along the horizontal axis, or the quantile function for men being shifted up by 20% or 30%. Almost all of the difference can be attributed to the discrimination/price effect. The composition effect is close to zero and is slightly negative due to the fact that men’s characteristics distribution is slightly worse than women’s characteristic distribution. The estimation of the counterfactual distributions can be done using the plug-in principle: Z Fˆhj |ki (y) = FˆYj |Xj (y|x)dFˆXk (x), y ∈ T, (2.2) where the conditional distribution estimator F̂Yj |Xj is the distribution regression estimator applied to observations of (Y, X) for group j, and the covariate distribution estimator F̂Xk is the empirical distribution of observations of X for group k. Just like other estimation problems we have seen, this estimator can be treated as the GMM estimator by simply stacking together various scores corresponding to different steps of the procedure. This also means that we can use the bootstrap as a practical inference tool. In practice, the map y 7→ F̂hj |ki (y ) may be non-monotone if y 7→ F̂Yj |Xj (y|x) is nonmonotone, but we can rearrange it into monotone function by simply sorting the values of function in a nondecreasing order. This typically improves the finite-sample properties of the estimator as we show in Section 3.2. 3. Generic Inference Method for Distributions and Quantiles 3.1. The Method. Let D denote the set of weakly increasing functions, mapping T to [0, 1]. We will call elements of this set “distribution functions”, albeit some of them need not be proper distribution functions. Let y 7→ F (y) in D denote some target distribution function. This target could be a conditional distribution function, an unconditional distribution function, or a counterfactual distribution function. Another target will be the left-inverse 6 a 7→ F ← (a) of F on T , which we can call the “quantile function” of F , where the leftinverse was defined in (1.1). We construct confidence sets for the distribution and quantile functions. The idea of the approach given here is as follows: 1. Construct confidence bands I = [L, U ] for F based upon bootstrapping an estimator F̂ . We can impose logical shape restrictions on these estimator and bands if necessary. 2. Convert these confidence bands into confidence bands for the quantile function defined as F ← by taking the inverse I ← = [U ← , L← ] of the bands I = [L, U ]. Consider a confidence band I = [L, U ] for F , with lower and upper bands L and U . Specifically, given two functions y 7→ U (y) and y 7→ L(y) in the set D such that L ≤ U , pointwise, we define a band I = [L, U ] as the collection of intervals I(y) = [L(y), U (y)], y ∈ T. We say that I covers F if F ∈ I pointwise namely F (y) ∈ I(y) for all y ∈ T . If U and L are some data-dependent bands, we say that I = [L, U ] is a confidence band for F of level p, if I covers F with probability at least p. Below we provide a bootstrap algorithm for computing a confidence band. Such method works if F̂ , the estimator of F , is sufficiently regular. Examples include estimators discussed in the previous sections. Our ultimate goal is to construct a confidence set for F ← from a generic confidence set [L, U ] for F . The following result provides a confidence set I ← . Here, we say that a collection of sets I ← = {I ← (a), a ∈ [0, 1]} covers F ← if F ← (a) ∈ I ← (a) for each a ∈ [0, 1]. Theorem 1 (Generic Bands for Quantile Functions). Consider a distribution function F and band functions L and U in the class D. (1) If F is covered by the band I := [L, U ], then the quantile function F ← is covered by the band I ← defined by I ← (a) := [U ← (a), L← (a)]. (2) Thus if the distribution function F is covered by I with probability at least p, then the quantile function F ← is covered by I ← with probability at least p. Proof. The result is immediate from the definition of the left inverse: For any a ∈ A, since L(y) ≤ F (y) for each y ∈ T , and F and L are in D, F ← (a) = inf{y ∈ T : F (y) ≥ a} ∧ sup T ≤ inf{y ∈ T : L(y) ≥ a} ∧ sup T = L← (a). 7 Analogously, conclude that F ← (a) ≥ U ← (a).  The band I ← can be narrowed without affecting coverage by exploiting the support restrictions. Suppose that T is the support of the distribution function F (and not merely a set of grid points at which we measure F ). This is relevant, for example, when outcomes are discrete or counts. Then it makes sense to exploit the support restriction that F ← (a) ∈ T by intersecting the confidence sets for F ← (a) with T . Clearly this won’t affect the coverage properties of the sets. Corollary 1 (Imposing Support Restrictions). Consider the set I˜← defined by pointwise intersection of I ← with T , namely I˜← (a) := I ← (a) ∩ T. Then, I˜← ⊆ I ← pointwise, and if I ← covers F ← then so does I˜← . The corollary is immediate because pointwise intersection of I ← with the set T does not change the coverage property, since F ← only takes values in T . Figure 1 illustrates the construction of bands using Theorem 1. It shows an F : [0, 10] 7→ [0, 1] covered by a band I = [L, U ]. It also shows that the inverse map F ← : [0, 1] 7→ [0, 10] is covered by the inverted band I ← = [U ← , L← ]. The band I ← is easy to obtain but does not exploit the fact that the support of the distribution F in this example is the set T = {0, 1, . . . , 10}. By intersecting I ← with T for each a ∈ A = [0, 1] we obtain the band I˜← which reflects such support restrictions. 3.2. Imposing Logical Shape Constraints. In many applications the point estimates F̂ and interval estimates [L0 , U 0 ] for the target distribution F do not satisfy the logical monotonicity or range restrictions, namely they don’t take values in the set D. Given such an ordered triple L0 ≤ F̂ ≤ U 0 , we can always transform it into another ordered triple L ≤ F̌ ≤ U that obey the logical monotonicity and shape restrictions. For example, we can set F̌ = S (Fˆ ), L = S(L0 ), U = S(U 0 ), (3.1) where S is the shaping operator that given a function f yields a mapping t 7→ S(f )(t) ∈ D with S(f ) = M(0 ∨ f ∧ 1), where the maximum and minimum are taken pointwise, and M is the rearrangement operator that given a function f : T 7→ [0, 1] yields a map t 7→ M(f )(t) ∈ D.1 The rearrangement operator is defined as follows. Let T be a countable set (in practice, a finite set). Given f : T 7→ [0, 1], we first consider Mf as a vector of sorted values of the set {f (t) : t ∈ T }, where the sorting is done in a non-decreasing order. Since T is an ordered 1Other monotonization operators, such as the projection on the set of weakly increasing functions, can also be used, as we remark further below. 0.0 0.4 F 0.8 8 2 4 6 8 10 6 0 2 4 InvF 8 10 y 0.0 0.2 0.4 0.6 0.8 1.0 10 p ●●●●●●●● ●●●●●●●●●● 8 ● ●●● ● 6 ●●●●●●●●●● ●● ●●● ●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● 4 InvF ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●● 2 ●● ●●●●●●●● ●●●●●●●●●● 0 ●●●●●●●●●●●●●●●●●●●●●●●●● 0.0 0.2 0.4 0.6 0.8 1.0 p Figure 1. Top: the distribution function F defined over T = {0, 1, ..., 10} and a confidence band I = [L, U ]. Both the distribution function and the confidence bands are interpolated by piecewise constant outside T . Middle: the quantile function F ← and the confidence band I ← = [U ← , L← ]. Bottom: The support-restricted confidence band I˜← = I ← ∩ T shown by circles and the inverse F ← shown by red balls; the red balls are encircled because F ← ∈ I˜← . 9 set of the same cardinality as Mf , we can assign the elements of Mf to T in one-to-one manner: to the k-th smallest element of T we assign the k-th smallest element of Mf . The resulting mapping t 7→ Mf (t) is the rearrangement map. The following lemma shows that shape restrictions improve the finite-sample properties of the estimators and confidence bands. Lemma 1 (Shaping Improves Point and Interval Estimates). The shaping operator S (a) is weakly contractive under the max distance: kS(A) − S(B)k∞ ≤ kA − Bk∞ , for any A, B: T → [0, 1], (b) is shape-neutral, S(F ) = F for any F ∈ D, (c) preserves the partial order: A ≤ B =⇒ S(A) ≤ S(B), for any A, B: T → [0, 1]. Consequently, (1) the re-shaped point estimate constructed via (3.1) is weakly better than the initial estimate under the max distance: kF̌ − F k∞ ≤ kF̂ − F k∞ , (2) the re-shaped confidence band constructed via (3.1) has weakly better coverage than the initial confidence band: P(L0 ≤ F ≤ U 0 ) ≤ P(L ≤ F ≤ U ), (3) and re-shaped confidence band constructed via (3.1) is weakly shorter than the original confidence band under the max norm, kU − Lk∞ ≤ kU 0 − L0 k∞ . The band [L, U ] is weakly better than the original band [L0 , U 0 ], in the sense that coverage is preserved while the width of the confidence band is weakly shorter. Remark 1 (Isotonization is Another Option). An alternative to rearrangement is to use the isotonization, which projects a given function on the set of weakly increasing functions that map T to [0, 1]. This also has the improving properties stated in Lemma 1. In fact any convex combination between isotonization and rearrangement has the improving properties stated in Lemma 1.  Remark 2 (Shape Restrictions on Confidence Bands by Intersection). An alternative way for imposing shape restrictions on the confidence band, is to intersect the initial band 10 [L0 , U 0 ] with the set of weakly increasing functions WI that map T to [0, 1]. That is, we simply set [LI , U I ] = WI ∩ [L0 , U 0 ] = {w ∈ WI : L0 (y) ≤ w(y) ≤ U 0 (y), ∀y ∈ T }. U0 UI Thus, is the greatest weakly increasing minorant of 0 ∨ ∧ 1 and LI is the smallest 0 majorant of 0 ∨ L ∧ 1. This approach gives the tightest confidence bands, in particular [LI , U I ] ⊆ [L, U ]. However, this construction might mis-behave under misspefication, where rearrangement continues to deliver meaningful confidence sets. Imagine that we have constructed an estimator F̂ that is consistent for the target F , which is not monotone, i.e. F 6∈ D. Then if the confidence bands [L0 , U 0 ] are sufficiently tight, then we can end up with empty intersection bands, [LI , U I ] = ∅. By contrast [L, U ] is non-empty and covers the probability limit F ∗ = S(F ) of Fˇ , where the limit target F ∗ does belongs to D.  3.3. The bootstrap algorithm. The following algorithm implements the ideas presented above. Algorithm 1 (Bootstrap Algorithm for Confidence Bands for F ). Given the ingredients above, we provide an explicit algorithm for the bootstrap construction of the joint confidence band I = [L, U ] for (F (y))y∈T based on the bootstrappable estimators (Fˆ (y))y∈T : (1) Obtain many bootstrap draws of the estimator (Fˆ (y))y∈T , (Fˆ ∗(j) (y))y∈T , j = 1, . . . , B where the index j enumerates the bootstrap draws. (2) For each y in T compute the bootstrap variance estimate ŝ2 (y) = B −1 B X (Fˆ (y)∗(j) − Fˆ (y))2 , j=1 (or use the estimate based on interquartile range). (3) Compute the critical value B  ∗(j) c(1 − α) = (1 − α)-quantile of max |F̂ (y) − Fˆ (y)|/ŝ(y) . y ∈T j=1 (4) Construct a preliminary (1 − α)-level confidence band for (F (y))y∈T as [L0 (y), U 0 (y)] = [Fˆ (y) ± c(1 − α)ŝ(y)], y ∈ T. (5) Impose the shape restrictions on Fˆ and L0 and U 0 obtaining Fˇ , and L and U . Report I = [L, U ] as a (1 − α)-level confidence band for (F (y))y∈T . 11 4. Generic Confidence Bands for Quantile Effects Our next goal is to construct a confidence band for the quantile effect function a 7→ ∆(a) defined by ∆(a) := F1← (a) − F0← (a), corresponding to the difference of quantile functions of the two distribution functions F0 and F1 with support sets T . The basic idea is as follows: 1. Construct joint confidence bands for F0 and F1 . 2. Convert these into confidence bands for F0← and F1← by inversion. 3. Contruct the confidence band for the quantile effect ∆ by taking the pointwise Minkowski difference between the confidence bands for F1← and F0← . Specifically, suppose we have the confidence bands I0← for F0← and I1← for F1← , which jointly cover F0← and F1← with probability p. For example, we can construct these sets using Theorem 1 in conjunction with the Bonferroni inequality.2 Algorithm 2 provides a construction of the confidence bands that has joint coverage property and is less conservative than using Bonferroni. Then we can convert these bands to confidence bands for ∆ by taking the pointwise Minkowski difference of the pairs of intervals at each a treated as sets of points. Recall that the Minkowski difference between two subsets V and U of a vector space is V U := {v − u : v ∈ V, u ∈ U }. Note that if V and U are intervals [v1 , v2 ] and [u1 , u2 ], then V U = [v1 , v2 ] [u1 , u2 ] = [v1 − u2 , v2 − u1 ]. Theorem 2 (Generic Bands for Quantile Effect Functions). Consider the distribution functions F0 and F1 and the band functions L0 , U0 , L1 and U1 in the class D. (1) If Fk is covered by Ik := [Lk , Uk ] for k = 0 and k = 1, then the quantile effect function ← = [U ← , L← ] − [U ← , L← ], where the minus ∆ = F1← − F0← is covered by the band I∆ 1 1 0 0 operator is defined by a pointwise Minkowski difference: ← I∆ (a) := [U1← (a), L← 1 (a)] [U0← (a), L← 0 (a)]. 2The joint coverage of two confidence sets with marginal coverage probabilities p̃ is at least p = 2p̃ − 1 by Bonferroni inequality. 12 (2) If the distribution functions F0 and F1 are jointly covered by I0 and I1 with probability at ← with probability least p, then the quantile effect function ∆ = F1← − F0← is covered by I∆ at least p. Proof. The results is immediate from the definition of the Minkowski sum.  ← without affecting coverage by imposing As in Theorem 1, we can narrow the band I∆ support restrictions on the bands for the quantile functions. ← = I˜← − I˜← defined by: Corollary 2 (Imposing Support Restrictions). Consider the band I˜∆ 1 0 ← I˜∆ (a) := I˜1← (a) I˜0← (a), I˜k← (a) := {[Uk← (a), L← k (a)] ∩ T }, k ∈ {0, 1}. ← ⊆ I ← , and if I ← covers ∆ then so does I˜← . Then I˜∆ ∆ ∆ ∆ Figure 2 shows bands for a pair of quantile functions together with bands for the quantile effect function constructed using Theorem 2. The top plot shows the bands I0← and I1← for the quantile functions F0← and F1← . The middle plot shows the band I∆ for the quantile effect function ∆ = F1← − F0← , obtained by taking the Minkowski difference of I1← and I0← . The bottom plot shows the confidence band I˜∆ for the quantile effect function ∆ resulting from imposing the support constraints. As the Theorem 2 predicts, the quantile function ∆ is covered by the band I∆ . In what follows we write down an explicit algorithm that implements the proposal above. Algorithm 2 (Bootstrap Algorithm for Confidence Bands for Quantile Effects). (1) Obtain many bootstrap draws of the estimator {(Fˆk (y))y∈T }k∈{0,1} , ∗(j) {(Fˆk (y))y∈T }k∈{0,1} , j = 1, . . . , B where the index j enumerates the bootstrap draws. (2) For each y in T and k ∈ {0, 1} compute the bootstrap variance estimate ŝk2 (y) = B −1 B X ∗(j) (Fˆ (y)k − Fˆk (y))2 , j=1 (or use the estimate based on the interquartile range). (3) Compute the critical value  B ∗(j ) max |Fˆk (y) − Fˆk (y)|/ŝk (y) c(1 − α) = (1 − α)-quantile of . y∈T,k∈{0,1} j=1 6 4 0 2 Quantiles 8 10 13 0.0 0.2 0.4 0.6 0.8 1.0 0.6 0.8 1.0 8 6 4 0 2 Quantile Effect 10 p 0.0 0.2 0.4 8 ●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●● 6 ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● 4 ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● 2 Quantile Effect 10 p 0 ●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●● 0.0 0.2 0.4 0.6 0.8 1.0 p Figure 2. Top: the two quantiles functions F0← and F1← and confidence bands I0← and I1← . Middle: the quantile effect F1← − F0← and the confidence ← . Bottom: the support restricted confidence band I˜← . band I∆ ∆ 14 (4) Construct a preliminary (1 − α)-level joint confidence region for {(Fk (y))y∈T }k∈{0,1} as [L0k (y), Uk0 (y)] = [Fˆk (y) ± c(1 − α)ŝk (y)], y ∈ T, k ∈ {0, 1}. (5) Impose the shape restrictions on F̂k and L0k and Uk0 by setting: Fˇk = S(Fˆk ), [Lk , Uk ] = [S(L0k ), S(Uk0 )]. (6) Report Ik = [Lk , Uk ] for k ∈ {0, 1} as a (1 − α)-level joint confidence band for (Fk (y))y∈T for k ∈ {0, 1}. Report Ik← = [Uk← , L← k ] for k ∈ {0, 1} as a (1−α)-level joint confidence band for (Fk← (y))y∈T for k ∈ {0, 1}. (7) Report ← ← ← I∆ = [U1← , L← 1 ] − [U0 , L0 ] as a (1 − α)-level confidence band for quantile effect. (8) If necessary, impose support restrictions in the previous two steps, following Corollaries 1 and 2. 5. Gender Wage Gap in 2012 To illustrate the use of the bands in practice we consider an application to gender wage gap using data from the U.S. March Supplement of the Current Population Survey (CPU) in 2012. We select white non-hispanic individuals, aged 25 to 64 years, and working more than 35 hours per week during at least 50 weeks of the year. We exclude self-employed workers; individuals living in group quarters; individuals in the military, agricultural or private household sectors; individuals with inconsistent reports on earnings and employment status; and individuals with allocated or missing information in any of the variables used in the analysis. The resulting sample consists of 29, 217 workers including 16, 690 men and 12, 527 of women. The variable of interest Y is the logarithm of the hourly wage rate constructed as the ratio of the annual earnings to the total number of hours worked, which is constructed in turn as the product of number of weeks worked and the usual number of hours worked per week.3 In this application F0 and F1 correspond to observed and counterfactual distributions of wages for women and men. Following Section 2, we denote these distributions by Fhj|ki , where Z Fhj |ki (y) = FYj |Xj (y|x)dFXk (x), y ∈ T, for j, k ∈ {m, w}, where m and w refer to men and women, FYj |Xj is the wage structure in group j, and FXk is the distribution of worker characteristics in group k. The worker 3This sample selection criteria and the variable construction follow [9]. 15 characteristics X include 5 marital status indicators (widowed, divorced, separated, never married, and married); 6 educational attainment indicators (0-8 years of schooling completed, high school dropouts, high school graduates, some college, college graduate, and advanced degree); 4 region indicators (midwest, south, west, and northeast); and a quartic in potential experience constructed as the maximum of age minus years of schooling minus 7 and zero, i.e., experience = max(age − education − 7, 0), interacted with the educational attainment indicators. We estimate Fhw|wi and Fhm|mi using the empirical distributions of Y for women and men. We estimate Fhm|wi by (2.2), where FˆYm |Xm is the distribution regression estimator with a logit link and a linear index in X in the sample of men. We use the empirical distribution of X for women to estimate FXw . All the estimators use sampling weights to account for nonrandom sampling in the March CPS. Table 1. Descriptive Statistics All Men Women log wage 2.79 2.90 female 0.43 0.00 married 0.66 0.69 widowed 0.01 0.00 divorced 0.12 0.10 separated 0.02 0.02 never married 0.19 0.19 0-8 years completed 0.00 0.01 high school dropout 0.02 0.03 high school graduated 0.25 0.27 some college 0.28 0.27 college graduated 0.28 0.28 advanced degree 0.15 0.14 northeast 0.20 0.20 midwest 0.27 0.27 south 0.35 0.35 west 0.18 0.19 potential experience 18.96 19.01 2.65 1.00 0.63 0.02 0.15 0.02 0.18 0.00 0.02 0.23 0.30 0.29 0.17 0.19 0.28 0.35 0.18 18.90 Source: March Supplement CPS 2012 Table 1 reports descriptive statistics for the variables used in the analysis. Working women are more highly educated than working men, but have less experience. The unconditional average gender wage gap is 25%. Figure 3 shows that the log hourly wage variable 16 Y has multiple mass points due to rounding in the reporting of the earnings and labor supply variables. In particular, we observe 2, 633 different values in the 16, 690 observations for men and 2, 091 different values in the 12, 527 observations for women. We note here that the inference methods of Section 3 are fully robust to the process that generates mass points in the variable Y . 500 60 Frequency 1000 1500 80 100 120 140 Men log hourly wages 0 0 20 40 Frequency Women log hourly wages −6 −2 2 4 Log of hourly wage 6 −8 −4 0 2 4 6 Log of hourly wage Figure 3. Histograms of log hourly wages for women and men. Figures 4 and 5 show point estimates and 90% joint confidence bands for the observed and counterfactual distribution and quantile functions of wages. We construct the confidence bands for the distributions using Algorithm 2 by empirical bootstrap with 200 repetitions. Following Theorem 1, we apply left-inverse to these bands to obtain the confidence bands for the quantile functions. We do not report bands that impose the support constrains because the variable Y takes on many values. Here we find that the wage distribution for men first order stochastically dominates the wage distribution for women. ← The counterfactual quantile function Fhm |wi is very similar to the observed quantile function for men, suggesting that most of the wage gap is due to gender discrimination. This observation is confirmed in Figure 6, which plots point estimates and 90% simultaneous confidence bands for the decomposition of the quantile gender wage gap into composition and discrimination effects of (2.1). We follow Theorem 2 to construct the bands by the Minkowski difference of the bands for the quantiles, as described in Algorithm 2. The gender wage gap ranges from 20% to 30% and is increasing with the quantile index. Most of the gap is explained by gender wage discrimination with the composition having a small 17 negative effect non statistically significant for most quantiles. This low contribution of the composition effect is explained by the similarity in the observable characteristics of man and women in Table 1. 6. Racial differences in mental ability of young children As a second empirical application, we analyze the racial IQ test score gap. We use data from the US Collaborative Perinatal Project (CPP) obtained from [7]. These data contain information on children from 30,002 women who gave birth in 12 medical centers between 1959 and 1965. Our main outcome of interest, Y , is the standardized test scores at the age of seven years (both Stanford-Binet and Wechsler Intelligence Test). In addition to the test score measure, the dataset contains a rich set of background characteristics for the children, X, including information on age, gender, region, socioeconomic status, home environment, prenatal conditions, and interviewer fixed effects. [7] provide a comprehensive description of the dataset and extensive descriptive statistics. A key feature of the test score variable is the discrete nature of its distribution. We observe only 128 different values for the standardized test score. Figure 7 presents the corresponding histogram. Note that each bar corresponds to exactly one value. For instance, more than 4% of the observations have exactly the same score. This is a common feature of test scores, which are necessarily discrete because they are based on a finite number of questions. In this application F0 and F1 correspond to observed and counterfactual distributions of test scores for black and white children. Following Section 2, we denote these distributions by Fhj|ki , where Z Fhj |ki (y) = FYj |Xj (y|x)dFXk (x), y ∈ T, (6.1) for j, k ∈ {w, b}, where w and b refer to white and black, FYj |Xj is the conditional distribution of test scores in group j, and FXk is the distribution of background characteristics in group k. With these counterfactual test score distributions it is possible to decompose the observed black-white test score gap into ← ← ← ← ← ← Fhw|w i − Fhb|bi = [Fhw|wi − Fhw|bi ] + [Fhw|bi − Fhb|bi ], where the first term in brackets corresponds is the composition effect due to differences in observable background characteristics and the second term is the unexplained difference. We estimate Fhw|wi and Fhb|bi by the empirical test score distributions for white and black children, respectively. We estimate the counterfactual distribution Fhw|bi by the sample analog of (6.1) replacing FYw |Xw by the DR estimator for white children, and FXb by the empirical distribution of X for black children. We use the logistic link function for the DR, but the results using the linear link function or the normal link function are similar. 18 1.0 Observed and Counterfactual Distributions (90% CI) 0.6 0.4 0.0 0.2 Probability 0.8 Observed women distribution Observed men distribution Counterfactual distribution 1.5 2.0 2.5 3.0 3.5 4.0 Log of hourly wage Figure 4. The distribution functions of observed wages for women and men, and the distribution function of counterfactual wages for women under men’s wage structure, with 90% joint confidence bands. Confidence bands obtained by empirical bootstrap with 200 repetitions. 19 Observed women quantiles Observed men quantiles Counterfactual quantiles 3.0 2.5 1.5 2.0 Log of hourly wage 3.5 4.0 Observed and Counterfactual Quantiles (90% CI) 0.2 0.4 0.6 0.8 Quantile index Figure 5. The quantile functions of observed wages for women and men, and the quantile function of counterfactual wages for women under men’s wage structure, with 90% joint confidence bands. 20 0.5 0.1 0.3 QE −0.1 Difference in log of hourly wage Gender wage gap (90% CI) 0.2 0.4 0.6 0.8 Quantile index 0.5 0.1 0.3 QE −0.1 Difference in log of hourly wage Discrimination (90% CI) 0.2 0.4 0.6 0.8 Quantile index 0.5 0.1 0.3 QE −0.1 Difference in log of hourly wage Composition (90% CI) 0.2 0.4 0.6 0.8 Quantile index Figure 6. Decomposition of observed gender wage gap into composition and discrimination effects, with 90% confidence bands. The composition effect is the difference between the quantile functions of observed wages for men and counterfactual wages for women under the men’s wage structure. The discrimination effect is the difference between the quantile functions of counterfactual wages for women under the men’s wage structure and observed wages for women. 21 0.03 0.02 0.00 0.01 Fraction 0.04 0.05 Test scores for seven−year−olds −3 −2 −1 0 1 2 3 Test score Figure 7. Histogram of test scores of seven year old children. Each bind corresponds to a unique value of the test score. Figure 8 reports the results of the decomposition. The first panel shows the observed ← ← ← and counterfactual quantile functions, Fhw|w i , Fhb|bi and Fhw|bi . The second panel shows the ← ← difference between the observed quantile functions, Fhw|w i − Fhb|bi . The third and fourth ← panels decompose these observed differences into the composition effect (Fh← w|wi − Fhw|bi ) ← ← and the unexplained component (Fhw|bi − Fhb|bi ). The point estimates are shown with their respective 95% simultaneous confidence bands constructed by Algorithm 2 using weighted bootstrap with standard exponential weights and B = 1, 000 replications. The bands impose the restrictions that the supports of the test scores correspond to the observed values in the sample. We find a large and statistically significant positive raw black-white gap. A formal test based on the uniform bands rejects the null hypothesis of a zero or a negative racial test score gap at all quantiles. The estimated QE function is increasing in the quantile index ranging from below 0.6 standard deviation units at the lower tail up to over one standard deviation unit at the upper tail of the distribution. The quantile differences at the tails 22 Quantile functions Black white observed gap 1.0 0.5 0.0 0 −2 −1 Test scores 1 Difference in test scores 2 Observed black quantiles Observed white quantiles Counterfact. quantiles 0.0 0.2 0.4 0.6 0.8 1.0 0.0 Composition effect 0.4 0.6 0.8 1.0 0.2 0.4 0.6 Probability 0.8 1.0 0.5 1.0 0.0 0.0 0.5 Difference in test scores 1.0 Unexplained difference 0.0 Difference in test scores 0.2 0.0 0.2 0.4 0.6 0.8 1.0 Probability Figure 8. Decomposition of observed racial differences in mental ability of seven year old children. Quantile functions, raw difference, composition effect, and unexplained difference including support restricted 95% confidence bands. substantially differ from the mean difference of 0.85 standard deviation units reported in [7]. In fact, we can formally reject the null hypothesis of a constant raw test score gap across the distribution because we can not draw a horizontal line at any value of the difference 23 of test scores, which is covered by the confidence band of the QE function at all quantile indexes. Our decomposition analysis shows that about two third of this gap can be explained by differences in the distribution of observable characteristics. Nevertheless, the remaining unexplained difference is significant, both in economic and in statistical terms. Looking at the quantile effect function, we can see that there is substantial effect heterogeneneity along the distribution. Interestingly, the increase in the test score gap at the upper quantiles can be fully explained by differences in background characteristics between black and white children. The resulting unexplained difference is maximized in the center of the distribution. Finally, our simultaneous confidence bands allow for testing several interesting hypothesis’ about the whole quantile effect function. For instance, we can reject the null hypothesis that the composition effect and the unexplained difference are zero, negative, or constant at all quantiles but we cannot reject that they are positive everywhere. Notes The distribution regression model was developed in [11], [5], and [3]. Ronald Oaxaca [10] and Alan Blinder [1] pioneered the use of least squares methods to carry out decompositions of mean wages. The decomposition was extended to distributions in [4], [8], [6], and [3], among others. The generic inference method of Section 3 was developed in [2]. Appendix A. Problems (1) Demonstrate that estimation of the counterfactual distribution can be formulated as a GMM framework. (2) Replicate the results of the empirical section. Provide a brief explanation for what you are doing. You can use the R code provided, but you need to write comments in the code explaining what each block of the code is doing. References [1] Alan S. Blinder. Wage discrimination: Reduced form and structural estimates. The Journal of Human Resources, 8(4):436–455, 1973. [2] V. Chernozhukov, I. FernandezVal, B. Melly, and K. Wut ´ ¨ hrich. Generic Inference on Quantile and Quantile Effect Functions for Discrete Outcomes. ArXiv e-prints, August 2016. [3] Victor Chernozhukov, Ivan Val, and Blaise Melly. Inference on counterfac´ Fernandez´ tual distributions. Econometrica, 81(6):2205–2268, 2013. 24 [4] John DiNardo, Nicole M. Fortin, and Thomas Lemieux. Labor market institutions and the distribution of wages, 1973-1992: A semiparametric approach. Econometrica, 64(5):1001–1044, 1996. [5] Silverio Foresi and Franco Peracchi. The conditional distribution of excess returns: An empirical analysis. Journal of the American Statistical Association, 90(430):451–466, 1995. [6] Nicole Fortin, Thomas Lemieux, and Sergio Firpo. Decomposition Methods in Economics, volume 4 of Handbook of Labor Economics, chapter 1, pages 1–102. Elsevier, 2011. [7] Jr. Fryer, Roland G. and Steven D. Levitt. Testing for racial differences in the mental ability of young children. American Economic Review, 103(2):981–1005, April 2013. [8] José A. F. Machado and José Mata. Counterfactual decomposition of changes in wage distributions using quantile regression. J. Appl. Econometrics, 20(4):445–465, 2005. [9] Casey B. Mulligan and Yona Rubinstein. Selection, investment, and women’s relative wages over time. The Quarterly Journal of Economics, 123(3):1061–1110, 2008. [10] Ronald Oaxaca. Male-female wage differentials in urban labor markets. International Economic Review, 14(3):693–709, 1973. [11] O Dale Williams and James E Grizzle. Analysis of contingency tables having ordered response categories. Journal of the American Statistical Association, 67(337):55–63, 1972.