Network and Panel Quantile Effects Via Distribution Regression

Network and Panel Quantile Effects Via Distribution Regression

{\bf Keywords:} Quantile Effects, Counterfactual Distributions, Fixed Effects, Incidental Parameter Problem, Long Panels

Introduction

Standard regression analyzes average effects of covariates on outcome variables. In many applications it is equally important to consider distributional effects. For example, a policy maker might be interested in the effect of an education reform not only on the mean but also the entire distribution of test scores or wages. Availability of panel data is very useful to identify ceteris paribus average and distributional effects because it allows the researcher to control for multiple sources of unobserved heterogeneity that might cause endogeneity or omitted variable problems. The idea is to use variation of the covariates over time for each individual or over individuals for each time period to account for unobserved individual and time effects. In this paper we develop inference methods for distributional effects in nonlinear models with two-way unobserved effects. They apply not only to traditional panel data models where the unobserved effects correspond to individual and time fixed effects, but also to models for other types of data where the unobserved effects reflect some grouping structure such as unobserved sender and receiver effects in network data models. The unobserved effects will be treated as fixed effects, i.e. parameters to be estimated, leaving their relation to observed covariates unrestricted.

We develop inference methods for quantile functions and effects. The quantile function corresponds to the marginal distribution of the outcome in a counterfactual scenario where the treatment covariate of interest is set exogenously at a desired level and the rest of the covariates and unobserved effects are held fixed, extending the construction of Chernozhukov Fernandez-Val and Melly ChernozhukovFernandezValMelly2013 for the cross section case. The quantile effect is the difference of quantile functions at two different treatment levels. Our methods apply to continuous and discrete treatments by appropriate choice of the treatment levels, and have causal interpretation under standard unconfoundedness assumptions for panel data. The inference is based upon the generic method of Chernozhukov, Fernandez-Val, Melly, and Wuthrich CFVMW2015 that projects joint confidence bands for distributions into joint confidence bands for quantile functions and effects. This method has the appealing feature that applies without modification to any type of outcome, let it be continuous, discrete or mixed.

The key input for the inference method is a joint confidence band for the counterfactual distributions at the treatment levels of interest. We construct this band from fixed effects distribution regression (FE-DR) estimators of the conditional distribution of the outcome given the observed covariates and unobserved effects. In doing so, we extend the distribution regression approach to model conditional distributions with unobserved effects. This version of the DR model is semiparametric because not only the DR coefficients can vary with the level of the outcome as in the cross section case, but also the distribution of the unobserved effects is left unspecified. We show that the FE-DR estimator can be obtained as a sequence of binary response fixed effects estimators where the binary response is an indicator of the outcome passing some threshold. To deal with the incidental parameter problem associated with the estimation of the unobserved effects (Neyman and Scott NeymanScott1948), we extend the analytical bias corrections of FernandezValWeidner2016 for single binary response estimators to multiple (possibly a continuum) of binary response estimators. In particular, the main technical contribution is to establish functional central limit theorems for the fixed effects estimators of the DR coefficients and associated counterfactual distributions, and show the validity of the bias corrections under asymptotic sequences where the two dimensions of the data set pass to infinity at the same rate. As in the single binary response model, the bias corrections remove the asymptotic bias of the fixed effects estimators without increasing their asymptotic variances.

We implement the inference method using multiplier bootstrap gz-84. This version of bootstrap constructs draws of an estimator as weighted averages of its influence function, where the weights are independent from the data. Compared to empirical bootstrap, multiplier bootstrap has the computational advantage that it does not involve any parameter reestimation. This advantage is particularly convenient in our setting because the parameter estimation require multiple nonlinear optimizations that can be highly dimensional due to the fixed effects. Multiplier bootstrap is also convenient to account for data dependencies. In network data, for example, it might be important to account for reciprocity or pairwise clustering. Reciprocity arises because observational units corresponding to the same pair of agents but reversing their roles as sender and receiver might be dependent even after conditioning on the unobserved effects. By setting the weights of these observational units equal, we account for this dependence in the multiplier bootstrap. In addition to the previous practical reasons, there are some theoretical reasons for choosing multiplier bootstrap. Thus, cck-16 established bootstrap functional central limit theorems for multiplier bootstrap in high dimensional settings that cover the network and panel models that we consider.

The methods developed in this paper apply to models that include unobserved effects to capture grouping or clustering structures in the data such as models for panel and network data. These effects allow us to control for unobserved group heterogeneity that might be related to the covariates causing endogeneity or omitted variable bias. They also serve to parsimoniously account for dependencies in the data. We illustrate the wide applicability with an empirical example to gravity models of trade. In this case the outcome is the volume of trade between two countries and each observational unit corresponds to a country pair indexed by exporter country (sender) and importer country (receiver). We estimate the distributional effects of gravity variables such as the geographical distance controlling for exporter and importer country effects that pick up unobserved heterogeneity possibly correlated with the gravity variables. We uncover significant heterogeneity in the effects of distance and other gravity variables across the distribution, which is missed by traditional mean methods. We also find that the Poisson model, which is commonly used in the trade literature to deal with zero trade in many country pairs, does not provide a good approximation to the distribution of the volume of trade due to heavy tails.

\paragraph{Literature review} Unlike mean effects, there are different ways to define distributional and quantile effects. For example, we can distinguish conditional effects versus unconditional or marginalized effects, or quantile effects versus quantiles of the effects. Here we give a brief review of the recent literature on distributional and quantile effects in panel data models emphasizing the following aspects: (1) type of effect considered; (2) type of unobserved effects in the model; and (3) asymptotic approximation. For the unobserved effects, we distinguish models with one-way effects versus two-way effects. For the asymptotic approximation we distinguish short panels with large $N$ and fixed $T$ versus long panels with large $N$ and large $T$ , where $N$ and $T$ denote the dimensions of the panel. We focus mainly on fixed effects approaches where the unobserved effects are treated as parameters to be estimated, but also mention some correlated random effects approaches that impose restrictions on the distribution of the unobserved effects. This paper deals with inference on marginalized quantile effects in large panels with two-way effects, which has not been previously considered in the literature. Indeed, to the best of our knowledge, it is the first paper to provide inference methods for quantile treatment effects from panel and network models with two-way fixed effects.

Koenker Koenker2004 introduced fixed effects quantile regression estimators of conditional quantile effects in large panel models with one-way individual effects using shrinkage to control the variability in the estimation of the unobserved effects. Lamarche2010 discussed the optimal choice of a tuning parameter in Koenker's method. In the same framework, Kato, Galvao, and Montes-Rojas KatoGalvaoMontesRojas2012, Galvao, Lamarche and Lima GalvaoLamarcheLima2013, Galvao and Kato GalvaoKato2016 and Arellano and Weidner ArellanoWeidner2016 considered fixed effects quantile regression estimators without shrinkage and developed bias corrections. All these papers require that $T$ pass to infinity faster than $N$, making it difficult to extend the theory to models with two-way individual and time effects. Graham, Hahn and Powell GrahamHahnPowell2009 found a special case where the fixed effects quantile regression estimator does not suffer of incidental parameter problem. Machado and Santos Silva mss18 has recently proposed a method to estimate conditional quantile effects in a location-scale model via moments.

In short panels, Rosen Rosen2012 showed that a linear quantile restriction is not sufficient to point identify conditional effects in a panel linear quantile regression model with unobserved individual effects. Chernozhukov, Fernandez-Val, Hahn and Newey ChernozhukovFernandezValHahnNewey2013 and chernozhukov2015nonparametric discussed identification and estimation of marginalized quantile effects in nonseparable panel models with unobserved individual effects and location and scale time effects under a time homogeneity assumption. They showed that the effects are point identified only for some subpopulations and characterized these subpopulations. Graham, Hahn, Poirier and Powell GrahamHahnPoirierPowell2015 considered quantiles of effects in linear quantile regression models with two-way effects. Finally, Abrevaya and Dahl AbrevayaDahl2008 and Arellano and Bonhomme ArellanoBonhomme2016 developed estimators for conditional quantile effects in linear quantile regression model with unobserved individual effects using correlated random effects approaches. None of the previous quantile regression based methods apply to discrete outcomes.

Finally, we review previous applications of panel data methods to network data. These include candelaria16, charbonneau17, CFW-16, Dzemski2017, FernandezValWeidner2016, gao20, Graham2016, Graham2017, jochmans18, toth17, and Yan2016statistical, which developed methods for models of network formation with unobserved sender and receiver effects for directed and undirected networks.\footnote{We refer to paula19 for an excellent up to date review on this topic.} None of these papers consider estimation of quantile effects as the outcome variable is binary, whether or not a link is formed between two agents.

