The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
62,016 characters
Panel Data Quantile Regression for Treatment Effect Models
\doublespacing
\author{Takuya Ishihara\footnote{Tohoku University, Graduate School of Economics and Management, email: [email removed]. I would like to thank the editor, associate editor, and anonymous referees for their careful reading and comments. I would also like to thank Hidehiko Ichimura, Hiroyuki Kasahara, Masayuki Sawada, Katsumi Shimotsu, and the seminar participants at the University of Tokyo, Kobe University, and Tohoku University. This research was supported by a Grant-in-Aid for JSPS Fellows (20J00900) from the JSPS.}}
\title{Panel Data Quantile Regression for Treatment Effect Models}
\date{\today}
\maketitle
\begin{abstract}
In this study, we develop a novel estimation method for quantile treatment effects (QTE) under rank invariance and rank stationarity assumptions.
\cite{ishihara2020identification} explores identification of the nonseparable panel data model under these assumptions and proposes a parametric estimation based on the minimum distance method.
However, when the dimensionality of the covariates is large, the minimum distance estimation using this process is computationally demanding.
To overcome this problem, we propose a two-step estimation method based on the quantile regression and minimum distance methods.
We then show the uniform asymptotic properties of our estimator and the validity of the nonparametric bootstrap.
The Monte Carlo studies indicate that our estimator performs well in finite samples.
Finally, we present two empirical illustrations, to estimate the distributional effects of insurance provision on household production and TV watching on child cognitive development.
\if0
Keywords: continuous treatments, nonseparable models, heterogenous treatment effects
\fi
\end{abstract}
\section{Introduction}
In the literature on program evaluation, it is important to learn about the distributional effects beyond the average effects of the treatment.
Policymakers are more likely to prefer a policy that tends to increase outcomes in the lower tail of the outcome distribution to one that tends to increase outcomes in the middle or upper tail of the outcome distribution.
Such effects can be captured by comparing the quantiles of the treated and control potential outcomes.
The parameter of interest here is the quantile treatment effects (QTE) or the quantile treatment effects on the treated (QTT).
For example, \cite{abadie2002instrumental} estimated the distributional impact of the Job Training Partnership Act (JTPA) program on earnings.
They showed that, for women, the JTPA program had the largest proportional impact at low quantiles.
However, the training impact for men was largest in the upper half of the distribution, with no significant effect on the lower quantiles.
This result could not have been achieved using a mean impact analysis.
Empirical researchers have estimated distributional effects, such as the QTE or QTT, in many areas of empirical economic researches.
For example, \cite{chernozhukov2004effects} estimated the QTE of participation in a 401(k) plan on several measures of wealth; \cite{james2006mean} estimated the QTE of welfare reforms on earnings, transfers, and income; \cite{martincus2010beyond} estimated the QTE of trade promotion activities; and \cite{havnes2015universal} and \cite{kottelenberg2017targeted} estimated the QTT of universal child care.
There is also a rich literature on the identification and estimation of the QTE and QTT in various contexts.
\cite{firpo2007efficient} explored the identification and estimation of the QTE under unconfoundedness.
\cite{abadie2002bootstrap}, \cite{chernozhukov2005iv}, \cite{chernozhukov2006instrumental}, and \cite{frolich2013unconditional} showed how instrumental variables can be used to identify the QTE.
\cite{athey2006identification}, \cite{melly2015changes}, and \cite{callaway2019quantile} provided the identification and estimation results for the QTT in a difference-in-differences (DID) setting by using repeated cross-sections or panel data.
Further, \cite{d2013nonlinear} studied the identification of nonseparable models with continuous treatments using repeated cross sections.
In this study, we use use panel data to develop a novel estimation method for the QTE under rank invariance and rank stationarity assumptions.
We propose a two-step estimator, based on the quantile regression and minimum distance methods.
The rank invariance assumption is used in many nonseparable models, such as those in \cite{matzkin2003nonparametric}, \cite{chernozhukov2005iv}, \cite{d2015identification}, \cite{torgovitsky2015identification}, \cite{feng2020estimation}, and \cite{ishihara2020identification}.
This assumption implies that a scalar unobserved factor determines the potential outcomes across treatment status.
The rank stationarity assumption implies that the conditional distribution of the unobserved factor, given explanatory variables and covariates, does not change over time.
In the literature on nonseparable panel data models, similar assumptions were employed by \cite{athey2006identification}, \cite{hoderlein2012nonparametric}, \cite{graham2012identification}, \cite{d2013nonlinear}, \cite{chernozhukov2013average}, \cite{chernozhukov2015nonparametric}, and \cite{ishihara2020identification}.
\cite{ishihara2020identification} also explores the identification of the nonseparable panel data model under the rank invariance and rank stationarity assumptions.
In this work, the structural function depends on the time period in an arbitrary way and does not require the existence of "stayers" - individuals with the same regressor values in two time periods.
It is important to consider nonlinear time trends when modeling the quantile function.
In this case, additive time trends may be restrictive.
For example, if the quantile function of $Y_t$ is written as $q_t(\tau) = g(\tau) + \mu_t$, the distribution of $Y_t$ is the same across time, up to the location.
However, such an assumption is not valid for many empirical applications.
In contrast, the nonseparable panel data model proposed by \cite{ishihara2020identification} captures nonlinear time effects.
Many nonseparable panel data models require the existence of stayers; this is included in \cite{evdokimov2010identification}, \cite{hoderlein2012nonparametric}, and \cite{chernozhukov2015nonparametric}.
In particular, \cite{evdokimov2010identification} requires the existence of stayers for any value of the treatment variable.
However, many empirically important models do not satisfy this assumption.
For example, in standard DID models, no individuals are treated during both time periods.
The identification approach of \cite{ishihara2020identification} does not require the existence of stayers and allows the support conditions that are employed in standard DID models.
\cite{ishihara2020identification} also proposes a parametric estimation based on the minimum distance method.
However, when the dimensionality of the covariates is large, the minimum distance estimator is computationally demanding.
Hence, when we add many covariates into the model, it is difficult to compute the estimator.
To overcome this problem, we propose a two-step estimation method based on the quantile regression and minimum distance methods.
Using quantile regression, we can obtain an estimator of the QTE by optimizing the objective function over a low-dimensional parameter.
This two-step estimation method is similar to the instrumental variable quantile regression method proposed by \cite{chernozhukov2006instrumental}.
In a DID setting, our model is similar to the changes-in-changes (CIC) model.
\cite{athey2006identification} suggest the CIC model as an alternative to the DID model.
The CIC model allows QTT estimation.
Their model is less restrictive than our model because their approach does not require the rank invariance assumption.
However, their approach does not work when the treatment is a continuous variable or a discrete variable with many different possible values.
There exist many empirical applications in which the treatment variable is continuous, such as when many researchers use panel data to estimate the effect of class size on children's test scores.
Our estimation method, contrary to the CIC model, works when the treatment variable is continuous.
\cite{d2013nonlinear} study the identification of nonseparable models with continuous treatments using repeated cross sections.
They allow for nonlinear time effects by assuming that the structural function $g_t(x,u)$ can be written as $m_t(h(x,u))$, where $m_t$ is a monotonic transformation.
Their study proposes a nonparametric estimation method of the QTT.
However, if we add many covariates into the model, their estimation method does not work because of the curse of dimensionality.
\cite{melly2015changes}, \cite{kottelenberg2017targeted}, and \cite{sawada2019noncompliance} also consider the estimation of the CIC model in the presence of covariates.
\cite{melly2015changes} suggest a flexible semiparametric estimator based on a quantile regression analysis.
They estimate the conditional distribution of outcomes for both treatment and control groups and both periods by using quantile regression, and then apply the changes-in-changes transformations.
Similar to \cite{melly2015changes}, \cite{sawada2019noncompliance} proposes a semiparametric estimator based on distribution regressions.
\cite{kottelenberg2017targeted} rely on the Firpo's (2007) extension to quantiles of the inverse propensity scores method.
However, none of them allow for continuous treatments.
An alternative approach estimates the distributional effects using panel data.
\cite{callaway2019quantile} provide identification and estimation results for the QTT under a straightforward extension of the most common DID assumption.
To identify the QTT, they employ two key assumptions: the distributional difference-in-differences assumption and the copula stability assumption.
The first assumption means that the distribution of the change in potential untreated outcomes does not depend on whether the individual belongs to the treatment or control group.
The second assumption means that the copula between the change in the untreated potential outcomes for the treated group and the initial untreated outcome for the treated group is stable over time.
The rest of the paper is organized as follows.
Section 2 introduces the model and assumptions and demonstrates that our model is nonparametrically identified.
In Section 3, we review the minimum distance estimator, as suggested by \cite{ishihara2020identification}.
We then propose a two-step estimator and show the uniform asymptotic properties of our estimator and the validity of the nonparametric bootstrap.
Section 4 contains the results of several Monte Carlo simulations and we illustrate our estimation method in two empirical settings in Section 5.
The paper concludes in Section 6.
The proofs of the theorems and auxiliary lemmas are provided in the Appendix.
\section{Assumptions and nonparametric identification}
We consider the following potential outcome framework.
The potential outcomes are indexed against the potential values $x$ of the treatment variable $X_{it} \in \mathbb{R}^{d_X}$ and denoted by $Y_{it}(x)$.
We cannot observe $Y_{it}(x)$ directly and the observed outcome is given by $Y_{it}\equiv Y_{it}(X_t)$.
Furthermore, we observe a vector of covariates, $Z_{it}$.
We define $\mathbf{Y}_i \equiv (Y_{i1}, \cdots , Y_{iT})'$, $\mathbf{X}_i \equiv (X_{i1}', \cdots , X_{iT}')'$, $\mathbf{Z}_i \equiv (Z_{i1}', \cdots, Z_{iT}')'$, $W_{it} \equiv (Y_{it},X_{it}',Z_{it}')'$ and $\mathbf{W}_i \equiv (W_{i1}',\cdots,W_{iT}')'$.
Let $\mathcal{X}_t$, $\mathcal{X}_{1:T}$, and $\mathcal{Z}$ denote the supports of $X_{it}$, $\mathbf{X}_i$, and $\mathbf{Z}_i$.
We assume that the potential outcome can be expressed as
\begin{eqnarray}
Y_{it}(x) &=& q_t\left( x,Z_{it},U_{it}(x) \right), \ \ \ \ i = 1, \cdots , n, \ \ t = 1, \cdots, T, \label{potential_outcome}
\end{eqnarray}
where $q_t(x,z_t,\tau)$ is the conditional $\tau$-th quantile of $Y_{it}(x)$ conditional on $\mathbf{Z}_i=(z_1, \cdots, z_T)'$ and $U_{it}(x)$ is uniformly distributed conditional on $\mathbf{Z}_i$.
This implies that the conditional distribution of $Y_{it}(x)$ conditional on $\mathbf{Z}_i$ depends only on $Z_{it}$.
When all covariates are time-invariant, this condition does not restrict the conditional distribution and expression (\ref{potential_outcome}) is known as the Skorohod representation.
Following \cite{chernozhukov2005iv}, we refer to $U_{it}(x)$ as the rank variable.
Additionally, we allow $q_t$ to depend on the time period in an arbitrary manner, similar to the work of \cite{ishihara2020identification}.
First, we impose the rank invariance assumption.
\begin{Assumption}
(i) For all $x$, we have $U_{it}(x) = U_{it}$.
(ii) For all $t$, $U_{it}$ is uniformly distributed on $[0,1]$ conditional on $\mathbf{Z}_i$.
(iii) For all $t$, $\mathbf{x}$, and $\mathbf{z}$, the support of $U_{it}|\mathbf{X}_i=\mathbf{x},\mathbf{Z}_i=\mathbf{z}$ is $[0,1]$.
\end{Assumption}
Assumption 1 (i) is referred to as the rank invariance assumption.
For example, \cite{matzkin2003nonparametric}, \cite{chernozhukov2005iv}, \cite{d2015identification}, \cite{torgovitsky2015identification}, \cite{feng2020estimation}, and \cite{ishihara2020identification} also employ similar assumptions.
This model is restrictive because the potential outcomes $\{Y_{it}(x)\}_{x \in \mathcal{X}_t}$ are not truly multivariate and, are jointly degenerate.
As discussed in \cite{chernozhukov2005iv}, we can relax the rank invariance assumption to the rank similarity assumption. That is, $U_{it}(x)|\mathbf{X}_i,\mathbf{Z}_i \overset{d}{=} U_{it}(\tilde{x})|\mathbf{X}_i,\mathbf{Z}_i$ for all $x$ and $\tilde{x}$.
Under the rank invariance assumption, the observed outcome can be written as
\begin{eqnarray}
Y_{it} &=& q_t \left( X_{it}, Z_{it}, U_{it} \right), \ \ \ \ i = 1, \cdots , n, \ \ t = 1, \cdots, T. \label{observed_outcome}
\end{eqnarray}
This is the nonseparable model with a scalar unobserved variable and model (\ref{observed_outcome}) is the same as the model proposed by \cite{ishihara2020identification} when there are no covariates.
If $U_{it}$ is independent of $X_{it}$ and $Z_{it}$, then this model is identical with the usual quantile regression model.
However, our model allows for correlation between $U_{it}$ and the treatment variable.
Hence, to achieve point identification, we require additional assumptions.
Next, we impose the rank stationarity assumption.
\begin{Assumption}
For all $t \neq s$, $\mathbf{x}$, and $\mathbf{z}$, the conditional distribution of $U_{it}|\mathbf{X}_i=\mathbf{x},\mathbf{Z}_i =\mathbf{z}$ is the same as that of $U_{is}|\mathbf{X}_i=\mathbf{x}, \mathbf{Z}_i = \mathbf{z}$.
\end{Assumption}
Assumption 2 implies that the rank variable is stationary across the time period.
In the literature on nonseparable panel data models, similar assumptions were employed by \cite{athey2006identification}, \cite{hoderlein2012nonparametric}, \cite{graham2012identification}, \cite{d2013nonlinear}, \cite{chernozhukov2013average}, \cite{chernozhukov2015nonparametric}, and \cite{ishihara2020identification}.
\cite{chernozhukov2013average} referred to Assumption 2 as ``time is randomly assigned'' or ``time is an instrument.''
Assumption 2 can be viewed as a quantile version of the identification condition of the following conventional linear panel data model:
$$
Y_{it} = X_{it}'\alpha + A_i + \epsilon_{it}, \ \ \ \ E[X_{is} \epsilon_{it}]=0 \ \text{for all $t$ and $s$,}
$$
where $A_i$ is a fixed effect and $\epsilon_{it}$ is a time-variant unobserved variable.
Let $\bar{E}[\cdot|\mathbf{X}_i]$ denote the linear projection on $\mathbf{X}_i$, as in \cite{chamberlain1982multivariate}.
\cite{chernozhukov2013average} show that the above equation is satisfied if and only if there is $\tilde{\epsilon}_{it}$ with
$$
Y_{it} = X_{it}'\alpha + \tilde{\epsilon}_{it}, \ \bar{E}[\tilde{\epsilon}_{it}|\mathbf{X}_i]=\bar{E}[\tilde{\epsilon}_{is}|\mathbf{X}_i] \ \text{for all $t$ and $s$.}
$$
In contrast, if the conditional quantile function is linear in $X_{it}$ and there are no covariates, then we can rewrite model (\ref{observed_outcome}) as
\begin{equation}
Y_{it} = X_{it}'\alpha(\tau) + \epsilon_{it}(\tau), \nonumber
\end{equation}
where $\epsilon_{it}(\tau) = X_{it}'(\alpha(U_{it})-\alpha(\tau))$.
Then, under Assumption 2, $\epsilon_{it}(\tau)$ satisfies $F_{\epsilon_t(\tau)| \mathbf{X}}(0|\mathbf{x}) = F_{\epsilon_s(\tau)| \mathbf{X}}(0|\mathbf{x})$ for all $t\neq s$ and $\mathbf{x}$.
Hence, the rank stationarity assumption can be viewed as a quantile version of the identification condition of the conventional linear panel data model.
Under Assumptions 1 and 2 and additional assumptions in Appendix 1, we can show that $q_t(x,z_t,\tau)$ is nonparametrically identified.
The following proposition is essentially the same as Corollary 1 in \cite{ishihara2020identification}.
\begin{Proposition}
Under Assumptions 1, 2, A.1, and A.2, the conditional quantile function $q_t$ is point identified.
\end{Proposition}
From the proof of Proposition 1, for any $t \neq s$, we have
$$
F_{Y_t|\mathbf{X},\mathbf{Z}}\left( q_t(x_t,z_t,\tau) | \mathbf{x},\mathbf{z} \right) = F_{Y_s|\mathbf{X},\mathbf{Z}}\left( q_s(x_s,z_s,\tau) | \mathbf{x},\mathbf{z} \right),
$$
where $\mathbf{x} = (x_1, \cdots , x_T)'$ and $\mathbf{z} = (z_1, \cdots , z_T)'$.
\cite{ishihara2020identification} demonstrates that this condition provides point identification when the support of $\mathbf{X}_i$ satisfies Assumption A.2.
In Section 3, we propose an estimation method based on this condition.
\begin{Remark}
From the proof of Proposition 1, we can identify $q_t$ from the conditional distribution of $Y_{it}|\mathbf{X}_i$.
This implies that we do not need to observe $(Y_{i1}, \cdots , Y_{iT})$ simultaneously and $q_t$ is identified from repeated cross-sections.
Even when we do not have panel data, we can sometimes observe $(Y_{it},\mathbf{X}_i)$ from repeated cross-sections.
For example, when $X_{it}$ is the minimum wage at time $t$ in the county where the unit $i$ lives, we can observe $(Y_{it},\mathbf{X}_i)$ from repeated cross-sections if we know the county where the unit $i$ lives.
Although the identification results do not require panel data, we need the existence of panel data in the estimation part.
Hence, in this study, we assume that panel data is obtained.
\end{Remark}
To illustrate our model, we consider the following two examples:
\begin{Example}[The CIC model]
In standard DID models, the support of $(X_{i1},X_{i2})$ becomes $\{(0,0),(0,1)\}$.
Then, $G_i \equiv \mathbf{1} \{ X_{i2} = 1 \}$ denotes an indicator for the treatment group.
In this setting, only individuals in group 1 in period 2 are treated.
Our model is then similar to the CIC model proposed by \cite{athey2006identification}.
Under the assumptions of Proposition 1, we can obtain
\begin{eqnarray}
F_{Y_2(0)|G=1}(y) &=& F_{Y_1|G=1}\left( F_{Y_1|G=0}^{-1}\left( F_{Y_2|G=0}(y) \right) \right), \label{Identification_AI1} \\
F_{Y_2(1)|G=0}(y) &=& F_{Y_1|G=0}\left( F_{Y_1|G=1}^{-1}\left( F_{Y_2|G=1}(y) \right) \right). \label{Identification_AI2}
\end{eqnarray}
\cite{athey2006identification} proved (\ref{Identification_AI1}) without the rank invariance assumption.
Hence, if the target parameter is the QTT, the rank invariance assumption is not required; whereas, if we focus on the QTE, the rank invariance assumption is required.
In this setting, Assumption 2 implies that $U_{i1}|G_i=g,\mathbf{Z}_i=\mathbf{z} \ \overset{d}{=} \ U_{i2}|G_i=g, \mathbf{Z}_i=\mathbf{z}$ for all $g$ and $\mathbf{z}$.
This allows the treatment and control groups to differ in terms of unobservable ability because it does not assume that $U_{it}|G_i=0,\mathbf{Z}_i=\mathbf{z} \ \overset{d}{=} \ U_{it}|G_i=1, \mathbf{Z}_i=\mathbf{z}$.
Hence, the rank stationarity assumption allows the treatment group to contain more high-ability people than the control group, implying that the conditional distribution of $Y_{i2}(0)|G_i = 1, \mathbf{Z}_i=\mathbf{z}$ may be different from that of $Y_{i2}(0)|G_i = 0, \mathbf{Z}_i=\mathbf{z}$.
\end{Example}
\begin{Example}[The TV effect on test scores]
Let $X_{it}$ denote the daily TV watching hours and $Y_{it}$ denote the test score of student $i$ in year $t$.
We assume that $Y_{it}$ can be written as
\[
Y_{it} = q_t(X_{it},Z_{it},U_{it}),
\]
where $Z_{it}$ is a vector of observed characteristics.
We assume that the unobserved factor $U_{it}$ can be decomposed into time-variant and time-invariant parts.
Let $U_{it} = U(A_i,\epsilon_{it})$, where $A_i$ and $\epsilon_{it}$ represent the students' ability and idiosyncratic shocks, respectively.
Here, Assumption 2 is satisfied when we have
\[
\epsilon_{it} | \mathbf{X}_i , \mathbf{Z}_i, A_i \ \overset{d}{=} \ \epsilon_{is} | \mathbf{X}_i , \mathbf{Z}_i, A_i.
\]
Hence, the rank stationarity assumption does not impose any restrictions on the dependence between unobserved ability and TV watching.
In this example, it is important to model nonlinear time effects; for example, if $q_t$ does not change over time, it follows from Assumption 2 that the conditional distribution of $Y_{it}$ is the same across time.
However, this is not plausible because the difficulty level of the test changes over time.
\end{Example}
\section{Estimation and inference}
In this section, we consider the estimation method of the QTE.
First, in Section 3.1, we review the minimum distance method proposed by \cite{ishihara2020identification} and show that the minimum distance estimator does not work when there are many covariates.
Second, in Section 3.2, we propose a two-step estimator based on the quantile regression and minimum distance methods and show that our estimator is computationally convenient.
Finally, in Sections 3.3 and 3.4, we demonstrate the consistency and uniform asymptotic normality of our estimator.
\subsection{The minimum distance estimator}
In this section, for simplicity, we assume that $T=2$ and there are no covariates.
\cite{ishihara2020identification} considers the following parametric model:
\[
q_t(x_t,\tau) = g_t(x_t,\tau ; \theta_0).
\]
The structural functions are parameterized by $\theta \in \Theta \subset \mathbb{R}^{d_{\theta}}$, where $\theta_0 \in \Theta$ is the true parameter.
Then, from Assumption 2, we obtain
\begin{eqnarray}
E\left[ \mathbf{1} \left\{ Y_{i1} \leq g_1(X_{i1},\tau ; \theta_0) \right\} | \mathbf{X}_i \right] = E\left[ \mathbf{1} \left\{ Y_{i2} \leq g_2(X_{i2},\tau ; \theta_0) \right\} | \mathbf{X}_i \right]. \label{conditional_moment_condition}
\end{eqnarray}
Thus, \cite{ishihara2020identification} proposes a minimum distance estimator based on (\ref{conditional_moment_condition}).
Let $\| \cdot \|_{\mu}$ denote the $L_2$-norm with respect to a probability measure $\mu$ with support $[0,1] \times \mathcal{V}$.
The minimum distance estimator $\hat{\theta}$ is then obtained from the following optimization:
\begin{eqnarray}
\hat{\theta} &=& \min_{\theta} \| \hat{D}_{\theta} \|_{\mu}, \label{MD_estimator} \\
\hat{D}_{\theta}(\tau,v) &\equiv & \frac{1}{n} \sum_{i=1}^n \left( \mathbf{1} \left\{ Y_{i1} \leq g_1(X_{i1},\tau ; \theta_0) \right\} - \mathbf{1} \left\{ Y_{i2} \leq g_2(X_{i2},\tau ; \theta) \right\} \right) \omega(\mathbf{X}_i,v), \nonumber
\end{eqnarray}
where $\omega(\mathbf{x},v)$ is a weight function.
Since $\hat{D}_{\theta}(\tau,v)$ is not continuous in $\theta$, the minimum distance estimator requires minimizing the discontinuous objective function over $\theta \in \Theta$.
If the dimension of $\theta$ is large, the optimization (\ref{MD_estimator}) is computationally demanding.
Therefore, adding many covariates into the model makes it difficult to compute the minimum distance estimator.
\subsection{A two-step estimator}
For estimation, we focus on the following linear-in-parameter model:
\begin{equation}
q_t(x_t,z_t,\tau) = x_t'\alpha(\tau) + z_t'\beta_t(\tau). \label{Model}
\end{equation}
The observed outcome is then written as
$$
Y_{it} = X_{it}'\alpha(U_{it})+Z_{it}'\beta_t(U_{it}), \ \ \ U_{it}|\mathbf{Z}_i \sim U(0,1),
$$
where $Z_{it}$ contains a constant term.
Hereafter, we set $\mathbf{Z}_i \in \mathbb{R}^{d_z}$ as a vector of all the variables of $Z_{i1}, \cdots, Z_{iT}$.
For example, if all covariates are time-invariant, we have $Z_{i1} = \cdots = Z_{iT} = \mathbf{Z}_i$.
We assume that $\mathcal{X}_{t}$ and $\mathcal{Z}$ are bounded.
In this model, we have $\partial q_t(x,z,\tau) / \partial x = \alpha(\tau)$; hence, our target parameter is $\alpha(\tau)$.
Because $\beta_t(\tau)$ depends on the time period, this model captures nonlinear time effects.
This model is similar to the IV quantile regression model proposed by \cite{chernozhukov2006instrumental}.
Using Proposition 1, we can identify $\alpha(\tau)$ and $\beta_t(\tau)$ using the following conditions:
\begin{eqnarray}
F_{Y_t|\mathbf{X},\mathbf{Z}}(x_t'\alpha(\tau)+z_t'\beta_t(\tau)|\mathbf{x},\mathbf{z}) &=& F_{Y_s|\mathbf{X},\mathbf{Z}}(x_s'\alpha(\tau)+z_s'\beta_s(\tau)|\mathbf{x},\mathbf{z}) \label{moment_1} \\
F_{Y_t-X_t'\alpha(\tau)|\mathbf{Z}}(z_t'\beta_t(\tau)|\mathbf{z}) &=& \tau, \label{moment_2}
\end{eqnarray}
where $\mathbf{x} \equiv (x_1, \cdots ,x_T)'$ and $\mathbf{z} \equiv (z_1, \cdots ,z_T)'$, respectively.
Similar to (\ref{MD_estimator}), we can construct a minimum distance estimator using (\ref{moment_1}) and (\ref{moment_2}).
However, if the dimensionality of covariates is high, the minimum distance approach cannot be directly applied because the minimum distance estimator is computationally demanding.
We propose the following two-step estimator based on the quantile regression and minimum distance methods.
Fix $\tau \in (0,1)$.
In the first step, we define $\tilde{\beta}_t(a,\tau)$ as
\begin{eqnarray}
\tilde{\beta}_t(a,\tau) & \equiv & \text{arg} \min_{b_t \in \mathcal{B}_t} \frac{1}{n} \sum_{i=1}^n \rho_{\tau}\left( Y_{it} - X_{it}'a - Z_{it}'b_t \right) \nonumber \\
&=& \text{arg} \min_{b_t \in \mathcal{B}_t} \frac{1}{n} \sum_{i=1}^n R_{\tau}(W_{it};a,b_t), \label{first step estimator}
\end{eqnarray}
where $\rho_{\tau}(u) \equiv (\tau - \mathbf{1}\{u<0\})u$, $\mathcal{B}_t$ is the parameter space of $\beta_t(\tau)$, and $R_{\tau}(W_{it};a,b_t) \equiv \rho_{\tau}\left( Y_{it} - X_{it}'a - Z_{it}'b_t \right)$.
This is an ordinary quantile regression of $Y_{it}-X_{it}'a$ on $Z_{it}$.
Then, from (\ref{moment_2}), $\tilde{\beta}_t \left( \alpha(\tau),\tau \right)$ becomes a consistent estimator of $\beta_t(\tau)$.
In the second step, we construct an estimator of $\alpha(\tau)$ using the minimum distance approach.
We define
\begin{eqnarray}
g_t(\mathbf{W}_i;a,b,v) &\equiv & \left( \mathbf{1}\{Y_{it} \leq X_{it}'a + Z_{it}'b_t\} - \frac{1}{T} \sum_{s=1}^T \mathbf{1}\{Y_{is} \leq X_{is}'a + Z_{is}'b_s\} \right) \omega(\mathbf{X}_i,\mathbf{Z}_i,v), \nonumber \\
\omega(\mathbf{X}_i,\mathbf{Z}_i,v) &\equiv & \exp\left(v_{\mathbf{x}}'\tilde{\mathbf{X}}_i + v_{\mathbf{z}}'\tilde{\mathbf{Z}}_{i} \right), \nonumber
\end{eqnarray}
where $b = (b_1', \cdots , b_T')'$, $v = (v_{\mathbf{x}}', v_{\mathbf{z}}')'$, and $\tilde{\mathbf{X}}_i$ and $\tilde{\mathbf{Z}}_{i}$ are standardized versions of $\mathbf{X}_i$ and $\mathbf{Z}_i$, where each component has a mean of 0 and a standard deviation of 1.
It follows from (\ref{moment_1}) that we have
\begin{equation}
E\left[ g_t(\mathbf{W}_i;\alpha(\tau),\beta(\tau),v) \right] = 0 \ \ \text{for all $t$ and $v$,} \label{moment_second_step}
\end{equation}
where $\beta(\tau) \equiv (\beta_1(\tau)', \cdots, \beta_T(\tau)')'$.
As shown in \cite{stinchcombe1998consistent}, if (\ref{moment_second_step}) holds for all $t$ and $v \in \mathcal{V} \equiv [-0.5,0.5]^{d_{X}\cdot T + d_{Z}}$, the conditional moment condition (\ref{moment_1}) is satisfied.
Let $\|\cdot\|_{L_2}$ be the $L_2$-norm over a compact set $\mathcal{V}$; that is, $\|f(v)\|_{L_2}^2 = \int_{\mathcal{V}} f(v)^2 dv$.
Using this norm, we obtain the following estimator of $\alpha(\tau)$:
\begin{eqnarray}
\hat{\alpha}(\tau) &\equiv & \text{arg} \min_{a \in \mathcal{A}} \frac{1}{T} \sum_{t=1}^T \left\| \frac{1}{n} \sum_{i=1}^n g_t(\mathbf{W}_i;a,\tilde{\beta}(a,\tau),v) \right\|_{L_2}^2 \nonumber \\
& = & \text{arg} \min_{a \in \mathcal{A}} \frac{1}{T} \sum_{t=1}^T \left\| \hat{D}^t_n(v;a,\tilde{\beta}(a,\tau)) \right\|_{L_2}^2, \label{Estimator alpha}
\end{eqnarray}
where $\tilde{\beta}(a,\tau) \equiv (\tilde{\beta}_1(a,\tau)', \cdots , \tilde{\beta}_T(a,\tau))'$, $\hat{D}^t_n(v;a,b)\equiv \frac{1}{n} \sum_{i=1}^n g_t(\mathbf{W}_i;a,b,v)$, and $\mathcal{A}$ is the parameter space of $\alpha(\tau)$.
Finally, we estimate $\beta_t(\tau)$ by $\hat{\beta}_t(\tau) \equiv \tilde{\beta}_t(\hat{\alpha}_t(\tau),\tau)$.
We briefly explain our two-step estimation method.
As discussed above, because the covariates are independent of $U_{it}$, $\tilde{\beta}_t \left( \alpha(\tau),\tau \right)$ becomes a consistent estimator of $\beta_t(\tau)$.
Using this result, under regularity conditions, we obtain
\[
\hat{D}_n^t \left( v ; \alpha(\tau), \tilde{\beta}(\alpha(\tau),\tau) \right) \rightarrow_p E[g_t(\mathbf{W}_i; \alpha(\tau), \beta(\tau),v)] = 0,
\]
which implies that the objective function of (\ref{Estimator alpha}) converges to zero for $a = \alpha(\tau)$.
Hence, we expect to obtain a consistent estimator of $\alpha(\tau)$ by minimizing (\ref{Estimator alpha}).
In practice, we can implement this estimation procedure as follows:
\begin{enumerate}
\item For fixed $\tau$, we run the ordinary $\tau$-th quantile regression of $Y_{it}-X_{it}'a$ on $Z_{it}$ and calculate $\tilde{\beta}_t(a,\tau)$ as a function of $a$.
\item We approximate $\| \hat{D}^t_n(v;a,\tilde{\beta}(a,\tau)) \|_{L_2}^2$ using a numerical integration method; that is, we approximate the objective function as $\frac{1}{J} \sum_{j = 1}^J \hat{D}^t_n(v_j;a,\tilde{\beta}(a,\tau))^2$ for an appropriate sequence $\{v_j\}_{j=1}^J$.
\item Minimize $\frac{1}{T} \sum_{t=1}^T \| \hat{D}^t_n(v;a,\tilde{\beta}(a,\tau)) \|_{L_2}^2$ over $a \in \mathcal{A}$ and obtain $\hat{\alpha}(\tau)$.
The estimate $\hat{\beta}_t(\tau)$ is given by $\tilde{\beta}_t(\hat{\alpha}(\tau),\tau)$.
\end{enumerate}
Our estimator is similar to that proposed by \cite{chernozhukov2006instrumental}.
They consider the IV quantile regression for heterogeneous treatment effect models and simultaneous equation models with nonadditive errors.
Similarly, our estimator is attractive from a computational point of view.
As ordinary quantile regressions are obtained by convex optimization, our first step estimation (\ref{first step estimator}) is computationally convenient.
Our second step estimation (\ref{Estimator alpha}) requires non-convex optimization; hence, it seems to be computationally demanding.
However, we can obtain (\ref{Estimator alpha}) by optimizing the objective function over the $\alpha$ parameter (typically one-dimensional).
This fact makes our estimator computationally convenient.
\begin{Remark}
Although we assume that $\alpha(\tau)$ does not depend on the time period, we can relax this assumption.
Even if the QTE parameter is $\alpha_t(\tau)$, we can estimate $\alpha_t(\tau)$ in a similar manner.
However, in such a case, our second step estimation requires non-convex optimization with respect to $a_1, \cdots, a_T$.
Hence, if $T$ is large, this estimation method becomes computationally demanding.
\end{Remark}
\subsection{Identification}
In this section, we show that $\alpha(\tau)$ and $\beta(\tau) = (\beta_1(\tau)', \cdots , \beta_T(\tau)')'$ uniquely solve the limit problems.
We define
\begin{equation}
\beta_t(a,\tau) \equiv \text{arg} \min_{b_t \in \mathcal{B}_t} E[ R_{\tau}(W_{it};a,b_t)], \label{beta limit problem}
\end{equation}
and
\begin{equation}
\alpha^*(\tau) \in \text{arg} \min_{a \in \mathcal{A}} \frac{1}{T} \sum_{t=1}^T \left\|D^t(v;a,\beta(a,\tau)) \right\|_{L_2}^2, \nonumber
\end{equation}
where $\beta(a,\tau) \equiv (\beta_1(a,\tau)', \cdots , \beta_T(a,\tau)')'$ and $D^t(v;a,b) \equiv E[g_t(\mathbf{W}_{i};a,b,v)]$.
Hence, to prove consistency, we need to show that $\alpha^*(\tau)$ is unique and $\alpha^*(\tau) = \alpha(\tau)$.
We define $e_t(a,\tau,\mathbf{z}) \equiv P ( Y_{it} \leq X_{it}'a + Z_{it}'\beta_t(a,\tau) | \mathbf{Z}_i=\mathbf{z} )$ and impose the following assumptions.
\begin{Assumption}
(i) A matrix $E[Z_{it}Z_{it}']$ has full rank for all $t$, and $E[X_{it}X_{it}']$ has full rank for some $t$.
(ii) For all $t$, $E[|Y_{it}|]$ is finite.
(iii) For all $t$, $a$, and $\tau$, $\beta_t(a,\tau)$ uniquely solves (\ref{beta limit problem}).
\end{Assumption}
\begin{Assumption}
For all $t$ and $a \in \mathcal{A}$, $e_t(a,\tau,\mathbf{z}) = \tau$ for some $\mathbf{z} \in \mathcal{Z}$.
\end{Assumption}
When the support of $(X_{i1},X_{i2})$ is $\{(0,0),(0,1)\}$, we have $E[X_{i1}^2] = 0$ but $E[X_{i2}^2]$ is positive.
Hence, Assumption 3 (i) holds in standard DID settings.
Assumption 4 is a technical condition that is satisfied in many situations.
Using the proof of Theorem 2 in \cite{angrist2006quantile}, it follows from the first-order condition of (\ref{first step estimator}) that we have $E\left[ \left( \mathbf{1}\{Y_{it} \leq X_{it}'a + Z_{it}'\beta_t(a,\tau)\} - \tau \right) Z_{it} \right]=0$, which implies that $E\left[ \left( e_t(a,\tau,\mathbf{Z}_i) - \tau \right) Z_{it} \right]=0$.
Hence, we have $E[e_t(a,\tau,\mathbf{Z}_{i})]=\tau$ because $Z_{it}$ contains a constant.
When $\mathbf{Z}_i$ has continuous covariates and $e_t(a,\tau,\mathbf{z})$ is continuous in $\mathbf{z}$, $e_t(a,\tau,\mathbf{z}) = \tau$ holds for some $\mathbf{z} \in \mathcal{Z}$.
Even when all covariates are discrete, if $Z_{it}$ is time invariant and the model is saturated, that is, the cardinality of $\mathcal{Z}$ is equal to the dimension of $\beta_t(a,\tau)$, then we have $e_t(a,\tau,\mathbf{z})= \tau$ for all $\mathbf{z} \in \mathcal{Z}$.
\begin{Theorem}
Suppose that (\ref{Model}) and Assumptions 1--4, A.1, and A.2 hold.
Then, for all $\tau \in (0,1)$, $\alpha(\tau)$ and $\beta(\tau)$ uniquely solve the limit problems.
That is, we have $\beta_t(\alpha(\tau),\tau) = \beta_t(\tau)$ and
\begin{equation}
\frac{1}{T} \sum_{t=1}^T \left\| D^t(v;a,\beta(a,\tau)) \right\|_{L_2}^2=0, \ \ a \in \mathcal{A} \ \ \Leftrightarrow \ \ a = \alpha(\tau). \label{Identification Estimator}
\end{equation}
\end{Theorem}
Theorem 1 implies that $\alpha(\tau)$ minimizes $\frac{1}{T} \sum_{t=1}^T \left\|D^t(v;a,\beta(a,\tau)) \right\|_{L_2}^2$.
Hence, if the objective function of (\ref{Estimator alpha}) converges to $\frac{1}{T} \sum_{t=1}^T \left\|D^t(v;a,\beta(a,\tau)) \right\|_{L_2}^2$ uniformly, then we obtain the consistency of $\hat{\alpha}(\tau)$.
\subsection{Asymptotic distribution}
In this section, we show the uniform asymptotic normality of our estimator and prove the validity of the nonparametric bootstrap.
Our asymptotic result also implies that our estimator is consistent.
Let $\mathcal{T}$ be a closed subset of $[\epsilon, 1- \epsilon]$ for $\epsilon > 0$.
In addition, we define $J_t^b(a,\tau) \equiv E\left[ f_{Y_t-X_t'a|Z_t}(Z_{it}'\beta_t(a,\tau)|Z_{it}) Z_{it} Z_{it}' \right]$ and $J_t^b(\tau) \equiv J_t^b(\alpha(\tau),\tau)$.
The following assumption is sufficient for the consistency of $\hat{\alpha}(\tau)$ and $\hat{\beta}_t(\tau)$.
\begin{Assumption}
(i) The data $\{\mathbf{W}_i\}_{i=1}^n$ are independent and identically distributed.
(ii) For all $\tau \in \mathcal{T}$, $\alpha(\tau)$ and $\beta_{t}(\tau)$ are contained in the compact parameter spaces $\mathcal{A}$ and $\mathcal{B}_t$, respectively.
(iii) For all $t$, $E[|Y_{it}|]$ is finite, and $\mathcal{X}_{1:T}$ and $\mathcal{Z}$ are bounded.
(iv) For all $a$, $t$, and $\tau$, $\beta_t(a,\tau)$ uniquely solves (\ref{beta limit problem}).
(v) For all $a$ and $t$, the conditional density $f_{Y_t-X_t'a|Z_t}(y|z_t)$ exists and $f_{Y_t-X_t'a|Z_t}(y|z_t)$ is continuous in $y$ and bounded above.
(vi) For all $t$ and $\tau$, $J_t^b(a,\tau)$ has full rank for all $a \in \mathcal{A}$ and $\tau \in \mathcal{T}$, and $J_t^b(a,\tau)$ is continuous in $a$ at $\alpha(\tau)$.
(vii) For all $t$, $F_{Y_t|\mathbf{X},\mathbf{Z}}(y|\mathbf{x},\mathbf{z})$ is uniformly continuous in $y$.
(viii) For all $t$, $\beta_t(a,\tau)$ is continuously differentiable in $a$ and its derivative $B_t(a,\tau) \equiv \frac{\partial}{\partial a'} \beta_t(a,\tau)$ is bounded.
\end{Assumption}
Condition (iii) imposes that $X_{it}$ and $Z_{it}$ are bounded.
If $X_{it}$ and $Z_{it}$ are unbounded, then for $q_t(x_t,z_t,\tau)$ to be monotonically increasing in $\tau$, $\alpha(\tau) = \alpha$ and $\beta_t(\tau) = \beta_t$ must hold.
Hence, we assume the boundedness of $\mathcal{X}_{1:T}$ and $\mathcal{Z}$.
Condition (v) means that there exists a continuous density $f_{Y_t-X_t'a|Z_t}(y|z_t)$ for all $a \in \mathcal{A}$.
Because we have
\begin{eqnarray*}
F_{Y_t-X_t'a|Z_t}(y|z_t) &=& \int F_{Y_t|X_t,Z_t}(y + x_t'a |x_t, z_t) dF_{X_t}(x_t),
\end{eqnarray*}
we obtain $f_{Y_t-X_t'a|Z_t}(y|z_t) = \int f_{Y_t|X_t,Z_t}(y + x_t'a |x_t, z_t) dF_{X_t}(x_t)$ if $f_{Y_t|X_t,Z_t}(y|x_t,z_t)$ is bounded.
Hence, condition (v) holds if $f_{Y_t|X_t,Z_t}(y|x_t,z_t)$ is bounded and continuous in $y$.
In addition, this implies that $J_t^b(a,\tau)$ is continuous in $a$ if $\beta_t(a,\tau)$ is continuous in $a$.
We define
\begin{eqnarray}
\gamma_1^t(v;a,\tau) &\equiv & E \left[ f_{Y_t|\mathbf{X},\mathbf{Z}}(X_{it}'a + Z_{it}'\beta_t(a,\tau)|\mathbf{X}_i,\mathbf{Z}_i) \omega(\mathbf{X}_i,\mathbf{Z}_i,v) (X_{it} + B_t(a,\tau)'Z_{it}) \right] \nonumber \\
\Gamma_1^t(v;a,\tau) &\equiv & \gamma_1^t(v;a,\tau) - \frac{1}{T} \sum_{s=1}^T \gamma_1^s(v;a,\tau) \nonumber \\
\gamma_2^{t,s}(v;a,b) &\equiv & \begin{cases}
\frac{T-1}{T} E \left[ f_{Y_t|\mathbf{X},\mathbf{Z}}(X_{it}'a + Z_{it}'b_t|\mathbf{X}_i,\mathbf{Z}_i) \omega(\mathbf{X}_i,\mathbf{Z}_i,v) Z_{it} \right], & \text{if $s=t$} \\
-\frac{1}{T} E \left[ f_{Y_s|\mathbf{X},\mathbf{Z}}(X_{is}'a + Z_{is}'b_s|\mathbf{X}_i,\mathbf{Z}_i) \omega(\mathbf{X}_i,\mathbf{Z}_i,v) Z_{is} \right], & \text{if $s \neq t$}
\end{cases}, \nonumber \\
\Gamma_2^t(v;a,b) &\equiv & \left( \gamma_2^{t,1}(v;a,b)' , \cdots , \gamma_2^{t,T}(v;a,b)' \right)', \nonumber
\end{eqnarray}
$\Gamma_1^t(v;\tau) \equiv \Gamma_1^t(v;\alpha(\tau),\tau)$, and $\Gamma_2^t(v;\tau) \equiv \Gamma_2^t(v;\alpha(\tau),\beta(\tau))$.
Then, the following assumption is required to derive the asymptotic distribution of the estimator.
\begin{Assumption}
(i) For all $a$, $t$, and $\tau$, $\alpha(\tau)$ and $\beta_t(a,\tau)$ are the inner points of $\mathcal{A}$ and $\mathcal{B}_t$, respectively.
(ii) A family of functions $\{y \mapsto f_{Y_t-X_t'a|Z_t}(y|z): a \in \mathcal{A}\}$ is equicontinuous for all $z \in \mathcal{Z}$.
(iii) There exists the conditional density $f_{Y_t|\mathbf{X},\mathbf{Z}}(y|\mathbf{x},\mathbf{z})$ and $f_{Y_t|\mathbf{X},\mathbf{Z}}(y|\mathbf{x},\mathbf{z})$ is uniformly continuous in $y$ and bounded above.
(vi) For all $t$, a family of functions $\{a \mapsto B_t(a,\tau): \tau \in \mathcal{T}\}$ is equicontinuous.
(v) There exists $c>0$ such that $T^{-1} \sum_{t=1}^T\|\Gamma_1^t(v;\tau)'a\|_{L_2}^2 \geq c^2 \|a\|^2$ for all $a \in \mathbb{R}^{d_X}$ and $\tau \in \mathcal{T}$.
\end{Assumption}
To derive the asymptotic distribution of $\hat{\alpha}(\tau)$, we need to show that $\sqrt{n} (\tilde{\beta}_t(a,\tau) - \beta_t(a, \tau) )$ converges in distribution uniformly in $a \in \mathcal{A}$.
We use condition (ii) to show this result.
Similarly, we need condition (iv) to show a uniform approximation of $\hat{\alpha}(\cdot)$.
Condition (v) indicates that the rank condition holds uniformly in $\tau \in \mathcal{T}$.
\begin{Theorem}
Suppose that (\ref{Model}) and (\ref{Identification Estimator}) hold for all $\tau \in \mathcal{T}$.
Under Assumptions 5 and 6, uniformly in $\tau \in \mathcal{T}$ we obtain
\begin{eqnarray}
\sqrt{n} (\hat{\alpha}(\tau)-\alpha(\tau)) &=& - \frac{1}{\sqrt{n}} \sum_{i=1}^n \mathbb{A}(\mathbf{W}_i; \tau) + o_p(1) \label{Asymptotic Normality}
\end{eqnarray}
and
\begin{equation}
\sqrt{n}(\hat{\beta}_t(\tau)-\beta_t(\tau)) = - \frac{1}{\sqrt{n}} \sum_{i=1}^n \left\{ J_t^b(\tau)^{-1} r_{\tau}(W_{it};\alpha(\tau),\beta_t(\tau)) - \mathbb{A}(\mathbf{W}_i; \tau) \right\} + o_p(1), \label{Asymptotic Normality first step estimator}
\end{equation}
where
\begin{eqnarray}
\mathbb{A}(\mathbf{W}_i; \tau) &\equiv & \Delta_1(\tau)^{-1} \left\{ \xi(\mathbf{W}_i;\tau) - \Delta_{12}(\tau) l(\mathbf{W}_i;\tau) \right\}, \nonumber \\
\xi(\mathbf{W}_i;\tau) &\equiv & \frac{1}{T} \sum_{t=1}^T \left[ \int_{\mathcal{V}} \Gamma_1^t(v;\tau) \omega(\mathbf{X}_i,\mathbf{Z}_i,v) dv \right] \mathbf{1}\{Y_{it} \leq X_{it}'\alpha(\tau)+Z_{it}'\beta_t(\tau)\}, \nonumber \\
l(\mathbf{W}_i;\tau) &\equiv & \left( r_{\tau}(W_{i1};\alpha(\tau),\beta_t(\tau))' J_1^b(\tau)^{-1}, \cdots, r_{\tau}(W_{iT};\alpha(\tau),\beta_t(\tau))'J_T^b(\tau)^{-1} \right)', \nonumber \\
r_{\tau}(W_{it};a,b_t) &\equiv & \left( \tau - \mathbf{1}\{Y_{it} \leq X_{it}'a + Z_{it}'b_t\} \right)Z_{it}, \nonumber \\
\Delta_1(\tau) &\equiv & \frac{1}{T} \sum_{t=1}^T \int_{\mathcal{V}} \Gamma_1^t(v;\tau)\Gamma_1^t(v;\tau)' dv, \nonumber \\
\Delta_{12}(\tau) &\equiv & \frac{1}{T} \sum_{t=1}^T \int_{\mathcal{V}} \Gamma_1^t(v;\tau)\Gamma_2^t(v;\tau)' dv. \nonumber
\end{eqnarray}
\end{Theorem}
The proof of this theorem is based on arguments similar to those in \cite{brown2002weighted}, \cite{chen2003estimation}, and \cite{torgovitsky2017minimum}.
The following corollary follows immediately from Theorem 2.
\begin{Corollary}
We define $\Sigma(\tau,\tilde{\tau}) \equiv E[\mathbb{A}(\mathbf{W}_i; \tau) \mathbb{A}(\mathbf{W}_i; \tilde{\tau})']$.
Under the assumptions of Theorem 2, $\sqrt{n} (\hat{\alpha}(\cdot) - \alpha(\cdot) )$ converges weakly to a zero mean Gaussian process $\mathbb{Z}(\cdot)$ with covariance function $\Sigma(\tau,\tilde{\tau})$.
\end{Corollary}
We consider the case in which $T=2$, $X_{it}$ is scalar, and there are no covariates.
In this case, we have
\begin{eqnarray}
\xi(\mathbf{W}_i;\tau) &=& \frac{1}{4} \left[ \int_{\mathcal{V}} \gamma_1(v;\tau) \omega(\mathbf{X}_i,v) dv \right] (\mathbf{1}\{U_{i1} \leq \tau\} - \mathbf{1}\{U_{i2} \leq \tau\} ), \nonumber \\
\Delta_{12}(\tau) l(\mathbf{W}_i;\tau) &=& \frac{1}{4} \left[ \int_{\mathcal{V}} \gamma_1(v;\tau) \gamma_2^1(v;\tau) dv \right] J_1^b(\tau)^{-1}( \tau - \mathbf{1}\{U_{i1} \leq \tau\} ) \nonumber \\
& & - \frac{1}{4} \left[ \int_{\mathcal{V}} \gamma_1(v;\tau) \gamma_2^2(v;\tau) dv \right] J_2^b(\tau)^{-1}( \tau - \mathbf{1}\{U_{i2} \leq \tau\} ), \nonumber
\end{eqnarray}
where $Y_{it}(\tau) \equiv \alpha(\tau) X_{it} + \beta_t(\tau)$,
\begin{eqnarray}
\gamma_1(v;\tau) &\equiv & E\left[ \left\{ f_{Y_1|\mathbf{X}}(Y_{i1}(\tau)|\mathbf{X}_i) (X_{i1}+B_1(\alpha(\tau),\tau)) \right. \right. \nonumber \\
& & \left. \left. \ \ \ \ \ - f_{Y_2|\mathbf{X}}(Y_{i2}(\tau)|\mathbf{X}_i)(X_{i2}+B_2(\alpha(\tau),\tau)) \right\} \omega(\mathbf{X}_i,v) \right], \nonumber
\end{eqnarray}
and $\gamma_2^t(v;\tau) \equiv E\left[ f_{Y_t|\mathbf{X}}(Y_{it}(\tau)|\mathbf{X}_i)\omega(\mathbf{X}_i,v) \right]$.
Because $J_t^b(\tau)$, $\gamma_2^1(v;\tau)$, and $\gamma_2^2(v;\tau)$ are positive, the variances of $\xi(\mathbf{W}_i;\tau)$ and $\Delta_{12}(\tau) l(\mathbf{W}_i;\tau)$ become small when $U_{i1}$ and $U_{i2}$ are positively correlated.
Specifically, if $U_{i1}=U_{i2}$, then $\xi(\mathbf{W}_i;\tau)$ is exactly equal to zero.
\if0
\begin{Remark}[Standard errors]
As seen in Corollary 1, the asymptotic variance of $\hat{\alpha}(\tau)$ is $\Sigma(\tau,\tau)$.
To estimate $\Sigma(\tau,\tau)$, we define
\begin{eqnarray}
\hat{\xi}(\mathbf{W}_i;\tau) &\equiv & \frac{1}{T} \sum_{t=1}^T \left[ \int_{\mathcal{V}} \hat{\Gamma}_1^t(v;\tau) \omega(\mathbf{X}_i,\mathbf{Z}_i,v) dv \right] \mathbf{1}\{Y_{it} \leq X_{it}'\hat{\alpha}(\tau)+Z_{it}'\hat{\beta}_t(\tau)\}, \nonumber \\
\hat{l}(\mathbf{W}_i;\tau) &\equiv & \left( \hat{J}_1^b(\tau)^{-1}r_{\tau}(W_{i1};\hat{\alpha}(\tau),\hat{\beta}_t(\tau))', \cdots, \hat{J}_T^b(\tau)^{-1}r_{\tau}(W_{iT};\hat{\alpha}(\tau),\hat{\beta}_t(\tau))' \right)', \nonumber \\
\hat{\Delta}_1(\tau) &\equiv & \frac{1}{T} \sum_{t=1}^T \int_{\mathcal{V}} \hat{\Gamma}_1^t(v;\tau) \hat{\Gamma}_1^t(v;\tau)' dv, \nonumber \\
\hat{\Delta}_{12}(\tau) &\equiv & \frac{1}{T} \sum_{t=1}^T \int_{\mathcal{V}} \hat{\Gamma}_1^t(v;\tau)\hat{\Gamma}_2^t(v;\tau)' dv. \nonumber
\end{eqnarray}
Here, $\hat{\Gamma}_1^t(v;\tau)$, $\hat{\Gamma}_2^t(v;\tau)$, and $\hat{J}_t^b(\tau)$ are uniformly consistent estimates of $\Gamma_1^t(v;\tau)$, $\Gamma_2^t(v;\tau)$, and $J_t^b(\tau)$, respectively.
These estimates are provided in Appendix 3.
Under some regularity conditions, we obtain $\hat{\Sigma}(\tau,\tau) \equiv \frac{1}{n} \sum_{i=1}^n \hat{\mathbb{A}}(\mathbf{W}_i;\tau) \hat{\mathbb{A}}(\mathbf{W}_i;\tau)' \to_p \Sigma(\tau,\tau)$, where $\hat{\mathbb{A}}(\mathbf{W}_i;\tau) \equiv \hat{\Delta}_1(\tau)^{-1} \left\{ \hat{\xi}(\mathbf{W}_i;\tau) - \hat{\Delta}_{12}(\tau) \hat{l}(\mathbf{W}_i;\tau) \right\}$.
\end{Remark}
\fi
Let $\{\mathbf{W}_{i}^*\}_{i=1}^n$ denote a bootstrap sample drawn with replacement from $\{\mathbf{W}_i\}_{i=1}^n$.
That is, $\{\mathbf{W}_{i}^*\}_{i=1}^n$ are independently and
identically distributed from the empirical measure, conditional on the realizations $\{\mathbf{W}_i\}_{i=1}^n$.
We define $\hat{\alpha}^*(\tau)$ as the bootstrap counterpart to $\hat{\alpha}(\tau)$.
Then, we can obtain the following theorem.
\begin{Theorem}
Under the assumptions of Theorem 2, $\sqrt{n}(\hat{\alpha}^*(\cdot) - \hat{\alpha}(\cdot))$ converges weakly to the limit distribution of $\sqrt{n}(\hat{\alpha}(\cdot) - \alpha(\cdot))$ in probability.
\end{Theorem}
\begin{Remark}
Using Theorem 3, we can consider the following null hypothesis:
\begin{equation}
H_0: \, \alpha(\tau) = r(\tau) \ \ \text{for each $\tau \in \mathcal{T}$,}
\end{equation}
where $r(\cdot)$ is known or estimable.
Then, we can use the following test statistic:
\[
S_n \ \equiv \ n \int_{\mathcal{T}} \left\| \hat{\alpha}(\tau) - \hat{r}(\tau) \right\|^2 d \tau.
\]
In practice, we approximate this integration by using a grid $\mathcal{T}_n$ in place of $\mathcal{T}$.
If the null hypothesis is that $\exists \alpha, \, \alpha(\tau) = \alpha$ for each $\tau \in \mathcal{T}$, $r(\tau)$ can be estimated by $\hat{r}(\tau) = \hat{\alpha}(0.5)$ under the null hypothesis.
In this case, Theorem 3 implies that the critical value can be calculated using the quantile of
\[
S_n^* \ \equiv \ n \int_{\mathcal{T}} \left\| (\hat{\alpha}^*(\tau) - \hat{\alpha}(\tau)) - (\hat{\alpha}^*(0.5) - \hat{\alpha}(0.5)) \right\|^2 d\tau.
\]
\end{Remark}
\section{Simulations}
\textbf{Simulation 1.}\,
Suppose that the potential outcomes are given by
\begin{eqnarray}
Y_{i1}(x) &=& \left(1 + 0.5 \Phi^{-1}(U_{i1}) \right) x + \Phi^{-1}(U_{i1}) + Z_{i}, \nonumber \\
Y_{i2}(x) &=& \left(1 + 0.5 \Phi^{-1}(U_{i2}) \right) x + 1.2 \Phi^{-1}(U_{i2}) + 1.2 Z_{i}, \nonumber
\end{eqnarray}
where $Z_{i} \sim U(0,1)$ and $\Phi$ is the standard normal distribution function.
The observed outcomes are generated from $Y_{it}=Y_{it}(X_{it})$.
We assume that $X_{it} = \Phi(\tilde{X}_{it})$, $U_{it} = \Phi(A_i + \tilde{U}_{it})$, $(\tilde{X}_{i1},\tilde{X}_{i2},A_i)' \sim N(0,\Sigma_{\mathbf{X} A})$, and $\tilde{U}_{it} \sim N(0,1-\rho^2)$, where $\rho \in [0,1]$ and
\begin{eqnarray}
\Sigma_{\mathbf{X} A} = \left(
\begin{array}{ccc}
1 & 0.5 & 0.5\rho \\
0.5 & 1 & 0.5\rho \\
0.5\rho & 0.5\rho & \rho^2
\end{array}
\right). \nonumber
\end{eqnarray}
Then, $U_{it}$ is uniformly distributed and $\rho$ represents the dependence between $U_{i1}$ and $U_{i2}$.
When $\rho = 0$, $U_{i1}$ and $U_{i2}$ are uncorrelated and when $\rho = 1$, $U_{i1}$ and $U_{i2}$ are perfectly correlated.
Here, we have $\alpha(0.25) = 0.66$, $\alpha(0.5) = 1$, and $\alpha(0.75) = 1.34$.
Table 1 contains the results of this experiment for two different choices of the sample size, $1000$ and $2000$, and three different choices of $\rho^2$, $0.1$, $0.5$, and $0.9$.
The number of replications is set at $1000$ throughout.
Table 1 shows the bias, standard deviation, and MSE of the estimates of $\alpha(\tau)$ for $\tau = 0.25, 0.5$, and $0.75$.
For all settings, the bias is quite small.
Table 1 shows that the standard deviation and MSE decrease in all experiments as the sample size increases.
As expected, when the correlation between $U_{i1}$ and $U_{i2}$ is high (i.e. $\rho^2 = 0.9$), the standard deviation decreases.
\begin{table}[H]
\begin{center}
\caption{Results of Simulation 1}
\begin{tabular}{c c r r r r r r} \hline
& & \multicolumn{3}{c}{$n=1000$} & \multicolumn{3}{c}{$n=2000$} \\ \hline
& & $\rho^2=0.1$ & $\rho^2=0.5$ & $\rho^2=0.9$ & $\rho^2=0.1$ & $\rho^2=0.5$ & $\rho^2=0.9$ \\ \hline \hline
& bias & -0.017 & -0.018 & -0.013 & -0.003 & -0.015 & -0.014 \\
$\tau=0.25$ & std & 0.232 & 0.236 & 0.178 & 0.157 & 0.151 & 0.112 \\
& mse & 0.054 & 0.056 & 0.032 & 0.025 & 0.023 & 0.013 \\ \hline
& bias & -0.006 & -0.011 & -0.017 & 0.001 & -0.008 & -0.008 \\
$\tau=0.50$ & std & 0.204 & 0.204 & 0.149 & 0.140 & 0.133 & 0.099 \\
& mse & 0.042 & 0.042 & 0.022 & 0.020 & 0.018 & 0.010 \\ \hline
& bias & -0.013 & -0.018 & -0.018 & -0.002 & -0.011 & -0.013 \\
$\tau=0.75$ & std & 0.232 & 0.230 & 0.177 & 0.156 & 0.156 & 0.118 \\
& mse & 0.054 & 0.053 & 0.032 & 0.024 & 0.024 & 0.014 \\ \hline
\end{tabular}
\end{center}
\end{table}
We also verify that the nonparametric bootstrap procedure works for $(n,\rho^2) = (2000,0.9)$.
We calculate 90\% and 95\% confidence intervals of $\alpha(\tau)$ to obtain the coverage probabilities for $\tau = 0.25, 0.5$, and $0.75$.
Table 2 shows the nominal and actual coverage probabilities are close in all settings.
\begin{table}[H]
\begin{center}
\caption{Coverage probabilities of Simulation 1}
\begin{tabular}{c c c c} \hline
& $\tau = 0.25$ & $\tau = 0.50$ & $\tau = 0.75$ \\ \hline \hline
90\% & 0.894 & 0.898 & 0.892 \\
95\% & 0.938 & 0.946 & 0.940 \\ \hline
\end{tabular}
\end{center}
\end{table}
\noindent
\textbf{Simulation 2.}\,
To compare our estimation method with that of \cite{athey2006identification}, we consider the following model.
We assume that $\mathcal{X}_{1:2} = \{(0,0),(0,1)\}$ and the potential outcomes are given by
\begin{eqnarray}
Y_{i1}(0) &=& \Phi^{-1}(U_{i1}), \nonumber \\
Y_{i2}(x) &=& \left( 1+0.5\Phi^{-1}(U_{i2}) \right)x + 0.5 \Phi^{-1}(U_{i2}). \nonumber
\end{eqnarray}
The observed outcomes are generated from $Y_{i1} = Y_{i1}(0)$ and $Y_{i2}=Y_{i2}(X_{i2})$.
It is assumed that $X_{i2}= \mathbf{1}\{\tilde{X}_i+A_i \geq 0\}$, $U_{it} = \Phi(A_i + \tilde{U}_{it})$, $\tilde{X}_i \sim N(0,1)$, $A_i \sim N(0,\rho^2)$, and $\tilde{U}_{it} \sim N(0,1-\rho^2)$, where $\rho \in [0,1]$.
Then, $G_i \equiv \mathbf{1}\{X_{i2}=1\}$ denotes an indicator for the treatment group.
Since $U_{it}$ satisfies the rank stationarity assumption, we have
\[
E[Y_{i2}(0)-Y_{i1}(0)|G_i = g] \ = \ - 0.5 E\left[ \Phi(U_{i2}) | G_i =g \right].
\]
The conditional distribution of $U_{i2}|G_i=0$ is different from that of $U_{i2}|G_i=1$; therefore, this model does not satisfy the parallel trend assumption employed in standard DID models and we cannot estimate the average treatment effect on the treated (ATT) using the standard DID estimation method.
By contrast, using our estimation method, we can estimate the quantile functions of the potential outcomes and obtain an estimate of the ATT.
From (\ref{Identification_AI1}) and (\ref{Identification_AI2}), we can estimate $F_{Y_2(0)|G=1}(y)$ and $F_{Y_2(1)|G=0}(y)$ by
\begin{eqnarray}
\hat{F}_{Y_2(0)|G=1}(y) &\equiv & \hat{F}_{Y_1|G=1}\left( \hat{F}_{Y_1|G=0}^{-1}\left( \hat{F}_{Y_2|G=0}(y) \right) \right) \ \text{and} \nonumber \\
\hat{F}_{Y_2(1)|G=0}(y) &\equiv & \hat{F}_{Y_1|G=0}\left( \hat{F}_{Y_1|G=1}^{-1}\left( \hat{F}_{Y_2|G=1}(y) \right) \right), \nonumber
\end{eqnarray}
where $\hat{F}_{Y_t|G=g}(\cdot)$ and $\hat{F}_{Y_t|G=g}^{-1}(\cdot)$ are the empirical distribution and quantile functions, respectively.
The marginal distributions of the potential outcomes and QTE can be obtained using these estimators and the empirical distributions of $Y_{i2}|G_i=0$ and $Y_{i2}|G_i=1$.
We refer to this estimator as AI estimator.
Table 3 presents the bias, standard deviation, and MSE of our estimator and the AI estimator for $n = 500$ and three different choices of $\rho^2$, $0.1$, $0.5$, and $0.9$.
For all settings, the results of our estimator are similar to those of the AI estimator.
Hence, when there are no covariates, our estimator is not worse than the AI estimator.
\begin{table}[H]
\begin{center}
\caption{Results of Simulation 2}
\begin{tabular}{c c r r r r r r} \hline
& & \multicolumn{3}{c}{Our estimator} & \multicolumn{3}{c}{AI estimator} \\ \hline
& & $\rho^2=0.1$ & $\rho^2=0.5$ & $\rho^2=0.9$ & $\rho^2=0.1$ & $\rho^2=0.5$ & $\rho^2=0.9$ \\ \hline \hline
& bias & -0.011 & -0.026 & -0.015 & -0.003 & -0.012 & -0.027 \\
$\tau=0.25$ & std & 0.129 & 0.129 & 0.097 & 0.136 & 0.129 & 0.108 \\
& mse & 0.017 & 0.017 & 0.010 & 0.018 & 0.017 & 0.012 \\ \hline
& bias & -0.002 & -0.010 & -0.006 & 0.002 & -0.001 & -0.004 \\
$\tau=0.50$ & std & 0.116 & 0.100 & 0.069 & 0.115 & 0.097 & 0.065 \\
& mse & 0.013 & 0.010 & 0.005 & 0.013 & 0.009 & 0.004 \\ \hline
& bias & -0.005 & -0.017 & -0.010 & -0.003 & -0.004 & -0.006 \\
$\tau=0.75$ & std & 0.123 & 0.106 & 0.077 & 0.125 & 0.110 & 0.073 \\
& mse & 0.015 & 0.012 & 0.006 & 0.016 & 0.012 & 0.005 \\ \hline
\end{tabular}
\end{center}
\end{table}
\section{Empirical illustrations}
\subsection{The impact of insurance provision on household production}
In this section, we use our method to study the impact of an agricultural insurance program on household production.
We use the data employed by \cite{cai2016impact} to estimate the QTE of insurance provision on tobacco production.
This empirical analysis is based on data obtained from 12 tobacco production counties in the Jiangxi province of China.
Across these 12 counties, only tobacco farmers in the county of Guangchang were eligible to buy the tobacco insurance policy.
In 2003, the People's Insurance Company of China (PICC) designed and offered the first tobacco production insurance program to households in Guangchang.
Hence, we use this county as a treatment group.
The sample includes information on approximately 3,400 tobacco households during 2002 and 2003.
Table 4 provides summary statistics for 2002 and shows that treatment regions are quite different from control regions in terms of their observed characteristics.
For example, control regions include more educated people than treatment regions.
The proportion of high school- or college-educated people in the treatment regions is $0.025$, whereas that in the control regions is $0.257$.
Hence, controlling the observed characteristics is important for adjusting the differences between the treatment and control regions.
\begin{table}[H]
\begin{center}
\caption{Summary Statistics}
\begin{tabular}{l c c c c} \hline
& Treatment & Control & Diff & P-val on Diff \\ \hline
Number of households & 1260 & 2128 & & \\
Area of tobacco production (mu) & 5.578 & 4.874 & 0.705 & 0.000 \\
Age & 41.119 & 41.522 & -0.403 & 0.173 \\
Household size & 4.877 & 4.665 & 0.212 & 0.000 \\
Education (Primary) & 0.367 & 0.323 & 0.044 & 0.009 \\
Education (Secondary) & 0.602 & 0.338 & 0.263 & 0.000 \\
Education (High school or College) & 0.025 & 0.257 & -0.232 & 0.000 \\ \hline
\end{tabular}
\end{center}
\end{table}
We estimate the following linear-in-parameter model:
\begin{eqnarray}
Y_{i,2002} &=& Z_i'\beta_{2002}(U_{i,2002}), \nonumber \\
Y_{i,2003} &=& G_i \alpha(U_{i,2003}) + Z_i'\beta_{2003}(U_{i,2002}), \nonumber
\end{eqnarray}
where $Y_{it}$ is the tobacco production area (mu), $G_i$ is a treatment indicator equal to one for the treatment regions and zero for the control regions, and $Z_i$ is a vector of covariates including a constant term.
We estimate $\alpha(\tau)$ for $\tau = 0.1, ... , 0.9$.
Following \cite{cai2016impact}, we employ the age of the household head, household size, and education level indicators as control variables.
The main results from our method are presented in Figure 1.
The DID estimate is $0.239$, and the $95$ \% confidence interval is $[0.078,0.388]$.
We use the nonparametric bootstrap method to construct this confidence interval.
Figure 1 shows that the estimates of $\alpha(\tau)$ differ across $\tau$, and the QTE increases in $\tau$.
The impact of the insurance provision is nearly zero at the lower and middle quantiles and positive at the upper quantiles.
For $\tau = 0.8$ and $0.9$, the QTE is statistically significant.
In addition, we consider the null hypothesis that the QTEs are constant along $\tau$.
We conduct the test described in Remark 3 and calculate the test statistic and the critical value at the $0.05$ significance level.
These values are $571.9$ and $296.3$, respectively; hence, the null hypothesis is rejected.
\begin{figure}[h]
\centering
\includegraphics[width=15cm]{empirical_1.pdf}
\caption{The estimates of the QTEs and the 95 \% confidence intervals in Section 5.1. The horizontal axis measures the value of $\tau$ and the dashed line denotes the DID estimate.}
\end{figure}
\cite{cai2016impact} analyzes the welfare impact of the insurance program through the calibration.
The parameter values of the production function are chosen to match the DID (or triple difference) estimate.
From this analysis, she concludes that providing a heavily subsidized compulsory insurance program has a positive welfare impact on rural households.
However, our results show that the insurance program does not significantly change households' investment behavior at the lower and middle quantiles, and hence, may not affect household welfare at such quantiles.
\subsection{The TV effect on child cognitive development}
Next, we use our method to study the effect of TV on child cognitive development.
We use the data employed by \cite{huang2010dynamic} to estimate the QTE of TV watching on children's cognitive development.
This empirical analysis is based on a childhood longitudinal sample from NLSY79 (National Longitudinal Survey of Youth 1979).
Following \cite{huang2010dynamic}, we use a longitudinal sample of approximately 2,400 children and treat the Peabody Individual Achievement Test (PIAT) reading scores at ages 6--7 and 8--9 as $Y_{i1}$ and $Y_{i2}$, respectively.
The PIAT reading score at ages 6--7 has mean 103.0 and SD 11.7, and that at ages 8--9 has mean 104.3 and SD 14.6.
The outcome distribution at ages 8--9 is more dispersed than that at ages 6--7.
Hence, in these cases, additive time trends may not be plausible.
We use daily TV watching hours at ages 6--7 and 8--9 as the treatment variables $X_{i1}$ and $X_{i2}$, respectively and estimate the following linear-in-parameter model:
\begin{eqnarray}
Y_{it} &=& X_{it}'\alpha(U_{it}) + Z_{it}'\beta_t(U_{it}), \ \ \ \ t = 1, 2, \nonumber
\end{eqnarray}
where $Z_{it}$ is a vector of covariates including a constant term, dummy variables of race and gender, and an indicator of whether a child has 10 or more children's books at home.
In addition, we employ the Home Observation Measurement of the Environment variable (HOME), where is often used in child development research as an aggregate quality indicator of the home environment.
The main results from using our method are presented in Figure 2.
We find that the estimates of $\alpha(\tau)$ differ slightly across $\tau$ and the QTE is nearly zero at the lower and middle quantiles.
However, the impact of TV watching on the PIAT reading score is statistically significant at $\tau = 0.8$.
We consider the null hypothesis that the QTEs are zero at all quantiles.
We conduct the test described in Remark 3 and calculate the test statistic and the critical value at the $0.05$ significance level.
These values are $145.7$ and $290.5$, respectively; hence, the null hypothesis is not rejected.
Similar to \cite{huang2010dynamic}, the magnitude of the effect is quite small compared to the standard deviation of the PIAT reading score.
Therefore, the effect of TV on child cognitive development is neither statistically nor economically significant.
\begin{figure}[h]
\centering
\includegraphics[width=15cm]{empirical_2.pdf}
\caption{The estimates of the QTEs and the 95 \% confidence intervals in Section 5.2. The horizontal axis represents the value of $\tau$.}
\end{figure}
\section{Conclusion}
In this study, we developed a novel estimation method for the QTE under rank invariance and rank stationarity assumptions.
Although \cite{ishihara2020identification} also explores the identification and estimation of the nonseparable panel data model under these assumptions, the minimum distance estimation using this process is computationally demanding when the dimensionality of covariates is large.
To overcome this problem, we proposed a two-step estimation method based on the quantile regression and minimum distance methods.
We then showed the uniform asymptotic properties of our estimator and the validity of the nonparametric bootstrap.
The Monte Carlo studies indicated that our estimator performs well in finite samples.
Finally, we presented two empirical illustrations to estimate the distributional effects of insurance provision on household production, and TV watching on child cognitive development.
\clearpage
\setcounter{equation}{0}