\paragraph{Plan of the paper} Section (ref) introduces the distribution regression model with unobserved effects for network and panel data, and describes the quantities of interest including model parameters, distributions, quantiles and quantile effects. Section (ref) discusses fixed effects estimation, bias corrections to deal with the incidental parameter problem, and uniform inference methods. Section (ref) provides asymptotic theory for the fixed effects estimators, bias corrections, and multiplier bootstrap. Section (ref) and (ref) report results of the empirical application to the gravity models of trade and a Monte Carlo simulation calibrated to the application, respectively. The proofs of the main results are given in the Appendix, and additional technical results are provided in the Supplementary Appendix.

\paragraph{Notation} For any two real numbers $a$ and $b$, $a\vee b = \max\{a,b\}$ and $a\wedge b = \min\{a,b\}$. For a real number $a$, $\lfloor a \rfloor$ denotes the integer part of $a$. For a set $\mathcal{A}$, $| \mathcal{A}|$ denotes the cardinality or number of elements of $\mathcal{A}$.

Model and Parameters of Interest

Distribution Regression Model with Unobserved Effects

We observe the data set $\{(y_{ij}, x_{ij}) : (i,j) \in \mathcal{D} \}$, where $y_{ij}$ is a scalar outcome variable with region of interest $\mathcal{Y}$, and $x_{ij}$ is a vector of covariates with support $\mathcal{X} \subseteq \mathbb{R}^{d_x}$.\footnote{If $y_{ij}$ has unbounded support, then the region $\mathcal{Y}$ is usually a subset of the support to avoid tail estimation.} The variable $y_{ij}$ can be discrete, continuous or mixed. The subscripts $i$ and $j$ index individuals and time periods in traditional panels, but they might index other dimensions in more general data structures. In our empirical application, for example, we use a panel where $y_{ij}$ is the volume of trade between country $i$ and country $j$, and $x_{ij}$ includes gravity variables such as the distance between country $i$ and country $j$. Both $i$ and $j$ index countries as exporters and importers respectively. The set $\mathcal{D}$ contains the indexes of the pairs $(i,j)$ that are observed. It is a subset of the set of all possible pairs $\mathcal{D}_0 := \{(i,j) : i = 1,\dots,I; j = 1, \dots, J \}$, where $I$ and $J$ are the dimensions of the panel. We introduce $\mathcal{D}$ to allow for certain forms of missing data that are common in panel and network applications, see Assumption (ref)(v) in Section (ref). For example, in the trade application $I=J$ and $\mathcal{D} = \mathcal{D}_0 \setminus \{(i,i) : i = 1,\dots,I \}$ because we do not observe trade of a country with itself. We denote the total number of observed units by $n$, i.e. $n = |\mathcal{D}|$.

Let $v_i$ and $w_j$ denote vectors of unspecified dimension that contain unobserved random variables or effects that might be related to the covariates $x_{ij}$. In traditional panels, $v_i$ are individual effects that capture unobserved individual heterogeneity and $w_j$ are time effects that account for aggregate shocks. More generally, these variables serve to capture some forms of endogeneity and group dependencies in a parsimonious fashion. We specify the conditional distribution of $y_{ij}$ given $(x_{ij},v_i,w_{j})$ using the distribution regression (DR) model with unobserved effects

where $\Lambda_y$ is a known link function such as the normal or logistic distribution, which may vary with $y$, $x \mapsto P(x)$ is a dictionary of transformations of $x$ such us polynomials, b-splines and tensor products, $\beta(y)$ is an unknown parameter vector, which can vary with $y$, and $(v,y) \mapsto \alpha(v,y)$ and $(w,y) \mapsto \gamma(w,y)$ are unspecified measurable functions. This DR model is a semiparametric model for the conditional distribution because $y \mapsto \theta(y) := (\beta(y),\alpha(v_1,y), \ldots, \alpha(v_I,y), \gamma(w_1,y), \ldots, \gamma(w_J,y))$ is a function-valued parameter and the dimension of $\theta(y)$ varies with $I$ and $J$, although we do not make this dependence explicit. We shall treat the dimension of $P(x)$ as fixed and set $\Lambda_y$ equal to the logistic distribution for all $y$ in the asymptotic analysis.

When $y_{ij}$ is continuous, the model (ref) has the following representation as an implicit nonseparable model by the probability integral transform $$ \Lambda_{y_{ij}}(P(x_{ij})'\beta(y_{ij}) + \alpha(v_i,y_{ij}) + \gamma(w_{j}, y_{ij})) = u_{ij}, \ \ u_{ij} \mid x_{ij}, v_i, w_{j} \sim U(0,1), $$ where the error $u_{ij}$ represents the unobserved ranking of the observation $y_{ij}$ in the conditional distribution. The parameters of the model are related to derivatives of the conditional quantiles. Let $Q_{y_{ij}}(u \mid x_{ij}, v_{i}, w_{j})$ be the $u$-quantile of $y_{ij}$ conditional on $(x_{ij}, v_{i}, w_{j})$ defined as the left-inverse of $y \mapsto F_{y_{ij} }(y \mid x_{ij}, v_i, w_{j})$ at $u$, namely $$ Q_{y_{ij} }(u \mid x_{ij}, v_{i}, w_{j}) = \inf\{y \in \mathcal{Y} : F_{y_{ij} }(y \mid x_{ij}, v_i, w_{j}) \geq u\} \wedge \sup \{ y \in {\mathcal{Y}}\}, $$ and $x_{ij} = (x_{ij}^1,\ldots, x_{ij}^{d_x})$.\footnote{We use the convention $\inf \{\emptyset\} = + \infty$.} Then, it can be shown that if $y \mapsto F_{y_{ij} }(y \mid x_{ij}, v_i, w_{j})$ is strictly increasing in the support of $y_{ij}$, $\partial \Lambda_y(z)/\partial z > 0$ for all $y$ in the support of $y_{ij}$ and $x_{ij} \mapsto Q_{y_{ij} }(u \mid x_{ij}, v_{i}, w_{j})$ is differentiable,\footnote{Indeed, $\Lambda_y(P(x_{ij})'\beta(y) + \alpha(v_i,y) + \gamma(w_{j}, y)) =u$ at $y = Q_{y_{ij}}(u \mid x_{ij}, v_{i}, w_{j})$. Differencing this expression with respect to $x_{ij}^k$ yields $$\left. \partial_{x_{ij}^{k}} P(x_{ij})'\beta(y) \right|_{y = Q_{y_{ij}}(u \mid x_{ij}, v_{i}, w_{j})} = - \left. \frac{\partial \Lambda_y\left(P(x_{ij})'\beta(y) + \alpha(v_i,y) + \gamma(w_{j}, y) \right)/\partial y}{\lambda_y\left(P(x_{ij})'\beta(y) + \alpha(v_i,y) + \gamma(w_{j}, y)\right)} \right|_{y = Q_{y_{ij}}(u \mid x_{ij}, v_{i}, w_{j})} \partial_{x_{ij}^{k}} Q_{y_{ij} }(u \mid x_{ij}, v_{i}, w_{j}),$$ where $\lambda_y(z) = \partial \Lambda_y(z)/\partial z$. Note that the first term of the right hand side does not depend on $k$ and is positive because $y \mapsto F_{y_{ij} }(y \mid x_{ij}, v_i, w_{j}) = \Lambda_y(P(x_{ij})'\beta(y) + \alpha(v_i,y) + \gamma(w_{j}, y))$ is strictly increasing at $y = Q_{y_{ij}}(u \mid x_{ij}, v_{i}, w_{j})$.} $$ \left. \partial_{x_{ij}^{k}} P(x_{ij})'\beta(y) \right|_{y=Q_{y_{ij} }(u \mid x_{ij}, v_{i}, w_{j})} \propto - \partial_{x_{ij}^{k}} Q_{y_{ij} }(u \mid x_{ij}, v_{i}, w_{j}), \ \ k = 1,\ldots,d_x, \ \ \partial_{x_{ij}^{k}} := \partial/ \partial x_{ij}^k. $$ If $P(x_{ij}) = x_{ij}$, then $\partial_{x_{ij}^{k}} P(x_{ij})'\beta(y) = \beta_{k}(y)$ such that $$ \left. \frac{\beta_{\ell}(y)}{\beta_k(y)} \right|_{y = Q_{y_{ij}}(u \mid x_{ij}, v_{i}, w_{j})} = \frac{\partial_{x_{ij}^{\ell}} Q_{y_{ij} }(u \mid x_{ij}, v_{i}, w_{j})}{\partial_{x_{ij}^k} Q_{y_{ij} }(u \mid x_{ij}, v_{i}, w_{j})}, \ \ \ell,k = 1,\ldots,d_x, $$ provided that $\partial_{x_{ij}^k} Q_{y_{ij} }(u \mid x_{ij}, v_{i}, w_{j}) \neq 0.$ The DR coefficients therefore are proportional to (minus) derivatives of the conditional quantile function, and ratios of DR coefficients correspond to ratios of derivatives.

Estimands

In addition to the model parameter $\beta(y)$, we are interested in measuring the effect on the outcome of changing one of the covariates holding the rest of the covariates and the unobserved effects fixed. Let $x = (t, z')'$, where $t$ is the covariate of interest or treatment and $z$ are the rest of the covariates that usually play the role of controls. One effect of interest is the quantile (left-inverse) function (QF) $$ Q_k(\tau) = F_k^{\leftarrow}(\tau) := \inf \{y \in \mathcal{Y} : F_k(y) \geq \tau\} \wedge \sup \{ y \in {\mathcal{Y}}\}, \ \ \tau \in (0,1), $$ where $$ F_k(y) = n^{-1} \sum_{(i,j) \in \mathcal{D}} \Lambda_y(P(t_{ij}^k,z_{ij}')'\beta(y) + \alpha(v_i,y) + \gamma(w_j,y)), $$ $t_{ij}^k$ is a level of the treatment that may depend on $t_{ij}$, and $k \in \{0,1\}$. We provide examples below. Note that in the construction of the counterfactual distribution $F_k$, we marginalize $(x_{ij},v_i,w_j)$ using the empirical distribution. The resulting effects are finite population effects. We shall focus on these effects because conditioning on the covariates and unobserved effects is natural in the trade application.\footnote{The distinction between finite and infinite population effects does not affect estimation, but affects inference AAIW-14. The estimators of infinite population effects need to account for the additional sampling variation coming from the estimation of the distribution of $(x_{ij},v_i,w_j)$.} We construct the quantile effect function (QEF) by taking differences of the QF at two treatment levels $$ \Delta(\tau) = Q_1(\tau) - Q_0(\tau), \ \ \tau \in (0,1). $$

We can also obtain the average effect using the relationship between averages and distributions. Thus, the average effect is $$ \Delta = \mu_1 - \mu_0, $$ where $\mu_k$ is the counterfactual average obtained from $F_k$ as

The integral in (ref) is over the real line, but the formula nevertheless is applicable to the case where the support of $dF_k$ is discrete or mixed.

The choice of the levels $t_{ij}^0$ and $t_{ij}^1$ is usually based on the scale of the treatment:

For example, in the trade application we use the levels $t_{ij}^0 = 0$ and $t_{ij}^1 = 1$ for binary covariates such as the indicators for common legal system and free trade area, and $t_{ij}^0 = t_{ij}$ and $t_{ij}^1 = t_{ij} +\log 2$ for the logarithm of distance.

All the previous estimands have causal interpretation under the standard unconfoundedness or conditional independence assumption for panel data where the conditioning set includes not only the observed controls but also the unobserved effects.

Fixed Effects Estimation and Uniform Inference

To simplify the notation in this section we write $P(x_{ij}) = x_{ij}$ without loss of generality, and define $\alpha_i(y) := \alpha(v_i,y)$ and $\gamma_{j}(y) := \gamma(w_{j},y)$.

Fixed Effects Distribution Regression Estimator

The parameters of the DR model can be estimated from multiple binary regressions with two-way effects. To see this, note that the conditional distribution in (ref) can be expressed as $$ \Lambda_y(x_{ij}'\beta(y) + \alpha_i(y) + \gamma_{j}(y)) = {\mathbb{E}}[1\{y_{ij} \leq y \} \mid x_{ij}, v_i, w_{j}] . $$ Accordingly, we can construct a collection of binary variables, $$1\{y_{ij} \leq y\}, \quad (i,j) \in \mathcal{D}, \quad y \in \mathcal{Y},$$ and estimate the parameters for each $y$ by conditional maximum likelihood with fixed effects. Thus, $\widehat \theta(y) := (\widehat \beta(y),$ $\widehat \alpha_1(y), \ldots, \widehat \alpha_I(y),$ $\widehat \gamma_1(y), \ldots, \widehat \gamma_{J}(y))$, the fixed effects distribution regression estimator of $\theta(y) := (\beta(y),$ $\alpha_1(y), \ldots, \alpha_I(y),$ $\gamma_1(y), \ldots, \gamma_{J}(y))$, is obtained as

for $y \in \mathcal{Y}$. When the link function is the normal or logistic distribution, the previous program is concave and smooth in parameters and therefore has good computational properties. See FernandezValWeidner2016, CFW-16 and alpaca17 for a discussion on computation of logit and probit regressions with two-way effects and available software.

The quantile functions and effects are estimated via plug-in rule, i.e., $$ \widehat Q_k(\tau) = \widehat F_k^{\leftarrow}(\tau) \wedge \sup \{ y \in {\mathcal{Y}}\}, \quad \tau \in (0,1), \quad k \in \{0,1\}, $$ where $$ \widehat F_k(y) = n^{-1} \sum_{(i,j) \in \mathcal{D}} \Lambda_y((t_{ij}^k,z_{ij}')' \widehat \beta(y) + \widehat \alpha_i(y) + \widehat \gamma_{j}(y)), \ \ y \in {\mathcal{Y}}, $$ and $$ \widehat \Delta(\tau) = \widehat Q_1(\tau) - \widehat Q_0(\tau) \quad \tau \in (0,1). $$

{

Incidental Parameter Problem and Bias Corrections

Fixed effects estimators can be severely biased in nonlinear models because of the incidental parameter problem (NeymanScott1948). These models include the binary regressions that we estimate to obtain the DR coefficients and estimands. We deal with the incidental parameter problem using the analytical bias corrections of FernandezValWeidner2016 for parameters and average partial effects (APE) in binary regressions with two-way effects. We note here that the distributions $F_0(y)$ and $F_1(y)$ can be seen as APE, i.e., they are averages of functions of the data, unobserved effects and parameters.

The bias corrections are based on expansions of the bias of the fixed effects estimators as $I,J \to \infty$. For example, Theorem (ref) shows that

where $n R^{(F)}_k(y) = o(I \vee J)$.\footnote{FernandezValWeidner2016 considered the case where $n = IJ$, i.e., there is no missing data, so that $I/n = 1/J$ and $J/n = 1/I$.} In Section (ref) we establish that this expansion holds uniformly in $y \in \mathcal{Y}$ and $k \in \{0,1\}$, i.e., $$\sup_{k \in \{0,1\}, y \in \mathcal{Y}} \| n R^{(F)}_k(y) \| = o(I \vee J).$$ This result generalizes the analysis of FernandezValWeidner2016 from a single binary regression to multiple (possibly a continuum) of binary regressions. This generalization is required to implement our inference methods for quantile functions and effects.

The expansion (ref) is the basis for the bias corrections. Let $\widehat B^{(F)}_k(y)$ and $\widehat D^{(F)}_k(y)$ be estimators of $B^{(F)}_k(y)$ and $D^{(F)}_k(y)$, which are uniformly consistent in $y \in \mathcal{Y}$ and $k \in \{0,1\}$. Bias corrected fixed effects estimators of $F_k$ and $Q_k$ are formed as

We also use the corrected estimators $\widetilde F_k$ as the basis for inference and to form a bias corrected estimator of the average effect.

Uniform Inference

One inference goal is to construct confidence bands that cover the QF $\tau \mapsto Q_k(\tau)$ and the QEF $\tau \mapsto \Delta(\tau)$ simultaneously over a set of quantiles $\mathcal{T} \subseteq [\varepsilon,1-\varepsilon]$, for some $0 < \varepsilon < 1/2$, and treatment levels $k \in \mathcal{K} \subseteq \{0,1\}$. The set $\mathcal{T}$ is chosen such that $Q_k(\tau) \in [\inf \{y \in \mathcal{Y}\}, \sup \{y \in \mathcal{Y} \}]$, for all $\tau \in \mathcal{T}$ and $k \in \mathcal{K}$.

We use the generic method of CFVMW2015 to construct confidence bands for quantile functions and effects from confidence bands for the corresponding distributions. Let $\mathbb{D}$ denote the space of weakly increasing functions, mapping ${\mathcal{Y}}$ to $[0,1]$. Assume we have a confidence band $I_k = [L_k,U_k]$ for $F_k$, with lower and upper endpoint functions $y \mapsto L_k(y)$ and $y \mapsto U_k(y)$ such that $L_k, U_k \in \mathbb{D}$ and $L_k(y) \leq U_k(y)$ for all $y \in \mathcal{Y}$.\footnote{If $[L_k',U_k']$ is a confidence band for $F_k$ that does not obey the constraint $L_k',U_k' \in \mathbb{D}$, we can transform $[L_k',U_k']$ into a new band $[L_k,U_k]$ such that $L_k,U_k \in \mathbb{D}$ using the rearrangement method of ChernozhukovFernandezValGalichon2009.} We say that $I_k$ covers $F_k$ if $F_k \in I_k$ pointwise, namely $L_k(y) \leq F_k(y) \leq U_k(y)$ for all $y \in {\mathcal{Y}}$. If $U_k$ and $L_k$ are some data-dependent bands, we say that $I_k$ is a confidence band for $F_k$ of level $p$, if $I_k$ covers $F_k$ with probability at least $p$. Similarly, we say that the set of bands $\{I_k : k \in \mathcal{K}\}$ is a joint confidence band for the set of functions $\{F_k : k \in \mathcal{K}\}$ of level $p$, if $I_k$ covers $F_k$ with probability at least $p$ simultaneously over $k \in \mathcal{K}$. The index set $\mathcal{K}$ can be a singleton to cover individual confidence bands or $\mathcal{K} = \{0,1\}$ to cover joint confidence bands. In Section (ref) we provide a multiplier bootstrap algorithm for computing joint confidence bands based on the joint asymptotic distribution of the bias corrected estimators $\{\widetilde F_k : k \in \mathcal{K}\}$.

The following result provides a method to construct joint confidence bands for $\{Q_k = F_k^\leftarrow : k \in \mathcal{K}\}$, from joint confidence bands for $\{ F_k : k \in \mathcal{K} \}$.

This Lemma establishes that we can construct confidence bands for quantile functions by inverting the endpoint functions of confidence bands for distribution functions. The geometric intuition is that the inversion amounts to rotate and flip the bands, and these operations preserve coverage.

We next construct simultaneous confidence bands for the quantile effect function $\tau \mapsto \Delta(\tau)$ defined by $$\Delta(\tau) = Q_1(\tau) - Q_0(\tau) = F_1^\leftarrow(\tau) - F_0^\leftarrow(\tau), \quad \tau \in \mathcal{T}.$$ The basic idea is to take appropriate differences of the bands for the quantile functions $Q_1$ and $Q_0$ as the confidence band for the quantile effect. Specifically, suppose we have the set of confidence bands $\{I^\leftarrow_k = [U^\leftarrow_k, L^\leftarrow_k] : k = 0,1\}$ for the set of functions $\{F_k^\leftarrow: k = 0,1\}$ of level $p$. CFVMW2015 showed that a confidence band for the difference $Q_1 - Q_0$ of size $p$ can be constructed as $[U^\leftarrow_1 - L^\leftarrow_0, L^\leftarrow_1 - U^\leftarrow_0]$, i.e., $I^\leftarrow_1 \ominus I^\leftarrow_0$ where $\ominus$ is the pointwise Minkowski difference.

Asymptotic Theory

This section derives the asymptotic properties of the fixed effect estimators of $y \mapsto\beta(y)$ and $\{ F_k : k \in {\mathcal K} \}$, as both dimensions $I$ and $J$ grow to infinity. We focus on the case where the link function is the logistic distribution at all levels, $\Lambda_y = \Lambda$, where $\Lambda(\xi) = (1+ \exp(-\xi))^{-1}$. We choose the logistic distribution for analytical convenience. In this case the Hessian of the log-likelihood function does not depend on $y_{it}$, leading to several simplifications in the asymptotic expansions. In particular, there are various terms that drop out from the second order expansions that we use to characterize the structure of the incidental parameter bias of the estimators $\widehat \beta(y)$ and $\widehat F(y)$. For the case of single binary regressions, FernandezValWeidner2016 showed that the properties of fixed effects estimators are similar for the logistic distribution and other smooth log-concave distributions such as the normal distribution. Accordingly, we expect that our results can be extended to other link functions, but at the cost of more complicated proofs and derivations to account for additional terms.

We make the following assumptions:

Asymptotic Distribution of the Uncorrected Estimator

We introduce first some further notation. Denote the $q^{th}$ derivatives of the cdf $\Lambda$ by $\Lambda^{(q)}$, and define $\Lambda^{(q)}_{ij}(y) = \Lambda^{(q)}( x_{ij}'\beta(y) + \alpha_{i}(y) + \gamma_{j}(y) )$ and $\Lambda^{(q)}_{ij,k}(y) = \Lambda^{(q)}( \mathbbm{x}_{ij,k}'\beta(y) + \alpha_{i}(y) + \gamma_{j}(y) )$ with $\mathbbm{x}_{ij,k}:=(t_{ij}^k,z_{ij}')'$ and $q = 1,2, \ldots$. For $\ell \in \{1,\ldots, d_x\}$ define the following projections of the $\ell$'th covariate $x^\ell_{ij}$,

and let $ \alpha_{x,i}(y)$ and $\gamma_{x,j}(y)$ be the $d_x$-vectors with components $\alpha^\ell_{x,i}(y)$ and $\gamma^\ell_{x,j}(y)$, where $\alpha^\ell_{x,i}(y)$ is the $i$th component of $\alpha^\ell_{x}(y)$ and $\gamma^\ell_{x,j}(y)$ is the $j$th component of $\gamma^\ell_{x}(y)$. Also define $\widetilde x_{ij}(y) = x_{ij} - \alpha_{x,i}(y) - \gamma_{x,j}(y)$ and $\widetilde{\mathbbm{x}}_{ij,k}(y) = \mathbbm{x}_{ij,k} - \alpha_{x,i}(y) - \gamma_{x,j}(y)$. Notice that $\widetilde{\mathbbm{x}}_{ij,k}(y)$ is defined using projections of $x_{ij}$ instead of ${\mathbbm{x}}_{ij,k}$. Also, while the locations of $\alpha_{x,i}(y)$ and $\gamma_{x,j}(y)$ are not identified, $\widetilde x_{ij}(y)$ and $\widetilde{\mathbbm{x}}_{ij,k}(y)$ are uniquely defined. Analogous to the projection of $x_{ij}^{\ell}$ above, we define $ \Psi_{ij,k}(y) = \alpha^\Psi_{i}(y) + \gamma^\Psi_{j}(y) $, where

For example, if $\mathbbm{x}_{ij,k} = x_{ij}$, then $\Psi_{ij,k}(y)=1$. Furthermore, we define\footnote{ The FOC of problem (ref) imply that $\sum_{(i,j) \in \mathcal{D}} \Lambda^{(1)}_{ij,k}(y) \, \widetilde{x}_{ij}(y)^{\, \prime} =0$, and we can therefore equivalently write $ \partial_{\beta} F_k(y) = \frac 1 {n} \sum_{(i,j) \in \mathcal{D}} \Lambda^{(1)}_{ij,k}(y) \, \left[ \widetilde{\mathbbm{x}}_{ij,k}(y) - \widetilde{x}_{ij}(y) \right]^{\, \prime} = \frac 1 {n} \sum_{(i,j) \in \mathcal{D}} \Lambda^{(1)}_{ij,k}(y) \, \left[ {\mathbbm{x}}_{ij,k}(y) - {x}_{ij}(y) \right]^{\, \prime}. $ }

and

where $\mathcal{D}_i := \{(i',j') \in \mathcal{D} : i' = i\}$ and $\mathcal{D}_j := \{(i',j') \in \mathcal{D} : j' = j\}$ are the subsets of observational units that contain the index $i$ and $j$, respectively. In the previous expressions, $\partial_{\beta} F_k(y)$ is a $1 \times d_x$ vector for each $k \in {\mathcal K}$ that we stack in the $|{\mathcal K}| \times d_x$ matrix $\partial_{\beta} F(y) = [\partial_{\beta} F_{k}(y) \, : \, k \in {\mathcal K}]$. Similarly, $F_{k}(y)$, $B^{(\Lambda)}_{k}(y)$, $ D^{(\Lambda)}_{k}(y)$ and $\Psi_{ij,k}(y)$ are scalars for each $k \in {\mathcal K}$, that we stack in the $|{\mathcal K}| \times 1$ vectors $ F(y) = [F_{k}(y) \, : \, k \in {\mathcal K}]$, $ B^{(\Lambda)}(y) = [B^{(\Lambda)}_{k}(y) \, : \, k \in {\mathcal K}]$, $ D^{(\Lambda)}(y) = [D^{(\Lambda)}_{k}(y) \, : \, k \in {\mathcal K}]$, $ \Psi_{ij}(y) = [\Psi_{ij,k}(y) \, : \, k \in {\mathcal K}]$.

Let $\ell^{\infty}(\mathcal{Y})$ be the space of real-valued bounded functions on $\mathcal{Y}$ equipped with the sup-norm $\| \cdot \|_{\mathcal{Y}}$, and $\rightsquigarrow$ denote weak convergence (in distribution). We establish a functional central limit theorem for the fixed effects estimators of $y \mapsto \beta(y)$ and $y \mapsto F(y)$ in $\mathcal{Y}$. All stochastic statements are conditional on $ \{(x_{ij}, v_i, w_j): (i,j) \in \mathcal{D}\}$.

Assumption (ref)(vi) guarantees the invertibility of $W(y)$ and $\overline W(y) $. Notice that $\overline W(y) $ is equal to the limit of $W(y)$ because $\Lambda^{(1)}_{ij}(y) = \Lambda_{ij}(y) \left[ 1 - \Lambda_{ij}(y) \right]$ by the properties of the logistic distribution. This information equality follows by the correct specification condition in Assumption (ref)(ii). By Assumption (ref)(v), we could have used $\sqrt{IJ}$ instead of $\sqrt{n}$, $1/J$ instead of $I/n$, and $1/I$ instead of $J/n$. However, if the panel is not balanced, then we expect the expressions in the theorem to provide a more accurate finite-sample approximation, because the standard deviation of the estimates will generally be of order $1/\sqrt{n}$ for unbalanced panels, and the leading order incidental parameter biases are generally proportional to the number of incidental parameters ($I$ and $J$ here) divided by the total sample size $n$, see e.g. FernandezValWeidner2018.

Bias Corrections

Theorem (ref) shows that the fixed effects DR estimator has asymptotic bias of the same order as the asymptotic standard deviation under the approximation that we consider. The finite-sample implications are that this estimator can have substantial bias and that confidence regions constructed around it can have severe undercoverage. We deal with these problems by removing the first order bias of the estimator.

We estimate the bias components using the plug-in rule. Define $\widehat \Lambda^{(q)}_{ij}(y) = \Lambda^{(q)}( x_{ij}' \widehat \beta(y) + \widehat \alpha_{i}(y) + \widehat \gamma_{j}(y) )$ and $\widehat \Lambda^{(q)}_{ij,k}(y) = \widehat \Lambda^{(q)}( \mathbbm{x}_{ij,k}' \widehat \beta(y) + \widehat \alpha_{i}(y) + \widehat \gamma_{j}(y) )$. Replacing $\Lambda^{(1)}_{ij}(y)$ and $\Lambda^{(1)}_{ij,k}(y)$ by $\widehat \Lambda^{(1)}_{ij}(y)$ and $\widehat \Lambda^{(1)}_{ij,k}(y)$ in the definitions of $\alpha^\ell_x(y)$, $\gamma^\ell_x(y)$, $\alpha^\Psi(y)$, and $\gamma^\Psi(y)$ yields the corresponding estimators. We plug-in these estimators to obtain $\widehat x_{ij}(y) = x_{ij} - \widehat \alpha_{x,i}(y) - \widehat \gamma_{x,j}(y)$, $\widehat {\mathbbm x}_{ij,k}(y) = \mathbbm{x}_{ij,k} - \widehat \alpha_{x,i}(y) - \widehat \gamma_{x,j}(y)$, and $\widehat \Psi_{ij,k}(y) = \widehat \alpha^\Psi_{i}(y) + \widehat \gamma^\Psi_{j}(y) $. Then we construct

and

We also define the $|{\mathcal K}| \times d_x$ matrix $\partial_{\beta} \widehat F(y) = [(\partial_{\beta} \widehat F_{k}(y)) \, : \, k \in {\mathcal K}]$, and the $|{\mathcal K}| \times 1$ vectors $ \widehat B^{(F)}(y) = [\widehat B^{(F)}_{k}(y) \, : \, k \in {\mathcal K}]$, $ \widehat D^{(F)}(y) = [\widehat D^{(F)}_{k}(y) \, : \, k \in {\mathcal K}]$, $ \widehat \Psi_{ij}(y) = [\widehat \Psi_{ij,k}(y) \, : \, k \in {\mathcal K}]$. Finally, we also construct the estimator of the asymptotic variance of $\widehat F(y)$

where $\widehat \Xi(y) = \widehat \Psi_{ij}(y) + (\partial_{\beta} \widehat F(y)) \widehat W^{-1}(y) \, \widehat x_{ij}(y)$.

Lemma (ref) in the Appendix shows that the estimators of the asymptotic bias are consistent, uniformly in $y \in \mathcal{Y}$. Bias corrected estimators of $\beta(y)$ and $F(y)$ can then be formed as

and

Alternatively, we could define the bias corrected version of $ \widehat F(y) $ as

where $\widetilde \xi(y) := (\widetilde \alpha_1(y), \ldots, \widetilde \alpha_I(y), \widetilde \gamma_{1}(y), \ldots, \widetilde \gamma_J(y))$ is a solution to $$ \max_{\xi \in \mathbb{R}^{I+J}} \sum_{(i,j) \in \mathcal{D}} (1\{y_{ij} \leq y \} \log \Lambda(x_{ij}'\widetilde \beta(y) + \alpha_i + \gamma_{j}) + 1\{y_{ij} > y \} \log [1 - \Lambda(x_{ij}'\widetilde \beta(y) + \alpha_i + \gamma_{j})]). $$ It can be shown that $\sup_{y \in \mathcal{Y}} \sqrt{n} \left| \widetilde F^*_k(y) - \widetilde F_k(y) \right| = o_P(1)$, that is, the difference between those alternative bias corrected estimators is asymptotically negligible. There is no obvious reason to prefer one over the other, and we present result for $\widetilde F_k$ in the following, which equivalently hold for $\widetilde F^*_k$.\footnote{We use the estimator $\widetilde F^*_k$ in the numerical examples for computational convenience as the bias correction involves estimating less terms. }

The following main result establishes the functional central limit theorem for the bias corrected estimators and uniform consistency of the estimators of the variance function.

Uniform Confidence Bands and Bootstrap

We show how to construct pointwise and uniform confidence bands for $y \mapsto \beta(y)$ and $y \mapsto F(y)$ on $ \mathcal{Y}$ using Theorem (ref). The uniform bands for $F$ can be used as inputs in Lemmas (ref) and (ref) to construct uniform bands for the QFs $\tau \mapsto Q_k(\tau) = F_k^{\leftarrow}(\tau),$ $k \in \mathcal{K}$, and the QEF $\tau \mapsto \Delta(\tau)$ on $\mathcal{T}$.

Let $\mathcal{B} \subseteq \{1, \ldots, d_x\}$ be the set of indexes for the coefficients of interest. For given $y \in \mathcal{Y}$, $\ell \in \mathcal{B}$, $k \in \mathcal{K}$, and $p \in (0,1)$, a pointwise $p$-confidence interval for $\beta_{\ell}(y)$, the $\ell$'th component of $\beta(y)$, is

and a pointwise $p$-confidence intervals for $F_k(y)$ is $$[ \widetilde F_k(y) \pm \Phi^{-1}(1-p/2) \widehat{\sigma}_{F_k}(y)],$$ where $\Phi$ denotes the cdf of the standard normal distribution, $\widehat{\sigma}_{\beta_\ell}(y)$ is the standard error of $\widetilde \beta_{\ell}(y)$ given in (ref), and $\widehat{\sigma}_{F_k}(y)$ is the standard error of $\widetilde F_k(y)$ given in (ref). These intervals have coverage $p$ in large samples by Theorem (ref).

We construct joint uniform bands for the coefficients and distributions using Kolmogorov-Smirnov type critical values, instead of quantiles from the normal distribution. A uniform $p$-confidence band joint for the vector of functions $\{ \beta_{\ell}(y) : \ell \in \mathcal{B}, y \in \mathcal{Y}\}$ is

where $t_{\mathcal{B}, \mathcal{Y}}^{(\beta)}(p)$ is the $p$-quantile of the maximal $t$-statistic

where $\sigma^{(\beta)}_{\ell}(y) = [\overline{W}(y)^{-1}]_{\ell,\ell}^{1/2},$ the square root of the $(\ell,\ell)$ element of the matrix $\overline{W}(y)^{-1}$. Similarly, a uniform $p$-confidence band joint for the set of distribution functions $\{ F_k(y) : k \in \mathcal{K}, y \in \mathcal{Y}\}$ is

where $t_{\mathcal{K}, \mathcal{Y}}^{(F)}(p)$ is the $p$-quantile of the maximal $t$-statistic

where $\sigma^{(F)}_{k}(y) = [\overline{\Omega}(y)]_{k,k}^{1/2}$, the square root of the $(k,k)$ element of the matrix $\overline{\Omega}(y,y)$. The previous confidence bands also have coverage $p$ in large samples by Theorem (ref).

The maximal t-statistics used to construct the bands $I_{\beta}$ and $I_{F}$ are not pivotal, but their distributions can be approximated by simulation after replacing the variance functions of the limit processes by uniformly consistent estimators. In practice, however, we find it more convenient to use resampling methods. We consider a multiplier bootstrap scheme that resamples the efficient scores or influence functions of the fixed effects estimators $\widehat \beta(y)$ and $\widehat F(y)$. This scheme is computationally convenient because it does not need to solve the high dimensional nonlinear fixed effects conditional maximum likelihood program (ref) or making any bias correction in each bootstrap replication. In these constructions we rely on the uncorrected fixed effects estimators instead of the bias corrected estimators, because they have the same influence functions and the uncorrected estimators are consistent under the asymptotic approximation that we consider.

To describe the standard errors and multiplier bootstrap we need to introduce some notation for the influence functions of $\widehat \theta(y)$ and $\widehat F(y)$. Let $\theta = (\beta, \alpha_1, \ldots, \alpha_I, \gamma_1, \ldots, \gamma_J)$ be a generic value for the parameter $\theta(y)$, the influence function of $\widehat \theta(y)$ is the $(d_x+I+J)$-vector $\psi_{ij}^y(\theta(y))$, where $$ \psi_{ij}^y(\theta) = H(\theta)^{\dagger} [\boldsymbol{1}\{y_{ij} \leq y \} - \Lambda(x_{ij}'\beta + \alpha_i + \gamma_j)] w_{ij}, \ \ w_{ij} = (x_{ij}, e_{i,I}, e_{j,J}) , \ \ y \in \mathcal{Y}, $$ $e_{i,I}$ is a unit vector of dimension $I$ with a one in the position $i$, $e_{j,J}$ is defined analogously, $H(\theta)^{\dagger}$ is the Moore-Penrose pseudo-inverse of $H(\theta)$, and $$ H(\theta) = \frac{1}{n} \sum_{(i,j)\in \mathcal{D}} \Lambda^{(1)}(x_{ij}'\beta + \alpha_i + \gamma_j) w_{ij} w_{ij}', \ \ \Lambda^{(1)}(z) = \Lambda(z) \Lambda(-z), $$ is minus the Hessian of the log-likelihood with respect to $\theta$, which does not depend on $y$ in the case of the logistic distribution.\footnote{We use the Moore-Penrose pseudo-inverse because $H(\theta)$ is singular if we do not impose a normalization on the location of $\alpha_i(y)$ and $\gamma_j(y)$.} The influence function of $\widehat F_k(y)$ is $\varphi_{ij,k}^y(\theta(y))$, where $$ \varphi_{ij,k}^y(\theta) = J_k(\theta)' \psi_{ij}^y(\theta), $$ and $$ J_k(\theta) = \frac{1}{n} \sum_{(i,j)\in \mathcal{D}} \Lambda^{(1)}(\mathbbm{x}_{ij,k}'\beta + \alpha_i + \gamma_j) \mathbbm{w}_{ij,k}, \ \ \mathbbm{w}_{ij,k} = (\mathbbm{x}_{ij,k}, e_{i,I}, e_{j,J}). $$

The standard error of $\widetilde \beta_{\ell}(y)$ is constructed as

the square root of the $(\ell,\ell)$ element of the sandwich matrix $n^{-2} \sum_{(i,j) \in \mathcal{D}} \psi_{ij}^y(\widehat \theta(y)) \psi_{ij}^y(\widehat \theta(y))'$. Similarly, the standard error of $\widetilde F_k(y)$ is constructed as

The following algorithm describes a multiplier bootstrap scheme to obtain the critical values for a set of parameters indexed by $\ell \in \mathcal{B} \subseteq \{1, \ldots, d_x\}$ and a set of distributions indexed by $k \in \mathcal{K} \subseteq \{0,1\}$. This scheme is based on perturbing the first order conditions of the fixed effects estimators with random multipliers independent from the data.

The next result shows that the multiplier bootstrap provides consistent estimators of the critical values of the inferential statistics. The proof follows from Theorem 2.2 of cck-16.

Theorem (ref) together with Theorem (ref) guarantee the asymptotic validity of the confidence bands $I_{\beta}$ and $I_F$ defined in (ref) and (ref) with the critical values $t^{(\beta)}_{\mathcal{B}, \mathcal{Y}}(p)$ and $t^{(F)}_{\mathcal{K}, \mathcal{Y}}(p)$ replaced by the bootstrap estimators $\widehat{t}^{(\beta)}_{\mathcal{B}, \mathcal{Y}}(p)$ and $\widehat{t}^{(F)}_{\mathcal{K}, \mathcal{Y}}(p)$.

Pairwise Clustering Dependence or Reciprocity

The conditional independence of Assumption (ref)(i) can be relaxed to allow for some forms of conditional weak dependence. A form of dependence that is relevant for network data is pairwise clustering or reciprocity where the observational units with symmetric indexes $(i,j)$ and $(j,i)$ might be dependent due to unobservable factors not accounted by unobserved effects.\footnote{CM-14 consider other patterns of dependence in linear models for dyadic data.} In the trade application, for example, these factors may include distributional channels or multinational firms operating in both countries. Formally, pairwise clustering means that $(y_{ij},y_{ji})$ is independently distributed across $(i,j) \in \mathcal{D}$ with $i \leq j$, conditional on all the observed and unobserved covariates $ \mathcal{C}_B:= \{(x_{ij}, v_i, w_j): (i,j) \in \mathcal{D}\}$.

The presence of reciprocity does not change the bias of the fixed effects estimators, but affects the standard errors and the implementation of the multiplier bootstrap. The standard error of $\widetilde \beta_{\ell}(y)$ becomes

Similarly, the standard error of $\widetilde F_k(y)$ needs to be adjusted to

In the previous expressions we assume that if $(i,j) \in \mathcal{D}$ then $(j,i) \in \mathcal{D}$ to simplify the notation. The modified multiplier bootstrap algorithm becomes:

The clustered multiplier bootstrap preserves the dependence in the symmetric pairs $(i,j)$ and $(j,i)$ by assigning the same multiplier to each of these pairs.

Average Effect

A bias corrected estimator of the average effect can be formed as

where $$ \widetilde \mu_k = \int [1(y \geq 0) - \mathbf{C} \widetilde F_k(y)] dy, \ \ k \in \{0,1\}. $$ Here the integral is over the real line, and $\mathbf{C}$ is an operator that extends $\widetilde F_k(y)$ from $\mathcal{Y}$ to $\mathbb{R}$ as a step function, that is, it maps any $f: \mathcal{Y} \to \mathbb{R}$ to $\mathbf{C} f: \mathbb{R} \to \mathbb{R}$, where $\mathbf{C} f(y) = 0$ for $y\leq \inf \mathcal{Y}$, $\mathbf{C} f(y) = 1$ for $y\geq \sup \mathcal{Y}$, and $\mathbf{C} f(y) = f(\sup\{ y' \in \mathcal{Y} : y' \leq y \})$ otherwise. The following central limit theorem for the bias corrected estimator of the average effect is a corollary of Theorem (ref) together with the functional delta method.

We can construct confidence intervals for the average effect using Corollary (ref). Let $$ \widehat{\sigma}_{\Delta} = n^{-1} \left[\sum_{(i,j) \in \mathcal{D}} \widehat \varphi_{ij}^2 \right]^{1/2}, \ \ \widehat \varphi_{ij} = - \int \left[\mathbf{C}\varphi_{ij,1}^y(\widehat \theta(y)) - \mathbf{C}\varphi_{ij,0}^y(\widehat \theta(y)) \right] dy. $$ Then, $\widehat{\sigma}_{\Delta}$ is an estimator of $\sigma_{\Delta}$, the standard deviation of the limit process $Z^{(\Delta)}$ in (ref), and $$ I_{\Delta} = [\widetilde{\Delta} \pm \Phi^{-1}(1-p/2) \widehat{\sigma}_{\Delta}], $$ is an asymptotic $p$-confidence interval for $\Delta$. The normal critical value $\Phi^{-1}(1-p/2)$ can be replaced by a multiplier bootstrap critical value $\widehat{t}^{(\Delta)}(p)$ obtained from Algorithm (ref) as $$ \widehat{t}^{(\Delta)}(p) = p-\text{quantile of } \{t^{(\Delta),m} : 1 \leq m \leq M\} $$ where $ t^{(\Delta),m} = |\widehat \Delta^m- \widehat \Delta|/\widehat \sigma_{\Delta} $ and $ \widehat \Delta^m = \widehat \Delta + n^{-1} \sum_{(i,j) \in \mathcal{D}} \omega_{ij}^m \widehat \varphi_{ij}. $

The standard errors and critical values of the average effects can be adjusted to account for pairwise clustering following the procedure described in Section (ref). Thus, the pairwise clustering robust standard error is $$ \widehat{\sigma}_{\Delta} = n^{-1} \left[\sum_{(i,j) \in \mathcal{D}} \left\{ \widehat \varphi_{ij} + \widehat \varphi_{ji}\right\} \widehat \varphi_{ij} \right]^{1/2}. $$

Quantile Effects in Gravity Equations for International Trade

We consider an empirical application to gravity equations for bilateral trade between countries. We use data from Helpman01052008, extracted from the Feenstra's World Trade Flows, CIA's World Factbook and Andrew Rose's web site. These data contain information on bilateral trade flows and other trade-related variables for 157 countries in 1986.\footnote{The original data set includes 158 countries. We exclude Congo because it did not export to any other country in 1986.} The data set contains network data where both $i$ and $j$ index countries as senders (exporters) and receivers (importers), and therefore $I = J = 157$. The outcome $y_{ij}$ is the volume of trade in thousands of constant 2000 US dollars from country $i$ to country $j$, and the covariates $P(x_{ij}) = x_{ij}$ include determinants of bilateral trade flows such as the logarithm of the distance in kilometers between country $i$'s capital and country $j$'s capital and indicators for common colonial ties, currency union, regional free trade area (FTA), border, legal system, language, and religion. Following AndersonWincoop2003, we include unobserved importer and exporter country effects.\footnote{See Harrigan1994 for an earlier empirical international trade application that includes unobserved country effects.} These effects control for other country specific characteristics that may affect trade such as GDP, tariffs, population, institutions, infrastructures or natural resources. We allow for these characteristics to affect differently the imports and exports of each country, and be arbitrarily related with the observed covariates.

Table (ref) reports descriptive statistics of the variables used in the analysis. There are $157 \times 156 = 24,492$ observations corresponding to different pairs of countries. The observations with $i = j$ are missing because we do not observe trade flows from a country to itself. The trade variable in the first row is an indicator for positive volume of trade. There are no trade flows for 55% of the country pairs. The volume of trade variable exhibits much larger standard deviation than the mean. Since this variable is bounded below at zero, this indicates the presence of a very heavy upper tail in the distribution. This feature also makes quantile methods specially well-suited for this application on robustness grounds.\footnote{In results not reported, we find that estimates of average effects are very sensitive to the trimming of outliers at the top of the distribution.}

The previous literature estimated nonlinear parametric models such as Poisson, Negative Binomial, Tobit and Heckman-selection models to deal with the large number of zeros in the volume of trade (e.g., EatonKortum2001, SantosSilvaTenreyro2006, and Helpman01052008).\footnote{See HeadMayer2014 for a recent survey on gravity equations in international trade.} These models impose strong conditions on the process that generates the zeros and/or on the conditional heteroskedasticity of the volume of trade. The DR model deals with zeros and any other fixed censoring points in a very flexible and natural fashion as it specifies the conditional distribution separately at the mass point. In particular, the model coefficients at zero can be arbitrarily different from the model coefficients at other values of the volume of trade. Moreover, the DR model can also accommodate conditional heteroskedasticity.

Figure (ref) shows estimates and 95% pointwise confidence intervals for the DR coefficients of log distance and common legal system plotted against the quantile indexes of the volume of trade. We report uncorrected and bias corrected fixed effects estimates obtained from (ref) and (ref), respectively. The confidence intervals are constructed using (ref). The x-axis starts at .54, the maximum quantile index corresponding to zero volume of trade. The region of interest $\mathcal{Y}$ corresponds to the interval between zero and the $0.95$-quantile of the volume of trade. The difference between the uncorrected and bias corrected estimates is the same order of magnitude as the width of the confidence intervals for the coefficient of log distance. We find the largest estimated biases for both coefficients at highest quantiles of the volume of trade, where the indicators $1\{y_{ij} \leq y\}$ take on many ones. The signs of the DR coefficients indicate that increasing distance has a negative effect and having a common legal system has a positive effect on the volume of trade throughout the distribution. Recall that the sign of the effect in terms of volume of trade, $y_{ij},$ is the opposite to the sign of the DR coefficient.

Figures (ref) and (ref) show estimates and 95% uniform confidence bands for distribution and quantile functions of the volume of trade at different values of the log of distance and the common legal system. The left panels plot the functions when distance takes the observed levels ($\text{dist}$) and two times the observed values $(\text{2*dist})$, i.e. when we counterfactually double all the distances between the countries. The right panels plot the functions when all the countries have the same legal system (legal=1) and different systems (legal=0). The confidence bands for the distribution are obtained by Algorithm (ref) with 500 bootstrap replications and standard normal multipliers, and a grid of values $\bar \mathcal{Y}$ that includes the sample quantiles of the volume of trade with indexes $\{.54, .55, \ldots, .95\}$. The bands are joint for the two functions displayed in each panel. The confidence bands for the quantile functions are obtained by inverting and rotating the bands for the corresponding distribution functions using Lemma (ref).

Figure (ref) displays estimates and 95% uniform confidence bands for the quantile effects of the log of distance and the common legal system on the volume of trade, constructed using Lemma (ref). For comparison, we also include estimates from a Poisson model. Here, we replace the DR estimators of the distributions by

where $\lfloor y \rfloor$ is the integer part of $y$, $\lambda_{ij,k} = \exp(\mathbbm{x}_{ij,k}'\widehat \beta + \widehat \alpha_i + \widehat \gamma_j)$, and $\widehat \theta = (\widehat \beta, \widehat \alpha_1, \ldots, \widehat \alpha_I, \widehat \gamma_1, \dots, \widehat \gamma_J)$ is the Poisson fixed effects conditional maximum likelihood estimator $$ \widehat \theta \in \arg \max_{\theta \in \mathbb{R}^{d_x + I + J}} \sum_{(ij) \in \mathcal{D}} [y_{ij} (x_{ij}'\beta + \alpha_i + \gamma_j) - \exp (x_{ij}'\beta + \alpha_i + \gamma_j)]. $$ We find that distance and common legal system have heterogeneously increasing effects along the distribution. For example, the negative effects of doubling the distance grows more than proportionally as we move up to the upper tail of the distribution of volume of trade. Putting all the countries under the same legal system has little effects in the extensive margin of trade, but has a strong positive effect at the upper tail of the distribution. The Poisson estimates lie outside the DR confidence bands reflecting heavy tails in the conditional distribution of the volume of trade that is missed by the Poisson model.\footnote{This misspecification problem with the Poisson model is well-known in the international trade literature. The Poisson estimator is treated as a quasi-likelihood estimator and standard errors robust to misspecification are reported (SantosSilvaTenreyro2006).} Figure (ref) shows confidence bands of the quantile effects that account for pairwise clustering. The bands are constructed from confidence bands from the distributions using Algorithm (ref) with $500$ bootstrap draws and standard normal multipliers. Accounting for unobservables that affect symmetrically to the country pairs has very little effect on the width of the bands in this case.

Montecarlo Simulation

We conduct a Montecarlo simulation calibrated to the empirical application of Section (ref). The outcome is generated by the censored logistic process $$ y^s_{ij} = \max\{x_{ij}'\widehat \beta + \widehat \alpha_i + \widehat \gamma_{j} + \widehat{\sigma} \Lambda^{-1}(u^s_{ij})/\sigma_{L}, 0 \}, \ \ (i,j) \in \mathcal{D}, $$ where $\mathcal{D} = \{(i,j): 1 \leq i,j \leq 157, i \neq j\}$, $x_{ij}$ is the value of the covariates for the observational unit $(i,j)$ in the trade data set, $\sigma_L = \pi/\sqrt{3}$, the standard deviation of the logistic distribution, and $(\widehat \beta, \widehat \alpha_1, \ldots, \widehat \alpha_I, \widehat \gamma_{1}, \ldots, \widehat \gamma_{J}, \widehat \sigma)$ are Tobit fixed effect estimates of the parameters in the trade data set with lower censoring point at zero.\footnote{We upper winsorize the volume of trade $y_{ij}$ at the $95.5\%$ quantile to reduce the effect of outliers in the Tobit estimation of the parameters.} We consider two designs: independent errors with $u^s_{ij} \sim \text{ i.i.d } \mathcal{U}(0,1),$ and pairwise dependent errors with $u^s_{ij} = \Phi(0.75 e^s_{ij} + \sqrt{1-0.75^2} e^s_{ji}),$ where $e^s_{ij} \sim \text{ i.i.d } \mathcal{N}(0,1)$ and $\Phi$ is the standard normal CDF.\footnote{The Spearman rank correlation between $u^s_{ij}$ and $u^s_{ji}$ in the design with pairwise-dependent errors is $0.73$.} In both cases the conditional distribution function of $y_{ij}^s$ is a special case of the DR model (ref) with link function $\Lambda_y = \Lambda$, the logistic distribution, for all $y$, $$ \beta(y) = \sigma_L (e_1 y - \widehat \beta)/\widehat \sigma, \ \ \alpha_i(y) = -\sigma_L \widehat \alpha_i/\widehat \sigma, \ \text{and} \ \gamma_j(y) = - \sigma_L \widehat \gamma_j/\widehat \sigma, $$ where $e_1$ is a unit vector of dimension $d_x$ with a one in the first component. As in the empirical application, the region of interest $\mathcal{Y}$ is the interval between zero and the $0.95$-quantile of the volume of trade in the data set. All the results are based on 500 simulated panels $\{(y^s_{ij}, x_{ij}) : (i,j) \in \mathcal{D}\}$.

Figures (ref) and (ref) report the biases, standard deviations and root mean square errors (rmses) of the fixed effects estimators of the DR coefficients of log-distance and legal system as a function of the quantiles of $y_{ij}$ in the design with independent errors.\footnote{The design with pairwise dependent errors produces similar results, which are not reported for the sake of brevity.} All the results are in percentage of the true value of the parameter. As predicted by the large sample theory, the fixed effects estimator displays a bias of the same order of magnitude as the standard deviation. As in fig. (ref), the bias is more severe for the coefficient of log distance. The bias correction removes most of the bias and does not increase the standard deviation, yielding a reduction in rmse of about 5% for the coefficient of log distance at the highest quantile indexes.

Figure (ref) reports the biases, standard deviations and rmses of the estimators of the counterfactual distributions at two levels of log-distance as a function of the quantiles of $y_{ij}$ in the design with independent errors. The levels of distance in these distributions are the same as in the empirical application, i.e. $k=0$ and $k=1$ correspond to the observed values and two times the observed values, respectively. All the results are in percentage of the true value of the functions. In this case we find that the uncorrected and bias corrected estimators display small biases relative to their standard deviations, and have similar standard deviations and rmses at both treatment levels. Indeed the standard deviations and rmses are difficult to distinguish in the figure as they are almost superposed. In results not reported, we find very similar patterns in the design with pairwise dependent errors and for the estimators of the counterfactual distributions at the same two levels of legal as in the empirical application.

Table (ref) shows results on the finite sample properties of 95% confidence bands for the DR coefficients and counterfactual distributions in the design with independent errors. The confidence bands are constructed by multiplier bootstrap with 500 draws, standard normal weights, and a grid of values $\bar \mathcal{Y}$ that includes the sample quantiles of the volume of trade with indexes $\{.54, .55, \ldots, .95\}$ in the trade data set. For the coefficients, it reports the average length of the confidence bands integrated over threshold values, the average value of the estimated critical values, and the empirical coverages of the confidence bands. For the distributions, it reports the same measures averaged also over the two treatment levels and where the coverage of the bands is joint for the two counterfactual distributions.\footnote{The joint coverage of the bands for the quantile functions and quantile effect is determined by the joint coverage of the bands of the distribution functions in our construction. We refer to CFVMW2015 for a numerical analysis on the marginal coverage of the bands for the quantile effects.} For comparison, it also reports the coverage of pointwise confidence bands using the normal distribution, i.e. with critical value equal to 1.96. The last row computes the ratio of the standard error averaged across simulations to the simulation standard deviation, integrated over threshold values for the coefficients and over thresholds and treatment levels for the distributions. We consider standard errors and confidence bands with and without accounting for pairwise clustering. All the results are computed for confidence bands centered at the uncorrected fixed effects estimates and at the bias corrected estimates. For the coefficients, we find that the bands centered at the uncorrected estimates undercover the true coefficients, whereas the bands centered at the bias corrected estimates have coverages close to the nominal level. The joint coverage of the bands for the distributions is close to the nominal level regardless of whether they are centered at the uncorrected or bias corrected estimates. We attribute this similarity in coverage to the small biases in the uncorrected estimates of the distributions found in fig. (ref). As expected, pointwise bands severely undercover the entire functions. The standard errors based on the asymptotic distribution provide a good approximation to the sampling variability of both the uncorrected and bias corrected estimators. Accounting for pairwise clustering in this design where it is not necessary has very little effect on the quality of the inference.

Table (ref) reports the same results as table (ref) for the design with pairwise dependent errors. The bands that do not account for pairwise clustering undercover the functions because the standard errors underestimate the standard deviations of the estimators. Compared to the design with independent errors, the critical values are similar but the bands that account for clustering are wider due to the increase in the standard errors. To sum up, inference methods robust to pairwise clustering perform well in both designs, whereas inference methods that do not account for clustering undercover in the presence of pairwise dependence. The bias corrections are effective in reducing bias and bringing the coverage probabilities of the bands close to their nominal level for the coefficients, whereas they have little effect for the distributions.

Conclusion

We have constructed confidence bands for quantile functions and quantile effects in nonlinear network and panel models with two-way unobserved effects. Our construction relies on the generic method of CFVMW2015 to convert confidence bands for distributions into confidence bands for quantiles. The same method can be applied to more complicated models such as nonlinear models with interactive unobserved effects or factor structure, provided that confidence bands for distributions in these models are supplied. Such bands are not currently available, but could be obtained by extending the central limit theorem of Chen, Fern\'andez-Val and Weidner (in press) to a functional central limit theorem. We leave such extension to future work.

\phantom{cfw14}

\setcounter{page}{1} \pagenumbering{roman} \setcounter{section}{0} \setcounter{theorem}{0} \setcounter{corollary}{0} \setcounter{lemma}{0} \setcounter{equation}{0}

The following proof of Lemma (ref) also relies on the results of Lemma (ref) and Lemma (ref)$(i)$, whose proof is presented afterwards, without using Lemma (ref) of course.

\bf Intermediate results for the proof of of Lemma (ref)

The proof of Lemma (ref) requires several intermediate results from the theory of stochastic processes, which are presented in the following. The notation $a_n \lesssim b_n$ means that $a_n \leq C \, b_n$ for some constant $C$ that is independent of the sample size $n$. It is also convenient to define ${\bf I} := \{ 1,2,\ldots, I\}$ and ${\bf J} := \{ 1,2,\ldots, J\}$. In this section we assume that $\mathcal{Y}$ is a bounded interval, and that $y_{ij}$ is continuously distributed with density bounded away from zero. The results for the case where $y_{ij}$ is discrete follow directly from FernandezValWeidner2016. Results for a mixed distribution of $y_{ij}$ follow by combining the results for the continuous and discrete cases.

Bounds on sample averages over the score $\partial_{\pi} \ell_{y,ij}$

For every $i \in {\bf I}$ we define the empirical process

Here, for ease of notation, we use the subscript $J$ to denote the sample size, corresponding to the balanced panel case where $|\mathcal{D}_i|=J$. Following standard notation we write $\| \mathbb{G}_{J,i} \|_{\mathcal F} := \sup_{f \in {\mathcal F} } \left| \mathbb{G}_{J,i} f \right| $. Every element of ${\mathcal F}$ correspond to exactly one $y \in {\mathcal Y} $, and in the following we write $ f_y : \widetilde y \; \mapsto \; 1(\widetilde y \leq y) $ for that element. Since ${\mathbb{E}} \; 1\{y_{ij} \leq y\} = \Lambda_{y,ij}$,

where

Our goal is to show that $\max_{i \in {\bf I}} \| \mathbb{G}_{J,i} \|_{\mathcal F} = o_P\left( I^{1/6} \right)$ by using the following theorem. The theorem uses standard notation $\mathbb{G}_n$ for the empirical process, denoting the sample size by $n$ (not $J$ or $|\mathcal{D}_i|$), and without an extra index $i$. The definition of $J(\delta, {\mathcal F})$ is given on p.239 of van der Vaart and Wellner (1996). All that matters to us is that $J(1, {\mathcal F})$ only depends on ${\mathcal F}$ (not on the probability measure or on the empirical process) and that for the ${\mathcal F}$ defined in (ref) we have $J(1, {\mathcal F}) < \infty$, because ${\mathcal F}$ is VC class. An obvious envelope function for that ${\mathcal F}$ is $F: \tilde y \mapsto 1$, which satisfies $\| F \|_n=1$. The minimal measurable majorant of $\{ \mathbb{G}_{n} f : f \in {\mathcal F} \} $ is denoted by $\| \mathbb{G}_{n} \|^*_{\mathcal F}$, and is identical to $\| \mathbb{G}_{n} \|_{\mathcal F}$ for our purposes.

Bounds on weighted sample averages over the score $\partial_{\pi} \ell_{y,ij}$

We also need results on sample averages of the form e.g. $ \frac 1 {\sqrt{n}} \sum_{(i,j) \in \mathcal{D}} \, b^{(n)}_{y,hij} \, \partial_{\pi} \ell_{y,ij}$, where $b^{(n)}_{y,hij}$ are weights that also depend on the index $y$. The following lemma is useful for that purpose.

Using Lemma (ref) we obtain the following corollary.

FCLT for weighted sample averages over the score $\partial_{\pi} \ell_{y,ij}$

The following theorem will be used in the proof of part $(ii)$ and $(iii)$ of Lemma (ref).

A semi-metric $\rho$ is Gaussian-dominated if it is bounded above by a Gaussian semi-metric. Any semi-metric such that $\int_0^\infty \sqrt{\log N(\epsilon, \mathcal{F}, \rho)} d \epsilon < \infty$ is Gaussian dominated.

\bf Proof of Lemma (ref)

\bf Proof of Lemma (ref)