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.
81,560 characters
Revisiting Panel Data Discrete Choice Models with Lagged Dependent Variables
\title{Revisiting Panel Data Discrete Choice Models with Lagged Dependent Variables\footnote{This paper uses unit record data from Household, Income and Labour Dynamics in Australia Survey (HILDA) conducted by the Australian Government Department of Social Services (DSS). The findings and views reported in this paper, however, are those of the author[s] and should not be attributed to the Australian Government, DSS, or any of DSS’ contractors or partners. DOI: 10.26193/YP7MNU.}\vspace{0.25cm}}
\author[1]{Christopher R. Dobronyi\footnote{E-mail address: \href{[email removed]}{[email removed]}.}}
\author[2]{Fu Ouyang\footnote{E-mail address: \href{[email removed]}{[email removed]}. }}
\author[3]{Thomas Tao Yang\footnote{E-mail address: \href{[email removed]}{[email removed]}.}}
\affil[1]{Google}
\affil[2]{School of Economics, University of Queensland}
\affil[3]{Research School of Economics, Australian National University}
\date{\today}
\maketitle
\begin{abstract}
\noindent This paper revisits the identification and estimation of a class of semiparametric (distribution-free) panel data binary choice models with lagged dependent variables, exogenous covariates, and entity fixed effects. We provide a novel identification strategy, using an ``identification at infinity'' argument. In contrast with the celebrated \cite{honore-k}, our method permits time trends of any form and does not suffer from the ``curse of dimensionality''. We propose an easily implementable conditional maximum score estimator. The asymptotic properties of the proposed estimator are fully characterized. A small-scale Monte Carlo study demonstrates that our approach performs satisfactorily in finite samples. We illustrate the usefulness of our method by presenting an empirical application to enrollment in private hospital insurance using the Household, Income and Labor Dynamics in Australia (HILDA) Survey data.
\end{abstract}
\noindent
{\it Keywords:} Dynamic Binary Choice Model, Fixed Effects, Identification at Infinity, Maximum Score Estimation
\vfill
\newpage
\section{Introduction}
In this paper, we propose new identification and estimation methods for
panel data binary choice models with fixed effects and ``dynamics'' (lagged dependent variables). Specifically,
suppose that there are $n$ individuals and $T+1$ time periods, $\{0,1,...,T\}
$. In each time period $t\in \{1,...,T\}$, each individual $i$ makes a
choice $y_{it}\in \{0,1\}$ according to the following latent utility model:
\begin{equation}
y_{it}=\mathds{1}\{\alpha _{i}+\gamma y_{it-1}+x_{it}^{\prime }\beta +\varpi
z_{it}\geq \epsilon _{it}\}, \label{model}
\end{equation}
where $\alpha _{i}$ is an entity fixed effect absorbing all relevant
time-invariant factors, $y_{it-1}$ is the lagged dependent variable, $
(x_{it},z_{it})$ is a $(p+1)$-vector of time-varying covariates, and $
\epsilon _{it}$ is an idiosyncratic error term. We separate $z_{it}$ from other covariates because, as it will be clear in the next section, we assign it a crucial role in the identification at infinity. In panel data literature, $
y_{it-1}$ is often called ``state
dependence'', and $\alpha _{i}$ is referred to as
``unobserved heterogeneity'' or
``spurious'' state dependence (see \cite
{Heckman1981b, Heckman1981a}). In model (\ref{model}), $
(y_{it},x_{it},z_{it})$ along with the ``initial
status'' $y_{i0}$ are observed in the data, whereas $\alpha
_{i}$ and $\epsilon _{it}$ are not observable to the econometrician. Note that we do not specify model (\ref{model}) in the initial period 0. This paper studies the identification and
estimation of the preference parameter $\theta :=(\gamma ,\beta ,\varpi )\in
\mathbb{R}^{p+2}$ in ``short'' panel
settings, i.e., $n\rightarrow \infty $ and $T<\infty $.
In line with the vast literature on panel data models with entity fixed effects, we do not impose any parametric restrictions on the distribution of $\alpha_i$ conditional on the initial choice $y_{i0}$ and observed covariates in model (\ref{model}). The prevalent methods for such models assume that $\epsilon_{it}$ are independently and identically distributed (i.i.d.) with a logistic distribution. \citet{ArellanoHonore2001}, \citet{honore2021identification}, and \citet{Hsiaobook} review various conditional likelihood approaches based on these parametric assumptions on $\epsilon _{it}$. Recent advances in the literature focus on constructing moment conditions for variants of dynamic Logit models. Representative works include \citet{honore-w}, \citet{dobronyi-gu-kim}, \citet{kitazawa2022transformations}, and \citet{dano2023transition}, among others.
\begin{comment}
In line with the vast literature on panel data models with entity fixed effects, we do not impose any parametric restrictions on the distribution of $\alpha_i$ conditional on the initial choice $y_{i0}$ and observed covariates in model (\ref{model}). Assuming $\epsilon_{it}$ are independently and identically distributed (i.i.d.) with a Logistic distribution, \cite{rasch} and \cite{andersen} show that, in the absence of $y_{it-1}$, this model is identifiable and can be estimated by the conditional maximum likelihood method. \cite{Chamberlain2010} proves that the Logit assumption is crucial for such models to have a non-zero information bound and achieve point identification without needing at least one observed covariate with unbounded support. \cite{chamberlain-85} and \cite{magnac} extend this approach to a Logit model with $y_{it-1}$ but no covariates, provided that there are at least four observations ($T\geq 3$) per individual. \cite{honore-k} advance this approach by developing a conditional maximum likelihood estimator for model (\ref{model}) that includes covariates, and \cite{hahn} examines the semiparametric efficiency of this estimator. Recent advances in variations of the Logit model can be found in \cite{johnson}, \cite{BartolucciNigro2010,BartolucciNigro2012}, \cite{bonhomme}, \cite{honore-w}, \cite{muris2020dynamic}, \cite{dobronyi-gu-kim}, \cite{honore2021dynamic}, \cite{kitazawa2022transformations}, and \cite{dano2023transition}, to name a few. Chapter 7 in \cite{Hsiaobook} provides a detailed review for readers unfamiliar with this body of literature.
\end{comment}
Without making distributional assumptions on $\epsilon _{it}$, \citet{manski-87} establishes the semiparametric identification of model (\ref{model}) that includes covariates, but not $y_{it-1}$. \cite{honore-k} extend this approach to include both $y_{it-1}$ and covariates in the model, showing that model (\ref{model}) with $T\geq 3$ can be identified under exogeneity and serial dependence assumptions stronger than those in \citet{manski-87}. However, their proposed estimator requires element-by-element matching of observed covariates over time, which rules out covariates with non-overlapping supports over time (e.g., time trend or dummies) and has a convergence rate decreasing in the dimension of the covariate space. \cite{OuyangYang2024binary} demonstrates that this curse of dimensionality can be mitigated by imposing certain serial dependence conditions on the covariates and by observing an extra time period. To highlight the novelty and contributions of this paper, we present a thorough comparison of our method with \citet{honore-k} and \cite{OuyangYang2024binary} in Appendices \ref{appendix0_1} and \ref{appendix0_2}, respectively.
There are alternative semiparametric and nonparametric approaches to model (\ref{model}). \citet{hl} demonstrate that model (\ref{model}) can be point identified if $z_{it}$ satisfies certain exclusion restrictions. More recently, \citet{ChenEtal2019} revisit this method and discuss the sufficient conditions for such exclusion restrictions. \citet{Williams2019} studies the nonparametric identification of dynamic binary choice models that satisfy certain exclusion restrictions. In the absence of excluded regressors, \citet{arist} establishes informative partial identification of model (\ref{model}) under weak conditions. \citet{KhanEtal2020} offer a partial identification result under even milder restrictions and prove that point identification is attainable in many interesting scenarios.
\begin{comment}
There are alternative semiparametric and nonparametric approaches to model (\ref{model}).\footnote{
Our literature review here focuses on fixed effects methods. It is well known that model (\ref{model}) can be estimated by the random effects or correlated random coefficients approach. See, e.g., \cite{ArellanoCarrasco2003}, \cite{Wooldridge2005}, and \cite{honore-tamer}, among others. These approaches often allow the econometrician to calculate choice probabilities and marginal effects in addition to $\theta$ at the cost of imposing restrictions on the distribution of $y_{i0}$ (initial condition problem) and the statistical relation between observed covariates and $\alpha_{i}$.} \citet{hl} demonstrate that the model, either with or without $y_{it-1}$, can be point identified if $z_{it}$ satisfies certain exclusion restrictions, i.e., the existence of an excluded regressor, regardless of the endogeneity of the covariates. More recently, \cite{ChenEtal2019} revisit this method and discusses the sufficient conditions for the exclusion restriction to be satisfied in dynamic settings. See also \cite{Williams2019} for nonparametric identification of dynamic binary choice models. In the absence of excluded regressors, \citet{arist} establishes informative partial identification of model (\ref{model}) under weak conditions. \cite{KhanEtal2020} offer a unique partial identification result under even milder restrictions, proving that point identification is attainable in many interesting scenarios.
\end{comment}
This paper revisits the distribution-free identification and estimation of model (\ref{model}). We show that the overlapping support restrictions required by \cite{honore-k} can be removed if $z_{it}$ is a free-varying covariate with full support. Our identification employs an ``identification at infinity'' strategy, first introduced in \cite{chamberlain1986asymptotic} and \cite{heckman-90}, and then applied in more recent work such as \cite{tamer2003incomplete}, \cite{bajari2010identification}, \cite{wan2014semiparametric}, and \cite{ouyang2020semiparametric}, among others. The combination of this strategy and \citeauthor{manski-87}'s (\citeyear{manski-87}) insight yields an estimator in the spirit of \citeauthor{honore-k}'s (\citeyear{honore-k}) conditional maximum score (MS) estimator, but without the need to match observed covariates over time. As a result, our estimator can accommodate flexible time effects and escape from the curse of dimensionality, in contrast to \cite{honore-k}. Through extending \cite{KimPollard1990} and \cite{SeoOtsu2018}, we demonstrate that our estimator converges at a rate slower than cube-root-$n$, is independent of the number of observed covariates, and has a non-standard limiting distribution. The asymptotics share similarities with those in \cite{honore-k} and \cite{OuyangYang2024binary}, with an important difference: the rate of convergence for our estimator depends on unknown factors, while the convergence rates of their estimators are known. We evaluate the finite-sample performance and implementability of our proposed estimator using both simulated and real-world data.
The rest of this paper is organized as follows. Section \ref{sec:identification} establishes the identification of $\theta$, which serves as the basis for the MS estimator presented in Section \ref{Sec:Estimation}. We then derive asymptotic properties of the proposed estimator in Section \ref{Sec:Asymptotics}. Results of Monte Carlo experiments are reported in Section \ref{sec_simulation}. We present an empirical illustration using the HILDA data in Section \ref{sec_application}. Finally, Section \ref{Sec:Conlcusions} concludes the paper with a brief discussion on possible future research directions. All proofs, supplementary discussions, and additional simulation results are included in the Supplementary Appendix.
For ease of reference, we list the notations maintained throughout this
paper here.
\begin{notation}
We reserve letter $i\in \{1,...,n\}$ for indexing individuals, and letter $
t\in \{1,...,T\}$ for indexing time periods. $\mathbb{R}^{k}$ is a $k$
-dimensional Euclidean space equipped with the Euclidean norm $\Vert \cdot
\Vert $ and $\mathbb{R}_{+}^{k}:=\{x\in \mathbb{R}^{k}|x>0\}$. We use $
P(\cdot )$ and $\mathbb{E}[\cdot ]$ to denote probability and expectation,
respectively. $\mathds{1}\{\cdot \}$ is an indicator function that equals
one when the event in the brackets occurs, and zero otherwise. Following a
substantial panel literature, we use the notation $\xi ^{t}$ to denote the
history of $\xi $ from period $1$ to period $t.$ For example, $x^{t}:=\left( x_{1},...,x_{t}\right) $ and $y^{t}:=\left(y_{1},...,y_{t}\right)$. We use $\setminus$ to denote set difference. For example, $\left( x_{1},x_{2},...,x_{t}\right)\setminus x_{1}:= \left(x_{2},...,x_{t}\right) $.
For two random vectors, $u$ and $v$, the notation $u\overset{d}{=}v|\cdot $
means that $u$ and $v$ have identical distribution, conditional on $\cdot $,
and $u\perp v|\cdot $ means that $u$ and $u$ are independent conditional on $
\cdot $. We use $\overset{p}{\rightarrow }$ and $\overset{d}{\rightarrow }$
to denote convergences in probability and in distribution, respectively. For
any (random) positive sequences $\{a_{n}\}$ and $\{b_{n}\}$, $a_{n}=O(b_{n})$
($O_{p}(b_{n})$) means that $a_{n}/b_{n}$ is bounded (in probability), $
a_{n}=o(b_{n})$ ($o_{p}(b_{n})$) means that $a_{n}/b_{n}\rightarrow 0$ ($
a_{n}/b_{n}\overset{p}{\rightarrow }0$). $a_{n}\lesssim b_{n}$ and $a_{n}\asymp b_{n}$ respectively mean
that there exist two constants $0<c_{1}\leq c_{2}<\infty$ such that
$c_{1}a_{n}\leq b_{n}$ and $c_{1}a_n\leq b_{n}\leq c_{2}a_{n}$. $a_{n}\ll b_{n}$ means $a_n = o(b_n)$.
\end{notation}
\section{Identification}\label{sec:identification}
Suppose in model (\ref{model}) $\epsilon_{it}$'s are i.i.d. over
time and independent of observed covariates $(x_{i}^{T},z_{i}^{T})$
and the initial choice $y_{i0}$ conditioning on the fixed effect
$\alpha_{i}$. The conditional probability of $y_{it}=1$ is equal
to:
\begin{equation}
P\left(y_{it}=1|\alpha_{i},y_{i}^{t-1},x_{i}^{T},z_{i}^{T}\right)=F_{\epsilon|\alpha}\left(\alpha_{i}+\gamma y_{it-1}+x_{it}^{\prime}\beta+\varpi z_{it}\right) \nonumber
\end{equation}
for each $i=1,\dots,n$ and $t=1,\dots,T$, where $F_{\epsilon|\alpha}(\cdot)$
denotes the cumulative distribution function (CDF) of $\epsilon_{it}$
conditional on $\alpha_{i}$. Consequently, the probability of observing
the choice history $y_{i}^{T}$ conditional on $(\alpha_{i},y_{i0},x_{i}^{T},z_{i}^{T})$
is expressed as:
\begin{align}
& P\left(y_{i}^{T}|\alpha_{i},y_{i0},x_{i}^{T},z_{i}^{T}\right)\nonumber\\
= & P\left(y_{i2},y_{i3},...,y_{iT}|\alpha_{i},y_{i0},y_{i1},x_{i}^{T},z_{i}^{T}\right)P\left(y_{i1}|\alpha_{i},y_{i0},x_{i}^{T},z_{i}^{T}\right)\nonumber \\
= & P\left(y_{i3},...,y_{iT}|\alpha_{i},y_{i0},y_{i1},y_{i2},x_{i}^{T},z_{i}^{T}\right)P\left(y_{i2}|\alpha_{i},y_{i0},y_{i1},x_{i}^{T},z_{i}^{T}\right)P\left(y_{i1}|\alpha_{i},y_{i0},x_{i}^{T},z_{i}^{T}\right)\nonumber \\
= & \cdots=\prod_{t=1}^{T}P\left(y_{it}|\alpha_{i},y_{i}^{t-1},x_{i}^{T},z_{i}^{T}\right)\nonumber \\
= & \prod_{t=1}^{T}F_{\epsilon|\alpha}\left(\alpha_{i}+\gamma y_{it-1}+x_{it}^{\prime}\beta+\varpi z_{it}\right)^{y_{it}}\left[1-F_{\epsilon|\alpha}(\alpha_{i}+\gamma y_{it-1}+x_{it}^{\prime}\beta+\varpi z_{it})\right]^{1-y_{it}} \nonumber
\end{align}
for each individual $i=1,\dots,n$.
In what follows, we will restrict the illustration of our approach
to model (\ref{model}) with $T=3$ and $\varpi>0$ to ease the exposition.
The condition $\varpi>0$ implies that we must know that the covariate
$z_{it}$ is included in the model and that it has a positive effect
on the choice probability of $y_{it}=1$. Applying our method to longer panels is straightforward, and the case with $\varpi<0$ is
symmetric. In addition, we will omit the subscript $i$ in our notation
whenever the context makes clear that all variables pertain to each
individual. Finally, we assume a balanced panel for simplicity. Our methods
are applicable to models with unbalanced panels, provided the
unbalancedness is not due to endogenous attrition.
Consider two choice histories
\[
C =\{y_{0}=d_{0},y_{1}=0,y_{2}=d_{2},y_{3}=1\} \textrm{ and } D =\{y_{0}=d_{0},y_{1}=1,y_{2}=d_{2},y_{3}=0\},
\]
where $d_{0},d_{2}\in\{0,1\}$. The conditional probability of the
choice history $C$ is equal to
\begin{align*}
& P(C|\alpha,y_{0}=d_{0},x^{T},z^{T})\\
= & p_{0}(\alpha,x^{T},z^{T})^{d_{0}}(1-p_{0}(\alpha,x^{T},z^{T}))^{1-d_{0}}(1-F_{\epsilon|\alpha}(\alpha+\gamma d_{0}+x_{1}^{\prime}\beta+\varpi z_{1}))\\
& \times F_{\epsilon|\alpha}(\alpha+x_{2}^{\prime}\beta+\varpi z_{2})^{d_{2}}(1-F_{\epsilon|\alpha}(\alpha+x_{2}^{\prime}\beta+\varpi z_{2}))^{1-d_{2}}F_{\epsilon|\alpha}(\alpha+\gamma d_{2}+x_{3}^{\prime}\beta+\varpi z_{3}),
\end{align*}
where $p_{0}(\alpha,x^{T},z^{T})$ denotes the conditional probability
of $y_{0}=1$. In a similar fashion,
\begin{align*}
& P(D|\alpha,y_{0}=d_{0},x^{T},z^{T})\\
= & p_{0}(\alpha,x^{T},z^{T})^{d_{0}}(1-p_{0}(\alpha,x^{T},z^{T}))^{1-d_{0}}F_{\epsilon|\alpha}(\alpha+\gamma d_{0}+x_{1}^{\prime}\beta+\varpi z_{1})\\
& \times F_{\epsilon|\alpha}(\alpha+\gamma+x_{2}^{\prime}\beta+\varpi z_{2})^{d_{2}}(1-F_{\epsilon|\alpha}(\alpha+\gamma+x_{2}^{\prime}\beta+\varpi z_{2}))^{1-d_{2}}(1-F_{\epsilon|\alpha}(\alpha+\gamma d_{2}+x_{3}^{\prime}\beta+\varpi z_{3})).
\end{align*}
Here, we take $d_{2}=1$ to illustrate, and the case with $d_{2}=0$ is
symmetric. Suppose the support of $z_{2}$ is unbounded above. Then, for $z_2>\sigma$, where $\sigma$ is a sufficiently large positive number, these probabilities satisfy
\begin{align}
& \frac{P(C|\alpha,y_{0}=d_{0},x^{T},z^{T})}{P(D|\alpha,y_{0}=d_{0},x^{T},z^{T})} \nonumber\\
\approx & \frac{1-F_{\epsilon|\alpha}(\alpha+\gamma d_{0}+x_{1}^{\prime}\beta+\varpi z_{1})}{1-F_{\epsilon|\alpha}(\alpha+\gamma d_{2}+x_{3}^{\prime}\beta+\varpi z_{3})} \times\frac{F_{\epsilon|\alpha}(\alpha+\gamma d_{2}+x_{3}^{\prime}\beta+\varpi z_{3})}{F_{\epsilon|\alpha}(\alpha+\gamma d_{0}+x_{1}^{\prime}\beta+\varpi z_{1})}. \label{eq:prob_ratio_1}
\end{align}
The idea of the above is to make negligible the effect of $y_{1}$
on $y_{2}$, i.e., $F_{\epsilon|\alpha}(\alpha+\gamma+x_{2}'\beta+\varpi z_{2})\approx F_{\epsilon|\alpha}(\alpha+x_{2}'\beta+\varpi z_{2})$
($\approx1$), by letting $z_2$ be sufficiently large.
Suppose $F_{\epsilon|\alpha}(\cdot)$ is strictly increasing. Then,
when $d_{2}=1$, equation (\ref{eq:prob_ratio_1}) implies that
\begin{align}
& \text{sgn}\left\{ P(C|\alpha,y_{0}=d_{0},x^{T},z^{T})-P(D|\alpha,y_{0}=d_{0},x^{T},z^{T})\right\} \nonumber \\
=& \text{sgn}\left\{ \gamma(d_{2}-d_{0})+(x_{3}-x_{1})^{\prime}\beta+\varpi(z_{3}-z_{1})\right\} \label{eq:iden_ineq1}
\end{align}
holds for $z_2>\sigma$ as $\sigma\rightarrow +\infty$, where $\text{sgn}\{\cdot\}$ is the
\textit{sign function}, which is equal to 1 if the expression inside
the brackets is strictly positive, to 0 if the expression inside the
brackets is zero, and to $-1$ if the expression inside the brackets
is strictly negative.
Equation (\ref{eq:iden_ineq1}) reveals that when $z_{2}$ is sufficiently large and $d_{2}=1$, the likelihood of observing event $C$ exceeds that of observing event $D$ if and only if $\gamma(d_{2}-d_{0})+(x_{3}-x_{1})^{\prime}\beta+\varpi(z_{3}-z_{1})>0$. In other words, the sign of $\gamma(d_{2}-d_{0})+(x_{3}-x_{1})^{\prime}\beta+\varpi(z_{3}-z_{1})$ determines the rank order of the conditional probabilities of events $C$ and $D$. Our focus on the subsample with $y_{3}\neq y_{1}$ aligns with \cite{manski-87} in forming his MS estimator. The distinction lies in our additional conditioning event of $z_{2}$ being large. It is worth noting that the same identification equation holds true when $-z_{2}$ is sufficiently large and $d_{2}=0$.
A natural way to build a population objective function based on equation
(\ref{eq:iden_ineq1}) is to define
\begin{align}
\bar{Q}_{1}(r,b,w):=\lim_{\sigma\rightarrow+\infty}\mathbb{E} & \left[\left(P(C|\alpha,y_{0}=d_{0},x^{T},z^{T})-P(D|\alpha,y_{0}=d_{0},x^{T},z^{T})\right)\right.\nonumber \\
& \left.\times\text{sgn}\left(r(d_{2}-d_{0})+(x_{3}-x_{1})^{\prime}b+w(z_{3}-z_{1})\right)|z_{2}>\sigma\right]\label{eq:pop_obj}
\end{align}
with $w>0$\ for $d_{2}=1$. By a symmetric argument, define
\begin{align}
\bar{Q}_{2}(r,b,w):=\lim_{\sigma\rightarrow+\infty}\mathbb{E} & \left[\left(P(C|\alpha,y_{0}=d_{0},x^{T},z^{T})-P(D|\alpha,y_{0}=d_{0},x^{T},z^{T})\right)\right.\nonumber \\
& \left.\times\text{sgn}\left(r(d_{2}-d_{0})+(x_{3}-x_{1})^{\prime}b+w(z_{3}-z_{1})\right)|z_{2}<-\sigma\right]\label{eq:pop_obj_2}
\end{align}
with $w>0$ for $d_{2}=0$. Note that equation (\ref{eq:iden_ineq1}) implies that $\bar{Q}_{1}(\gamma,\beta,\varpi)\geq\bar{Q}_{1}(r,b,w)$
and $\bar{Q}_{2}(\gamma,\beta,\varpi)\geq\bar{Q}_{2}(r,b,w)$ for
all $(r,b,w)\neq(\gamma,\beta,\varpi)$. Establishing that $\theta:=(\gamma,\beta,\varpi)$ is the unique maximum of either objective function (\ref{eq:pop_obj}) or (\ref{eq:pop_obj_2}) would confirm the point identification of these coefficients. The following conditions are sufficient for this.
\begin{assumptionp}{A} For all $\alpha$ and $s,t\in\mathcal{T}:=\{1,2,3\}$
($T=3$), the following conditions hold:\par \begin{enumerate}\par
\item[A1] (i) $\epsilon^{T}\perp(x^{T},z^{T},y_{0})|\alpha$,
(ii) $\epsilon_{s}\perp\epsilon_{t}|\alpha$, (iii) $\epsilon_{s}\overset{d}{=}\epsilon_{t}|\alpha$,
and (iv) conditional on $\alpha$, the CDF of $\epsilon_{t}$ is absolutely
continuous with support $\mathbb{R}$.\par
\item[A2] $z_{2}$ has unbounded support conditional on $(\alpha,y_{0},x^{T},z_{1},z_{3})$. \par
\item[A3] One of the elements in $\left(x_{3}-x_1,z_{3}-z_1\right)$,
denoted as $\xi_{31}$, has a bounded Lebesgue density that is positive
almost everywhere (a.e.) on $\mathbb{R}$ conditional on $(\alpha,x_{3}-x_1,z_{3}-z_1)\left\backslash \xi_{31}\right.$
and $\{z_{2}>\sigma\}\cup\{z_{2}<-\sigma\}$ as $\sigma\rightarrow +\infty$.
Moreover, the coefficient before $\xi_{31}$ is non-zero. \par
\item[A4] As $\sigma\rightarrow +\infty$, (i) the support $\mathcal{S}$
of $(y_{2}-y_{0},x_{3}-x_{1},z_{3}-z_{1})$ conditional on $\{z_{2}>\sigma\}$
or $\{z_{2}<-\sigma\}$ is not contained in any proper linear subspace
of $\mathbb{R}^{p+2}$, and (ii) the joint probability density function
(PDF) of $(x_{3}-x_{1},z_{3}-z_{1})$ conditional on $\{z_{2}>\sigma\}$
or $\{z_{2}<-\sigma\}$ is non-degenerate and uniformly bounded. \par
\item[A5] Let $\Theta$ be the set $\{\vartheta:=(r,b,w)\in\mathbb{R}^{p+2}\ |\ \Vert\vartheta\Vert=1,w>0\}$. $\theta$ is an interior
point of $\Theta$. \end{enumerate} \end{assumptionp}
Assumptions A1(i)--(iii) place the same restrictions on the joint
distribution of $(\alpha,\epsilon^{T},x^{T},z^{T})$ as \citet{honore-k},
which implies that the unobserved heterogeneity (entity fixed effects)
$\alpha$ picks up both the autocorrelation in the unobservables and
the dependence between explanatory variables and unobservables. As
a result, $\epsilon_{t}$ is independent of $(x^{T},z^{T},y^{t-1})$
conditional on $\alpha$ for all $t\in\mathcal{T}$. Assumption A1(iv)
is a regularity condition to guarantee that any possible sequence
of $y^{T}$ has a positive probability to occur.
Assumption A2 is a pivotal assumption that enables the ``identification at infinity'' approach, and when combined with Assumption A1, it establishes the identification equation
(\ref{eq:iden_ineq1}). It is clear from the derivation of equation (\ref{eq:iden_ineq1}) that relaxing this assumption may require additional restrictions on the parameter space $\Theta$,
the support of $x_{2}$, and the distribution of $(\epsilon,\alpha)$.
The support and continuity restrictions on
$\xi_{31}$ imposed by Assumption A3 are common for the family of
MS-type estimators, which are required to achieve the point
identification instead of a set identification. See, e.g., \citet{manski-75,manski-85,manski-87},
\citet{horowitz}, \citet{honore-k}, \citet{Fox2007}, \citet{ShiEtal2018},
\citet{YanYoo2019}, and \citet{khan2021inference}, among others.
Given the significance of Assumptions A2 and A3 in both our theoretical results and empirical application, we provide further discussion on them in Appendix \ref{appendix0_3}.
Assumption A4(i) is a familiar full-rank condition. Note that Assumptions
A3 and A4(ii) require $(x_{3}-x_{1},z_{3}-z_{1})$ to have sufficient
variation conditional on $\alpha$ and event $\{z_{2}>\sigma\}$ or
$\{z_{2}<-\sigma\}$. Assumption A5 applies the scale normalization
and restricts the search of $\theta$ in a compact set, which also
facilitates the asymptotic analysis of our estimator proposed in the
next section.
Additionally, we compare the key identification assumptions imposed in \citet{honore-k} and \citet{OuyangYang2024binary}, along with other aspects, with our method in Appendices \ref{appendix0_1} and \ref{appendix0_2}, respectively.
\begin{remark} \label{remark_A5}
Assumption A5 applies scale normalization by restricting $\vartheta$ to lie on a unit sphere. Alternatively, one can normalize one nonzero element of $\vartheta$, such as $w$ in this paper, to be 1. Following this convention, we express the parameter space as:\par\medskip
\noindent \textbf{Assumption A}5': $\Theta:=\{\vartheta:=(r,b,1)\in\mathbb{R}^{p+2}\}\cap\Xi$,
where $\Xi\subset\mathbb{R}^{p+2}$ is a compact set. $\theta$ is
an interior point of $\Theta$. \par\medskip
\noindent These two methods for scale normalization are essentially equivalent when $w \neq 0$. Therefore, researchers often use either of these methods based on their convenience in exposition or derivation. For instance, \cite{ShiEtal2018} use both methods in different sections. In this paper, we insist on Assumption A5 in Appendices \ref{appendixA} and \ref{appendixB}, as we derive the asymptotic properties of our estimator. This is because the two primary references for doing this, \cite{KimPollard1990} and \cite{SeoOtsu2018}, both normalize the parameter space to a unit sphere. Following the same convention facilitates our use of their established asymptotic theory and makes it easier for interested readers to review our proofs. When we apply our method to simulation studies and an empirical application in Sections \ref{sec_simulation} and \ref{sec_application}, we switch to the normalization defined in Assumption A5’, which reduces one parameter to estimate and avoids imposing restrictions on the optimization algorithm. We thank one anonymous referee for highlighting this point.
\end{remark}
Our identification results are stated in the following theorem and
the proof of which is provided in Appendix \ref{appendixA}.
\begin{theorem} \label{T:identify} Suppose Assumption A holds. Then,
$\theta$ is identified. \end{theorem}
\section{Estimation}\label{Sec:Estimation}
Applying the analogy principle, the population
objective functions (\ref{eq:pop_obj}) and (\ref{eq:pop_obj_2}) translate into MS estimation procedures (\cite{manski-75,manski-85,manski-87}).
Assume a random sample of $n$ observations is drawn from model (\ref{model}) that satisfies Assumption A.
Let $\vartheta :=(r,b,w)\in \mathbb{R}^{p+2}$ and $
\sigma _{n}\rightarrow \infty $ as $n\rightarrow \infty $. When the support
of $z_{i2}$ is unbounded above, we propose the MS estimator $\hat{\theta}_{n}
$ of $\theta $ maximizing the following objective function over the
parameter space $\Theta $:
\begin{equation}
Q_{n1}(\vartheta ):=\frac{1}{n}\sum_{i=1}^{n}y_{i2}(y_{i3}-y_{i1})\cdot
\mathds{1}\{z_{i2}>\sigma _{n}\}\cdot \mathds{1}
\{r(y_{i2}-y_{i0})+(x_{i3}-x_{i1})^{\prime }b+w(z_{i3}-z_{i1})>0\}.
\label{eq:Qn1}
\end{equation}
When the support of $z_{i2}$ is unbounded below, one can instead define $
\hat{\theta}_{n}$ with objective function
\begin{equation}
Q_{n2}(\vartheta ):=\frac{1}{n}\sum_{i=1}^{n}(1-y_{i2})(y_{i3}-y_{i1})\cdot
\mathds{1}\{z_{i2}<-\sigma _{n}\}\cdot \mathds{1}
\{r(y_{i2}-y_{i0})+(x_{i3}-x_{i1})^{\prime }b+w(z_{i3}-z_{i1})>0\}.
\label{eq:Qn2}
\end{equation}
If the support of $z_{i2}$ is unbounded both above and below, the objective
function can be a combination of (\ref{eq:Qn1}) and (\ref{eq:Qn2}) such as
\begin{equation}
Q_{n}(\vartheta ):=Q_{n1}(\vartheta )+Q_{n2}(\vartheta ). \label{eq:Qn}
\end{equation}
Note that (3.3) puts the same weight on $Q_{n1}(\vartheta)$ and $Q_{n2}(\vartheta)$,
which is a generic choice and probably not optimal in specific applications.
In some cases, it might be preferable to put more weight on one side
if additional information, such as restrictions on the error distribution,
suggests that the identification at infinity is more effective on
that side, especially if $z_{2}$ has a relatively heavier tail against
the error term.
Since $\text{sgn}(u)=2\cdot \mathds{1}\{u>0\}-1$ almost surely for any
continuous variable $u$, objective functions (\ref{eq:Qn1}) and (\ref{eq:Qn2}) are sample
analogues to monotone transformations of population functions (\ref{eq:pop_obj}) and (\ref
{eq:pop_obj_2}), respectively.
It is clear from expressions (\ref{eq:Qn1}) and (\ref{eq:Qn2}) that the
effective sample size for the estimator $\hat{\theta}_{n}$ is controlled by
the tuning parameter $\sigma_{n}$, and being similar to
\citeauthor{manski-87}'s (\citeyear{manski-87}) and \citeauthor{honore-k}'s (
\citeyear{honore-k}) estimators, only ``switchers'' who change
choices in periods 1 and 3 are used in the estimation. Besides, estimating (identifying) $\gamma$ relies on the variation in $y_{i2} - y_{i0}$, which means that we need some observations with $y_{i2}\neq y_{i0}$ and some with $y_{i2}=y_{i0}$.
Our proposed estimator $\hat{\theta}_n$ has two advantages, compared with
\citeauthor{honore-k}'s (\citeyear{honore-k}) estimator: First, the
estimation only needs to condition on a single univariate covariate, rather
than a vector of covariates, and hence it does not encounter the curse of
dimensionality. This property makes the procedure proposed above more
practical when the number of covariates is large. More importantly, our
estimator does not require matching $(x_t,z_t)$ in different periods.
Consequently, it allows covariates with non-overlapping support over time,
such as age, time trends, time dummy variables, etc.
\section{Asymptotic Properties}
\label{Sec:Asymptotics}
\subsection{Consistency}
This section establishes the asymptotic properties of the MS estimator
proposed in Section \ref{Sec:Estimation}. Given that objective functions (\ref
{eq:Qn1}) and (\ref{eq:Qn2}) are symmetric, it suffices to only investigate
the estimator $\hat{\theta}_{n}$ obtained from maximizing objective function
(\ref{eq:Qn1}) requiring the support of $z_{2}$ to be unbounded above. The
derivation for $\hat{\theta}_{n}$ associated with objective functions (\ref
{eq:Qn2}) or (\ref{eq:Qn}) is analogous. Additionally, for the sake of simplicity, we focus on the case where $\xi_{31}=z_{31}$ in Assumption A3.
To ensure the consistency of $\hat{\theta}_n$, we need additional technical
conditions.
\begin{assumptionp}{B}
For all $t\in \mathcal{T}$, the following conditions hold:
\begin{enumerate}
\item[B1] The data $\{y_{i0},y_{i}^{T},x_{i}^{T},z_{i}^{T}\}_{i=1}^{n}$ are
i.i.d. across $i$.
\item[B2] $\sigma_{n}$ is a sequence of positive numbers such that as $
n\rightarrow\infty$: (i) $\sigma_{n}\rightarrow\infty$, and (ii) $
nP(z_{2}>\sigma_{n})/\log n\rightarrow\infty$.
\item[B3] Let $\Lambda (\vartheta ):=y_{2}(y_{3}-y_{1})\cdot \mathds{1}
\{r(y_{2}-y_{0})+(x_{3}-x_{1})^{\prime }b+w(z_{3}-z_{1})>0\}$ for $\vartheta
\in \Theta $. Then (i) $\lim_{\sigma \rightarrow +\infty }\mathbb{E}[\Lambda
(\vartheta )|z_{2}>\sigma ]$ exists for all $\vartheta \in \Theta $, and
(ii) there exists an absolute constant $L$ such that
\begin{equation}
|\mathbb{E}[\Lambda (\vartheta _{1})|z_{2}>\sigma ]-\mathbb{E}[\Lambda
(\vartheta _{2})|z_{2}>\sigma ]|\leq L\Vert \vartheta _{1}-\vartheta
_{2}\Vert \label{assumption_lip}
\end{equation}
holds for all $\vartheta _{1},\vartheta _{2}\in \Theta $ and $\sigma >0$.
\end{enumerate}
\end{assumptionp}
Assumptions B2 imposes mild restrictions on the tuning parameter $\sigma _{n}
$. It is worth noting that Assumption B2(ii) indicates that the choice of $
\sigma_{n}$ depends on the tail behavior of the distribution of $z_{2}$.
For example, if $z_{2}$ has a sub-exponential right tail with $P(z_{2} > \sigma_{n}) \asymp e^{-\sigma_{n}}$, then any $\sigma_n$ satisfying $1 \ll \sigma_n \leq (1-\varepsilon) \log(n)$, e.g., $\sigma_n = \log\log(n/\log n)$, meets Assumption B2(ii), for some $\varepsilon\in(1/4,1)$. However, when the distribution of $z_{2}$ has a
(too) thin right tail $P(z_{2}>\sigma _{n})\asymp e^{-e^{\sigma _{n}}}$, $\sigma _{n}=\log\log(n/\log n)$ gives $nP(z_{2}>\sigma _{n})/\log n=O(1)$, violating
Assumption B2(ii). Notably, $nP(z_{2}>\sigma _{n})$ essentially controls the
``effective sample size'' for our proposed
procedure. As demonstrated in Theorem \ref{T:limiting_dist}, the tail behavior of the distribution of $z_{2}$ and the
choice of $\sigma _{n}$ jointly determine the convergence rate of the
proposed estimator $\hat{\theta}_{n}$. Assumption B3(ii) is a Lipschitz
condition essential for proving the uniform convergence of the objective
function (\ref{eq:Qn1}) to its population analogue. We provide a set of more concrete sufficient conditions for it in Appendix \ref{appendix0_34}.
The theorem below states that the proposed procedure described in (\ref
{eq:Qn1})--(\ref{eq:Qn})
gives a consistent estimator of $\theta$, whose
proof is left to Appendix \ref{appendixA}.
\begin{theorem}
\label{T:consistency} Suppose Assumptions A and B hold. Then, $\hat{\theta}_n
\overset{p}{\rightarrow}\theta$.
\end{theorem}
\subsection{Asymptotic Distribution}\label{Sec:asymp_dist}
We proceed to study the asymptotic distribution of
the estimator $\hat{\theta}_n$. Before presenting additional technical
conditions and the main results, we introduce some new notation to
facilitate exposition:
\begin{itemize}
\item[-] $h_{n}:=P(z_{2}>\sigma_{n}|y_{2}=1)$.
\item[-] For generic vectors $\xi_t$ and $\xi_s$, denote $
\xi_{ts}=\xi_t-\xi_s$.
\item[-] $\chi:=(y_{0},y^{T},x^{T},z^{T})$ and $
\bar{\chi}:=(y_{20},x_{31},z_{31})$.
\item[-] $u(\vartheta ):=\mathds{1}\{ry_{20}+x_{31}^{\prime }b+wz_{31}>0\}$,
thus, $u(\vartheta )=\mathds{1}\{\bar{\chi}^{\prime }\vartheta>0\} $.
\item[-] $\bar{q}_{n1,\vartheta }(\bar{\chi}):=\mathbb{E}\left[ y_{31}\left(
u(\vartheta )-u(\theta )\right) |z_{2}>\sigma _{n},y_{2}=1,\bar{\chi}\right] $
, and $\bar{q}_{1\vartheta }^{+}(\bar{\chi}):=\lim_{n\rightarrow \infty }\bar{q
}_{n1,\vartheta }(\bar{\chi})$.
\item[-] $\kappa _{n}(\bar{\chi}):=\mathbb{E}\left[ y_{31}|z_{2}>\sigma
_{n},y_{2}=1,\bar{\chi}\right] $ and $\kappa
^{+}(\bar{\chi}):=\lim_{n\rightarrow \infty }\mathbb{E}\left[
y_{31}|z_{2}>\sigma _{n},y_{2}=1,\bar{\chi}\right] $. $\dot{\kappa}_{n}(\nu
):=\left. \frac{\partial \kappa _{n}(\bar{\chi})}{\partial \bar{\chi}}
\right\vert _{\bar{\chi}=\nu }$ and $\dot{\kappa}^{+}(\nu ):=\left. \frac{
\partial \kappa ^{+}(\bar{\chi})}{\partial \bar{\chi}}\right\vert _{\bar{\chi}=\nu
}$.
\item[-] $F_{\bar{\chi}}(\cdot |z_{2}>\sigma _{n},y_{2}=1)$ ($
f_{\bar{\chi}}(\cdot |z_{2}>\sigma _{n},y_{2}=1)$) denotes the joint CDF (PDF)
of $\bar{\chi}$ conditional on $\{z_{2}>\sigma _{n},y_{2}=1\}$. $
F_{\bar{\chi}}^{+}(\cdot |y_{2}=1):=\lim_{n\rightarrow \infty
}F_{\bar{\chi}}(\cdot |z_{2}>\sigma _{n},y_{2}=1)$.
\end{itemize}
\begin{assumptionp}{C}
Suppose the following conditions hold. $\text{ }$
\begin{enumerate}
\item[C1] The proposed estimator $\hat{\theta}_{n}$ satisfies $Q_{n1}(\hat{
\theta}_{n})\geq\sup_{\vartheta\in\Theta}Q_{n1}(\vartheta)-o_{p}
((nh_{n})^{-2/3})$.
\item[C2] $P\left( z_{2}>\sigma |y_{2}=1,y_{31},\bar{\chi}\right) >0$ for all $
\sigma >0$ and almost every $(y_{31},\bar{\chi})$.
\item[C3] (i) $F_{\bar{\chi}}^{+}(\cdot |y_{2}=1)$ is non-degenerate and has
an uniformly bounded PDF $f_{\bar{\chi}}^{+}(\cdot |y_{2}=1)$, and (ii) $
\sup_{\nu }|f_{\bar{\chi}}(\nu |z_{2}>\sigma
_{n},y_{2}=1)-f_{\bar{\chi}}^{+}(\nu |y_{2}=1)|=o(1)$.
\item[C4] (i) $\kappa _{n}(\bar{\chi})$ and $\kappa ^{+}(\bar{\chi})$ are
differentiable in $\bar{\chi}$, and (ii) $\sup_{\nu }|\dot{
\kappa}_{n}(\nu )-\dot{\kappa}^{+}(\nu )|=o(1)$.
\item[C5] (i) $\mathbb{E}[\bar{q}_{n1,\vartheta }(\bar{\chi})]$ and $\mathbb{E
}[\bar{q}_{1\vartheta }^{+}(\bar{\chi})]$ are twice continuously
differentiable at $\vartheta $ in a small neighborhood of $\theta $, (ii) $
(\alpha ,x^{T})$ has a compact support, and (iii) for any constant $
\varsigma $, $\sup_{\alpha }P(\epsilon _{2}\geq \varsigma +\sigma
_{n}|\alpha )=o((nh_{n})^{-1/3})$.
\item[C6] $h_{n}\gtrsim n^{-1+\varepsilon }$ for some small positive $\varepsilon$.
\end{enumerate}
\end{assumptionp}
Assumption C1 is standard in the literature (see, e.g., \cite{KimPollard1990} and \cite{SeoOtsu2018}). This assumption implies that the maximization of $Q_{n1}(\vartheta)$ need not be exact, and any approximate maximizer close enough to the exact one will be enough for the asymptotic analysis.
Assumption C2 is an implication of Assumptions A1 and A2. We list it as a
separate condition here mainly because it is more directly related to our
proof process presented in Appendix \ref{appendixB}. Assumption C3
strengthens Assumption A4(ii). Assumption C4 requires the two conditional
probabilities $\kappa_{n}(\bar{\chi})$ and $\kappa ^{+}(\bar{\chi})$ to be
smooth enough, which, together with Assumption C3, is important for calculating the expected value of the limiting distribution of the
estimator $\hat{\theta}_{n}$.
The smoothness conditions imposed in Assumption C5(i) are standard in the literature as well.
Assumption C5(ii) is made to simplify the proof process and can be relaxed to allow for unbounded support, albeit with more tedious discussions. The essential requirement here is to exclude the scenario in which $\alpha + x_2^{\prime} \beta \rightarrow -\infty$ as $z_2 \rightarrow +\infty$.
Assumption C5(iii) essentially places a restriction on the relative tail
behavior of the observed regressor $z_t$ and unobserved error $\epsilon_t$. As shown in
the proof of Theorem \ref{T:limiting_dist}, this assumption ensures that the bias of the estimator $\hat{
\theta}_{n}$ shrinks sufficiently fast. If this condition is violated, the bias term dominates the distribution, and inferences are not possible. It is worth noting that such
condition plays a crucial role in determining the rate of convergence of
estimators based on ``irregular identification'' strategies including the
``identification at infinity'' as a special case. See \cite
{khan2010irregular} for an in-depth investigation on this issue.
Assumption C6, together with Assumption C5(iii), guides the selection of the tuning parameter $\sigma_{n}$. These two conditions are in the same spirit of Assumptions 8 and 8* in \cite{andrews}. On one hand, since $
h_{n}=P(z_{2}>\sigma_{n}|y_{2}=1)$ controls the effective sample size of the
estimation procedure, Assumption C6 implies that $\sigma_{n}$ should not increase too rapidly as $n\rightarrow\infty$, ensuring enough effective observations to control the variance of $\hat{\theta}_{n}$. On the
other hand, Assumption C5(iii) suggests that $\sigma_{n}$ should grow
sufficiently fast as $n\rightarrow\infty$ to lower the bias of $\hat{\theta}
_{n}$.
However, there is no way to determine the optimal $
\sigma_{n}$ since this requires the knowledge of relative (unknown) tail
behavior of $z_{t}$ and $\epsilon_{t}$. This feature is well known to the ``identification at infinity'' type of estimators, see, e.g., \cite{andrews}. We suggest choices of $\sigma_n$ that satisfy both Assumptions C5(iii) and C6 in some special cases in Table \ref{T:sigman}. From the table, there are no valid $\sigma_n$ in case (I) when $\lambda^{\prime}>\lambda$, and in case (III). The valid choices of $\sigma_n$, if exists, differ from case to case. As expected, we prefer the cases where $z_2$ possesses heavier tails than $\epsilon$. \cite{andrews} share similar results, for details, see their discussions after Assumption 8*.
For practice, we propose to take $\sigma_{n}=\sqrt{\log n/2.95}$. This choice of $\sigma_n$ is valid for case (I) with $\lambda^{\prime}=2$ and $\lambda'<\lambda$, and case (II) with $\lambda=2$. Moreover, with this choice of $\sigma_n$, Assumption B2(ii) is satisfied for $z_{2}$ with $P(z_{2} > \sigma_{n}) \gtrsim e^{-\left(1-\varepsilon\right) 2.95\sigma_{n}^{2}}$ for some $\varepsilon\in (1/4,1)$. We show the finite sample properties of our estimator with this choice of $\sigma_n$ by means of simulations in Section \ref{sec_simulation}. This chosen $\sigma_n$ works well (the bias does not dominate the distribution) even in the situation that belongs to case (I) with $\lambda^{\prime}>\lambda$, where no valid $\sigma_n$ exists.
\begin{table}[ptb]
\caption{Choice of $\sigma_n$ that satisfies both Assumptions C5(iii) and C6}
\small
\label{T:sigman}\centering
\begin{tabular}
[c]{l|c|c}\hline\hline
& $ P\left( \epsilon>t\right) \asymp \exp\left( -t^{\lambda}\right) $ & $P\left( \epsilon>t\right) \asymp t^{-\lambda}$\\\hline
$\begin{array}
[c]{c}
P\left( z_{2}>t\right) \\
\asymp \exp( -t^{\lambda^{\prime}})
\end{array} $ & (I):
$
\begin{array}{lc}
\sigma_{n}=(c\log n)^{1/\lambda'}\textrm{ }\forall c\in(0,1-\varepsilon], & \textrm{if }\lambda'<\lambda\\
\sigma_{n}=(c\log n)^{1/\lambda'}\textrm{ }\forall c\in(1/4,1-\varepsilon], & \textrm{if }\lambda'=\lambda\\
\textrm{No }\sigma_{n}, & \textrm{if }\lambda'>\lambda
\end{array}
$ & (III): $\text{No }\sigma_{n}$\\\hline
$\begin{array}
[c]{c}
P\left( z_{2}>t\right) \\
\asymp t^{-\lambda^{\prime}}
\end{array}$ & (II): $\left( \log n/3\right) ^{1/\lambda}
<\sigma_{n}\lesssim n^{\frac{1-\varepsilon}{\lambda^{\prime}}}$ & (IV):
$n^{\frac{1}{3\lambda+\lambda^{\prime}}}\ll\sigma_{n}\lesssim
n^{\frac{1-\varepsilon}{\lambda^{\prime}}}$\\
\hline\hline
\multicolumn{3}{l}{Note: We focus on the right tail, and $\lambda^{\prime}, \lambda > 0$.} \\
\end{tabular}
\end{table}
The above conditions are sufficient to characterize the asymptotic distribution of the estimator obtained by maximizing (\ref{eq:Qn1})--(\ref{eq:Qn}), as presented in the following theorem, along with its proof in Appendix \ref{appendixB}.
\begin{theorem}
\label{T:limiting_dist} Suppose Assumptions A--C hold. Then, (i) $\hat{\theta}
_{n}-\theta =O_{p}\left( (nh_{n})^{-1/3}\right) $, and (ii)
\begin{equation*}
\left( nh_{n}\right) ^{1/3}(\hat{\theta}_{n}-\theta )\overset{d}{\rightarrow
}\arg \max_{s\in \mathbb{R}^{p+2}}Z(s),
\end{equation*}
where $Z(s)$ is a Gaussian process with continuous sample paths, expected
value $s^{\prime }Vs/2$, and covariance kernel $H(s_{1},s_{2})$ for $
s_{1},s_{2}\in \mathbb{R}^{p+2}$. $V$ and $H(\cdot ,\cdot )$ are defined in
Lemma \ref{lemma_M1} and expression (\ref{eq:cov_kernel}), respectively.
\end{theorem}
Note that Theorem \ref{T:limiting_dist} does not determine the exact rate of convergence of $\hat{\theta}_n$, as $h_n$ depends on the unknown tail probabilities of $z_2$. However, the lack of this knowledge does not render statistical inference infeasible. In Section 6, we will apply the $m$-out-of-$n$ bootstrap to conduct the inference in an empirical application. We choose this method for two reasons: it is comparatively easier to implement, and it provides an estimate of the convergence rate for our estimator.
In Remark \ref{remark_boot}, we discuss several sampling-based methods with the potential to enable statistical inference in the absence of knowledge of the exact convergence rate of the estimator.
\begin{remark}
It is worth noting that the rate of convergence of $\hat{\theta}_{n}$ depends
on the relative tail behavior of the distributions of $z_{t}$ and
$\epsilon_{t}$. To achieve a faster convergence rate of $\hat{\theta}_{n}$, it
is desirable for the distribution of $z_{2}$ to have heavier tails compared to
$\epsilon_{2}$. To see this, consider any eligible $\sigma_{n}$ satisfying
both Assumptions C5 and C6. Suppose $P\left( \epsilon_{2}>\sigma_{n}\right)
\asymp h_{n}^{\upsilon}$ (the bias term) for some $\upsilon>0$, and
$h_{n}^{\upsilon}\asymp(nh_{n})^{-1/3}$ for the fastest possible convergence
rate $n^{-\upsilon/\left( 1+3\upsilon\right) }$. Recall that $P\left(
z_{2}>\sigma_{n}\right) =h_{n}.$ When $z_{2}$ and $\epsilon
_{2}$ have the same tail, $\upsilon=1$ and the convergence rate is $n^{-1/4}$.$\ $Loosely
speaking$,$ $\upsilon\ $increases as the tail of $z_{2}$ becomes thicker, and
decreases otherwise. Therefore, for any eligible $\sigma_{n}$, the convergence rate
of $\hat{\theta}_{n}$ increases in $\upsilon$ (as the tail of $z_{2}$ becomes
thicker) and approaches $n^{-1/3}$ for large $\upsilon$ (as the tail of $z_{2}$ becomes much
thicker than $ \epsilon_{2}$).
\end{remark}
\begin{remark}\label{remark_boot}
Theorem \ref{T:limiting_dist} indicates that the proposed estimator $\hat{\theta}_n$ has a slower than cube-root-$n$ rate of convergence and its
asymptotic distribution is not Gaussian. As a result, standard inference methods
based on asymptotic normality no longer work here. Smoothing the objective
function in the sense of \cite{andrews} and \cite{horowitz} (See also \cite{kyriazidou1997estimation} and \cite{charlier1997limited}) may yield a
faster rate and regain an asymptotically normal estimator. However, this
involves choosing additional kernel functions and tuning parameters for the
two indicator functions in objective functions (\ref{eq:Qn1}) and (\ref{eq:Qn2}). A more practical alternative may be to consider sampling-based
inference methods. It is known that the naive nonparametric bootstrap is
typically invalid under the cube-root asymptotics (\cite{AbrevayaHuang2005}). For the ordinary MS estimator, valid inference can be conducted using
subsampling (\cite{DelgadoEtal2001}), the $m$-out-of-$n$ bootstrap (\cite{LeePun2006}), the numerical bootstrap (\cite{HongLi2020}), and a model-based bootstrap with modified objective function (\cite{cattaneo2020bootstrap}), among other procedures. \cite{OuyangYang2024binary} show that \citeauthor{LeePun2006}'s (\citeyear{LeePun2006}), \citeauthor{HongLi2020}'s (\citeyear{HongLi2020}),
and \citeauthor{cattaneo2020bootstrap}'s (\citeyear{cattaneo2020bootstrap})
methods, with certain modifications, are valid for kernel weighted MS
estimators with asymptotics similar to Theorem \ref{T:limiting_dist}.
Similar methods might apply to the estimator proposed in this paper. However, extending these bootstrap methods to cases with unknown convergence rates requires significant effort and is beyond
the scope of this paper. We, therefore, defer this task to future studies.
\end{remark}
\section{Monte Carlo Experiments}\label{sec_simulation}
In this section, we investigate the finite-sample performance of the
proposed estimators by means of Monte Carlo experiments. We examine
two designs, each with a less favorable scenario for our estimator.
In these scenarios, $z_{t}$ has a thinner tail than $\epsilon_{t}$,
and there is no theoretically valid $\sigma_{n}$. These are the first
scenarios in both Designs 1 and 2 presented below. Despite these challenges,
our estimator performs reasonably well, as the bias term does not
appear to dominate the distribution.
We start by considering a benchmark design similar to that used in
\citet{honore-k}, but we add an additional covariate $z_{it}$ and
a time trend that \citet{honore-k} cannot handle. Specifically, this
design (referred to as Design 1) is specified as follows:
\begin{align*}
y_{i0} & =\mathds{1}\left\{ \alpha_{i}+\delta\times\left(0-2\right)+\beta_{1}x_{i0,1}+z_{i0}\geq\epsilon_{i0}\right\} ,\\
y_{it} & =\mathds{1}\left\{ \alpha_{i}+\delta\times\left(t-2\right)+\gamma y_{it-1}+\beta_{1}x_{it,1}+z_{it}\geq\epsilon_{it}\right\} ,\text{ \ }t\in\left\{ 1,2,3\right\} ,
\end{align*}
where we set $\gamma=\beta_{1}=1$ and $\delta=1/2$. Following the discussion in Remark \ref{remark_A5}, we normalize the coefficient on $z_{it}$ to 1 for all designs investigated in this section and Appendix \ref{appendixC}. We consider
two scenarios. For each, we let $x_{it,1}\overset{d}{\sim}N\left(0,1\right),$
$\epsilon_{it}\overset{d}{\sim}(\pi^{2}/3)^{-1/2}\cdot$Logistic$\left(0,1\right)$
(the variance of $\epsilon_{it}$ is 1)$,$\ and $\alpha_{i}=\left(x_{i0,1}+x_{i1,1}+x_{i2,1}+x_{i3,1}\right)/4,$\ but
we consider $z_{it}$ with different tail behaviors. $x_{\cdot,1},z_{\cdot},$
and $\epsilon_{\cdot}$ are independent of each other, and all covariates
are i.i.d. across $i$ and $t.$ In the first scenario, we set $z_{it}\overset{d}{\sim}N(0,1),$
and denote it as ``Norm''. In the second scenario, we set $z_{it}\overset{d}{\sim}\text{Laplace}(0,\sqrt{2}/2)$
(with zero mean and unit variance), and denote it as ``Lap''. Note
that the density function of the Laplace distribution decays like
$e^{-\left\vert x\right\vert /c}$ for some constant $c$ at its tail$,$
which is heavier than the tail of the normal density.
In the second design (referred to as Design 2), the setup is the same
as that in Design 1, except that we add one more covariate to examine how
our estimators perform in a higher dimensional design. Specifically,
\begin{align*}
y_{i0} & =\mathds{1}\left\{ \alpha_{i}+\delta\times\left(0-2\right)+\beta_{1}x_{i0,1}+\beta_{2}x_{i0,2}+z_{i0}\geq\epsilon_{i0}\right\} ,\\
y_{it} & =\mathds{1}\left\{ \alpha_{i}+\delta\times\left(t-2\right)+\gamma y_{it-1}+\beta_{1}x_{it,1}+\beta_{2}x_{it,2}+z_{it}\geq\epsilon_{it}\right\} ,\text{ \ }t\in\left\{ 1,2,3\right\} ,
\end{align*}
where we set $\gamma=\beta_{1}=\beta_{2}=1$ and $\delta=1/2$. Random
covariates are generated as$\ x_{it,1},x_{it,2}\overset{d}{\sim}N\left(0,\sqrt{2}/2\right)$,
$\epsilon_{it}\overset{d}{\sim}(\pi^{2}/3)^{-1/2}\cdot\text{Logistic}\left(0,1\right)$,
and $\alpha_{i}=\sum_{t=0}^{3}(x_{it,1}+x_{it,2})/4$. Similarly,
we consider two scenarios with the same $z_{it}$ as in design 1$.$
Again, $x_{\cdot,1}$, $x_{\cdot,2}$, $z_{\cdot}$, and $\epsilon_{\cdot}$
are independent of each other. To investigate only the impact of higher dimension, we
set the variance of $x_{it,1}+x_{it,2}$ in Design 2 to be the same
as that of $x_{it,1}$ in Design 1.
As discussed in Section \ref{Sec:asymp_dist}, we set
$\sigma_{n}$ as
\[
\sigma_{n}=\widehat{\text{std}\left(z_{i2}\right)}\sqrt{\log n^{*}/2.95},
\]
where $\widehat{\text{std}\left(z_{i2}\right)}$ is the sample standard
deviation of $z_{2}$, and $n^{*}$ is the number of ``switchers'',
that is, observations with $y_{3}\neq y_{1}$. The usage of $n^{*}$
is intended to provide better control over the tuning parameters,
based on the features of the data. In practice, one may normalize
$z_{it}$ to mean 0 and variance 1 and set $\sigma_{n}=\sqrt{\log n^{*}/2.95}$.
We consider sample sizes of $n=5000,10000$, and $20000$. All the simulation results presented in this section are based on
1000 replications of each sample size. We implement MS
estimations in R, using the differential evolution (DE)
algorithm to attain a global optimum of the
objective function. The DE algorithm, developed by \citet{storn1997differential},
is capable of searching for the global optimum of a real-valued function
with real-valued parameters, even if the function lacks continuity
or differentiability. This algorithm has been effectively employed
in calculating MS-type estimators in the literature, including \citet{Fox2007}
and \citet{YanYoo2019}. \citet{mullen2011deoptim} provides a comprehensive
introduction to the R package $\texttt{DEoptim}$, which implements
the DE algorithm. We report the mean bias (MBIAS) and the root mean square errors (RMSE) of the estimates for Designs
1 and 2 in Tables \ref{T:D1} and \ref{T:D2}, respectively.
We summarize the findings in Tables \ref{T:D1} and \ref{T:D2}. First,
the RMSEs of all parameters decrease as the sample size increases,
but they converge to zero slower than the parametric rate. Second,
the convergence rate is faster with a thicker-tailed $z_{i2}$, as evidenced by comparing the RMSEs from Norm to Lap. Third, the RMSE does not appear to increase
for $\gamma$ and $\delta$ as we have one more covariate from Design
1 to Design 2. This confirms our theoretical findings. Note that the RMSE increases a bit for $\beta_{1}$, but this is probably due to
the lower variance of $x_{\cdot,1}$ in Design 2. To investigate the
sensitivity of the results to $\sigma_{n},$ we consider $\sigma_{n}=0.9\cdot\widehat{\text{std}\left(z_{i2}\right)}\sqrt{\log n^{*}/2.95}$
and $\sigma_{n}=1.1\cdot\widehat{\text{std}\left(z_{i2}\right)}\sqrt{\log n^{*}/2.95}$
(we need larger $\sigma_{n}$ to be in line with the discussion in
Section \ref{Sec:asymp_dist}), and report the corresponding results
in Tables \ref{T:D1_robust} and \ref{T:D2_robust} in Appendix \ref{appendixC}.
We note that the results are not sensitive to the choices of the tuning
parameters.
\begin{comment}
In Appendix \ref{appendixC}, we examine the impact of auto-correlations
of regressors on the performance of our estimators. After removing
the time trend term, we construct a new design and we compare the
performance of our estimator with the estimators in \citet{honore-k}
and \citet{OuyangYang2024binary}.\footnote{These two competing methods are not applicable in the presence of time
trends and dummies.} We briefly summarize the results in Appendix \ref{appendixC}. With
certain degrees of autocorrelations, our estimator performs reasonably
well, yet does not perform as well as before. Our estimator's performance
is comparable to that of the semiparametric estimators in \citet{honore-k}
and \citet{OuyangYang2024binary}. Notably, these competing methods are not valid in scenarios involving time trends and dummies, which are prevalent in empirical applications. In such contexts, our approach offers a valuable alternative.
\end{comment}
In Appendix \ref{appendixC}, we report additional results from supplementary simulation studies. Firstly, we investigate the impact of auto-correlations of the regressors on the performance of our estimator. Additionally, we compare the performance of our estimator with those proposed by \citet{honore-k} and \citet{OuyangYang2024binary} in designs without the time trend term. We direct interested readers to Appendix \ref{appendixC} for a more detailed discussion. Here, we provide a brief summary of these results. Our estimator still performs reasonably well with certain degrees of auto-correlations, but as expected, not as well as in Designs 1 and 2, where regressors are serially independent. Our estimator's performance is comparable to that of the semiparametric estimators proposed by \citet{honore-k} and \citet{OuyangYang2024binary}. It is essential to highlight that these alternative methods are not applicable in scenarios involving time trends or dummies, which are common in empirical applications. In such contexts, our approach offers a valuable alternative.
\begin{table}[H]
\caption{Simulation Results of Design 1}
\label{T:D1}\centering
\begin{tabular}{r|cc|cc|cc}
\hline \hline
& \multicolumn{2}{c|}{$\beta_{1}$} & \multicolumn{2}{c|}{$\gamma$} & \multicolumn{2}{c}{$\delta$}\tabularnewline
& MBIAS & RMSE & MBIAS & RMSE & MBIAS & RMSE \tabularnewline
\hline
$n_1$ & 0.120 & 0.407 & -0.027 & 0.549 & 0.068 & 0.228 \\
Norm $n_2$ & 0.075 & 0.299 & -0.015 & 0.427 & 0.053 & 0.171 \\
$n_3$ & 0.048 & 0.219 & -0.062 & 0.340 & 0.043 & 0.132 \\ \hline
$n_1$ & 0.039 & 0.249 & -0.021 & 0.386 & 0.030 & 0.154 \\
Lap $n_2$ & 0.024 & 0.185 & -0.049 & 0.306 & 0.019 & 0.116 \\
$n_3$ & 0.027 & 0.146 & -0.041 & 0.247 & 0.020 & 0.095 \\
\hline
\hline
\multicolumn{7}{l}{Note: $n_{1}=5000,n_{2}=10000,n_{3}=20000$.}\tabularnewline
\end{tabular}
\end{table}
\begin{table}[H]
\caption{Simulation Results of Design 2}
\label{T:D2}\centering
\begin{tabular}{r|cc|cc|cc|cc}
\hline \hline
& \multicolumn{2}{c|}{$\beta_{1}$} & \multicolumn{2}{c|}{$\beta_{2}$} & \multicolumn{2}{c|}{$\gamma$} & \multicolumn{2}{c}{$\delta$}\tabularnewline
& MBIAS & RMSE & MBIAS & RMSE & MBIAS & RMSE & MBIAS & RMSE \tabularnewline
\hline
$n_1$ & 0.144 & 0.471 & 0.157 & 0.475 & 0.018 & 0.542 & 0.093 & 0.235 \\
Norm $n_2$ & 0.087 & 0.344 & 0.085 & 0.355 & -0.032 & 0.430 & 0.057 & 0.172 \\
$n_3$ & 0.048 & 0.276 & 0.062 & 0.273 & -0.048 & 0.332 & 0.044 & 0.137 \\\hline
$n_1$ & 0.072 & 0.314 & 0.078 & 0.313 & -0.009 & 0.395 & 0.046 & 0.154 \\
Lap $n_2$ & 0.023 & 0.216 & 0.036 & 0.238 & -0.025 & 0.303 & 0.023 & 0.111 \\
$n_3$ & 0.025 & 0.182 & 0.022 & 0.178 & -0.035 & 0.236 & 0.021 & 0.091 \\
\hline
\hline
\multicolumn{9}{l}{Note: $n_{1}=5000;n_{2}=10000;n_{3}=20000$.}\tabularnewline
\end{tabular}
\end{table}
A final note is that when using observational data, the choice of $\sigma_n$ depends on the unknown tail behavior of the variable $z_2$. As there are no formal methods to determine the appropriateness of a specific $\sigma_n$, we suggest practitioners try different $\sigma_n$'s in estimation and check if the results are sensitive to different choices.
\section{Empirical Illustration} \label{sec_application}
In Australia, Medicare is the universal tax-funded public health insurance scheme that provides free access to public
hospitals. Medicare patients in public hospitals receive free treatment from
doctors nominated by hospitals and free (shared) accommodations. Patients
may opt to receive private care in either private or public hospitals (as
private patients) to have their choice of doctors and nurses, better
amenities (e.g., private rooms, family member accommodation, etc.), and
quicker access to treatment by avoiding long waiting time experienced by
many Medicare patients. Medicare does not cover private hospital care. On
top of a patient copayment, the cost is either afforded by private patients
themselves as out-of-pocket expenditure or covered by their private hospital
(insurance) cover (PHC), if any. Having PHC does not preclude using hospital
care as a Medicare patient. The institutional context for Australia's
Medicare and private health insurance schemes has been more thoroughly
described in the vast health economics literature, e.g., Section 2 of \cite
{cheng2014measuring}. We refer interested readers to \cite
{cheng2014measuring} and references therein for more detailed information.
In this section, we apply our MS estimator to analyze the state dependence
and the impacts of government incentives on the choice to purchase PHC,
using 10 waves (waves 11--20 corresponding to years 2011--2020) of the
Household, Income and Labor Dynamics in Australia (\href{https://melbourneinstitute.unimelb.edu.au/hilda}{HILDA}) Survey data. Since 2011, the HILDA survey has begun recording information about respondents' enrollment in PHC.
We denote the dependent variable, $y_{it}$, as whether individual $i$ has PHC in year $t$. We are interested in the effects of ``Lifetime Health
Cover'' (LHC) policy, ``Medicare Levy
Surcharge'' (MLS), and the state persistence $(y_{it-1})$
on one's purchasing PHC.
The age dummy variable $\text{Above30}_{it}$ indicates if individual $i$ is
30 years old or above in year $t$, namely, $\text{Above30}_{it}:=\mathds{1}
\{\text{Age}_{it}\geq30\}$. Following the insight of the (sharp)
``regression discontinuity'' design, its coefficient captures the
effects of Australia's LHC policy introduced in 2000 to encourage the uptake
of PHC. Loosely speaking, the LHC states that if an individual has not taken
out and maintained PHC from the year she turns 31, she will pay a 2\% LHC
loading on top of her premium for every year she is aged over 30 if she
decides to take out PHC later in life. If LHC is a strong incentive, we
would expect a significant ``jump" in the PHC enrollment rate at this age.
The MLS is a levy paid by Australian taxpayers who do not have PHC and earn
above a certain income threshold. In the sample years of our data, MLS
rates remain unchanged, while the thresholds have been raised yearly until
2014. It is worth noting that the 2014 rise in MLS thresholds was only 50\%
of previous years, and the thresholds have remained at the same level until
2022. The time dummy $D_{2014,t}$ is included in model (\ref{app_model}
) to examine whether this change in MLS policy would affect people's
willingness to purchase PHC. Note that \citeauthor{honore-k}'s (
\citeyear{honore-k}) estimators do not allow either age ($\text{Age}_{it}$)
or fixed time effects ($D_{2014,t}$) since they do not have overlapping
supports across time.
$I_{it}$ represents standardized annual household disposable income
using the entire sample in the survey, which serves as the continuous regressor
with rich enough support required for point identification (by Assumptions A2 and A3). The standardization is performed before dropping missing data by subtracting the sample mean from each individual value and then dividing the difference by the standard deviation.
We also include a (location) dummy variable $\text{GCC}_{it}$ that indicates whether individual $i$ lives in a major city/greater capital city in year $t$. This variable is included to control the accessibility to private hospital services. Tables \ref{tab_def} and \ref{tab_sum} provide definitions and summary statistics of all aforementioned variables, respectively. Note that some observations are excluded due to missing information in other variables, so in Table \ref{tab_sum}, $I_{it}$ does not have an exact zero mean and unit standard deviation.
With all these covariates, we specify our empirical model as follows:
\begin{equation}
y_{it}=\mathds{1}\{\alpha _{i}+\gamma y_{it-1}+\delta D_{2014,t}+\beta _{1}
\text{Abov30}_{it}+\beta _{2}\text{Age}_{it}+\beta _{3}\text{GCC}
_{it}+ I_{it}\geq \epsilon _{it}\}, \label{app_model}
\end{equation}
where $\epsilon _{it}$ and $\alpha _{i}$ are, respectively, the usual idiosyncratic error and unobserved heterogeneity in fixed effects panel data models.
In our analysis, we restrict the coefficient on $I_{it}$ to be 1, following the same convention for scale normalization as in Section \ref{sec_simulation}. This choice warrants justification; that is, household income enters the model with a significant positive coefficient, as implicitly required by Assumption A5'. We provide the following rationale for this based on common sense and evidence from exploratory regression. Practitioners seeking to justify normalizing the coefficient on $z_{it}$ to 1 can adopt similar argumentation method.
Firstly, in Australia, Medicare provides free access to public hospitals, and Medicare patients in public hospitals receive free treatment and accommodations. However, people can purchase private hospital insurance to cover faster and more premium services. Taking up or maintaining private insurance coverage requires a household to have sufficient disposable income. Besides, as income increases, the marginal utility of saving or other consumption may eventually become lower than that of enhanced private health care. In addition, Australia's tax system also gives considerable financial incentives for high-income households to buy private insurance. Therefore, common sense suggests that income should play a positive and significant role in private insurance purchases.
Secondly, we conduct a simple probit regression using one wave of the data and included income as the only regressor. The estimate is positive and significant, with $p$-value smaller than $10^{-15}$. This result holds true across all data waves, confirming our argument. This finding aligns with the results of more in-depth structural analyses in the health economics literature, such as \cite{cheng2014measuring}.
Note that Assumption A3 can be demanding. To address this concern, we relax Assumption A3 to Assumption A3' for a model closely resembling the current application and demonstrate identification under this relaxed condition, as detailed in Appendix \ref{appendix0_32}. Additionally, we justify our use of \( I_{it} \) as \( z_{it} \) under this modified condition in Appendix \ref{appendix0_33}, specifically by showing the kernel density and summary statistics of \( I_{it+1} - I_{it-1} \). For a more detailed discussion on Assumptions A2 and A3 and their roles in this empirical application, we refer interested readers to Appendix \ref{appendix0_3}.
\begin{table}[tpb]
\caption{Definition of Variables}
\label{tab_def}{\small \centering
\begin{tabularx}{1.0\textwidth}[c]{l|X}
\hline\hline
Variable & Description\\\hline
Private hospital cover ($y$) & 1 if has private hospital cover for the whole year, otherwise 0 \\\hline
Standardized income ($I$) & Standardized household's financial year disposable income (in the 2011 Australian dollar) \\\hline
Above 30 years old (Above30) & 1 if age 30 years old or above, otherwise 0 \\\hline
Age & Age \\\hline
Major city or greater capital city (GCC) & 1 if lives in a major city or greater capital city, otherwise 0 \\\hline
Year 2014 ($D_{2014}$) & 1 if in financial year 2014, otherwise 0 \\\hline\hline
\end{tabularx}
}
\end{table}
\begin{table}[tpb]
\caption{Summary Statistics}
\label{tab_sum}
\centering
\begin{tabular}{lccccc}
\hline\hline
Variable & $n\times T$ & Mean & Std.Dev. & Min & Max \\
\midrule $y_{it}$ & 65,603 & 0.527 & 0.499 & 0 & 1 \\
$I_{it}$ & 65,603 & -0.128 & 0.946 & -1.401 & 13.118 \\
$\text{Above30}_{it}$ & 65,603 & 0.855 & 0.353 & 0 & 1 \\
$\text{Age}_{it}$ & 65,603 & 50.091 & 17.484 & 17 & 99 \\
$\text{GCC}_{it}$ & 65,603 & 0.588 & 0.492 & 0 & 1 \\
$D_{2014,t}$ & 65,603 & 0.136 & 0.343 & 0 & 1 \\ \hline\hline
\end{tabular}
\end{table}
Let $x_{it}:=(\text{Above30}_{it},\text{Age}_{it},\text{GCC}_{it})$ and $
\beta :=(\beta _{1},\beta _{2},\beta _{3})$. We estimate $\theta :=(\delta
,\gamma ,\beta)$ through maximizing the objective function
\begin{equation}
Q_{n1}(\vartheta ):=\frac{1}{n}\sum_{i=1}^{n}
\sum_{t=2}^{T_{i}-1}y_{it}(y_{it+1}-y_{it-1})\cdot \mathds{1}\{I_{it}>\sigma
_{n}\}\cdot \mathds{1}\{u_{it}(\vartheta )>0\}, \label{app_Qn1}
\end{equation}
where $u_{it}(\vartheta
):=r(y_{it}-y_{it-2})+d(D_{2014,t+1}-D_{2014,t-1})+(x_{it+1}-x_{it-1})^{
\prime }b+(I_{it+1}-I_{it-1})$ and $\vartheta :=\left(d, r,b\right)$. Objective function (\ref{app_Qn1}) extends (\ref{eq:Qn1}) for longer and unbalanced panels in which the number of waves being observed varies across
individuals $i$ ($=:T_{i}$). We select the tuning parameter $\sigma_{n}$ using the same approach as described in Section \ref{sec_simulation}. It is important to note that the distribution of $I_{it}$ exhibits a significantly longer right tail compared to its left tail (skewness=3.78). Consequently, for sufficiently large $\sigma_{n}$, the objective function (\ref{app_Qn1}) has a much larger number of observations to use than the objective function
\begin{equation}
Q_{n2}(\vartheta ):=\frac{1}{n}\sum_{i=1}^{n}
\sum_{t=2}^{T_{i}-1}(1-y_{it})(y_{it+1}-y_{it-1})\cdot \mathds{1}
\{I_{it}<-\sigma _{n}\}\cdot \mathds{1}\{u_{it}(\vartheta )>0\}, \notag \label{app_Qn2}
\end{equation}
which extends (\ref{eq:Qn2}) for left tail observations. In fact, in this application, we set $\sigma_{n}=1.478$, which exceeds the absolute value of the lower bound of $I_{it}$ ($=1.401$ as shown in Table \ref{tab_sum}), thereby effectively using only objective function (\ref{app_Qn1}) and observations satisfying $\{I_{it}>\sigma_{n}\}$. Previous versions of this paper explored smaller values of $\sigma_n$ that allowed for the inclusion of left-tail observations (i.e., $\{I_{it}<-\sigma_{n}\}$), yielding similar results.
\begin{comment}
Let $x_{it}:=(\text{Above30}_{it},\text{Age}_{it},\text{GCC}_{it})$ and $
\beta :=(\beta _{1},\beta _{2},\beta _{3})$. We estimate $\theta :=(\delta
,\gamma ,\beta ,\varpi )$ through maximizing the objective function
\begin{equation}
Q_{n}(\vartheta ):=Q_{n1}(\vartheta )+Q_{n2}(\vartheta ), \label{app_Qn}
\end{equation}
where
\begin{equation}
Q_{n1}(\vartheta ):=\frac{1}{n}\sum_{i=1}^{n}
\sum_{t=2}^{T_{i}-1}y_{it}(y_{it+1}-y_{it-1})\cdot \mathds{1}\{I_{it}>\sigma
_{n}\}\cdot \mathds{1}\{u_{it}(\vartheta )>0\} \notag \label{app_Qn1}
\end{equation}
with $u_{it}(\vartheta
):=r(y_{it}-y_{it-2})+d(D_{2014,t+1}-D_{2014,t-1})+(x_{it+1}-x_{it-1})^{
\prime }b+w(I_{it+1}-I_{it-1}),$ $\vartheta :=\left( r,d,b,w\right) $ , and
\begin{equation}
Q_{n2}(\vartheta ):=\frac{1}{n}\sum_{i=1}^{n}
\sum_{t=2}^{T_{i}-1}(1-y_{it})(y_{it+1}-y_{it-1})\cdot \mathds{1}
\{I_{it}<-\sigma _{n}\}\cdot \mathds{1}\{u_{it}(\vartheta )>0\}, \notag
\label{app_Qn2}
\end{equation}
extending (\ref{eq:Qn1}) and (\ref{eq:Qn2}) respectively for longer and
unbalanced panels in which the number of waves being observed varies across
individuals $i$ ($=:T_{i}$). Since $\theta $ can only be identified up to
scale, we restrict the search of $\hat{\theta}_{n}$ on the unit sphere.
\end{comment}
By construction, procedure (\ref{app_Qn1}) only uses the subsample of
individuals who can be observed for at least four consecutive waves. After
dropping observations with missing values, our sample consists of 14,880
individuals satisfying this criterion. The panel is unbalanced with $3\leq
T_i\leq 9$ using the notation in previous sections. In total, we have 65,603
observations, among which about 7.36\% observations are ``switchers'' that
are useful for either ours or \citeauthor{honore-k}'s (\citeyear{honore-k})
estimators. As in Section \ref{sec_simulation}, we use $n^{*}$ to denote the number of ``switchers''.
We choose $\sigma _{n}=c\cdot \widehat{\text{std}(I_{it})}\sqrt{\log n^{*}/2.95}$ with $
c=1.0$ and $1.1$ to implement our MS estimation and
report the results in Table \ref{tab_res}. We provide summary statistics for the sub-sample of switchers with $I_{it}>\sigma_n$ in Table \ref{tab_sum_eff} of Appendix \ref{appendix0_33}. We also conducted estimations using $\sigma_n$ with $c=0.5, 0.7,$ and $0.9$. While these results show patterns similar to Table \ref{tab_res}, they highlight the bias-variance trade-off inherent in choosing the tuning parameter, a common challenge in semiparametric methods. These additional results and their discussion are included in Appendix \ref{appendix0_33}.
In addition to the estimates of $\theta$, we also try calculating the 90\% and 95\% confidence intervals
(CIs) for $\theta $ using the $m$-out-of-$n$ bootstrapping. Here we sample $n$ individuals (clusters) to create the bootstrap sample. The main
difficulty in implementing this (or alternative sampling-based) method is
that Theorem \ref{T:limiting_dist} does not give an analytical convergence
rate for the estimator $\hat{\theta}$ due to the unknown tail probabilities
of $I_{it}$. We apply the method proposed in Remark 3 of \cite{LeePun2006}
to solve this problem; that is, assume $\hat{\theta}_n$ has convergence rate of $
n^{\lambda }$ and obtain an estimate $\hat{\lambda}$ of $\lambda $ using a
double $m$-out-of-$n$ bootstrapping procedure with two bootstrap sample
sizes $m_{1}=n^{\rho _{1}}$ and $m_{2}=n^{\rho _{2}}$ for $\rho _{1},\rho
_{2}\in (0,1)$. The 90\% and 95\% CI reported in Table \ref{tab_res} are
calculated with $B=500$ bootstrap replications, $m=n^{7/8}$, and $\hat{
\lambda}=0.309$ (obtained with $\rho _{1}=6/7$ and $\rho _{2}=7/8$).
\begin{table}[htbp]
\caption{Estimates of Preference Coefficients}
\label{tab_res}\centering
\begin{tabular}{clLrcrc}
\hline\hline
& Variable & \multicolumn{1}{c}{Estimate} & \multicolumn{2}{c}{[90\% Conf.Int.]} & \multicolumn{2}{c}{[95\% Conf.Int.]} \\
\midrule
\multirow{5}[2]{*}{$c = 1.0$} & $y_{it-1}$ & 5.275^{\ast\ast} & 0.714 & 20.664 & 0.278 & 21.400 \\
& $\textrm{Above30}_{it}$ & 5.190^{\ast\ast} & 0.319 & 20.286 & 0.016 & 21.465 \\
& $\textrm{Age}_{it}$ & -0.465 & -11.384 & 7.298 & -11.953 & 9.039 \\
& $\textrm{GCC}_{it}$ & -0.317 & -11.140 & 9.127 & -11.645 & 9.738 \\
& $D_{2014,t}$ & -1.548 & -13.747 & 5.664 & -14.276 & 6.648 \\
\hline
\multirow{5}[2]{*}{$c = 1.1$} & $y_{it-1}$ & 5.140^{\ast\ast} & 1.008 & 20.042 & 0.370 & 20.927 \\
& $\textrm{Above30}_{it}$ & 5.227^{\ast\ast} & 1.000 & 19.524 & 0.455 & 20.581 \\
& $\textrm{Age}_{it}$ & -0.505 & -10.909 & 7.203 & -11.470 & 8.683 \\
& $\textrm{GCC}_{it}$ & -1.770 & -13.725 & 4.950 & -14.299 & 6.573 \\
& $D_{2014,t}$ & -1.578 & -13.274 & 5.303 & -13.598 & 6.582 \\
\hline\hline
\end{tabular}
\end{table}
We can see from Table \ref{tab_res} that the estimation results are similar
for the two tuning parameters. Therefore, the following discussion on the empirical
results will be mainly based on the estimates obtained with $c=1.0$. The insignificant coefficient indicates that living in GCC may not affect people's
willingness to buy PHC. The significant positive coefficient
on $y_{it-1}$ demonstrates the strong state persistence of PHC, which explains
why we can only observe a small percentage of switchers in the data.
Surprisingly, people's decision to buy PHC is hardly influenced by age. For
the two policy variables, the large positive coefficient on $\text{Above30}
_{it}$ confirms that the LHC policy is a strong incentive for people to buy
PHC, while the change in MLS income threshold does not exhibit a strong
impact represented by the coefficient on $D_{2014,t}$. An intuitive
explanation for the latter is that although MLS promoted PHC purchases when
it was introduced in 1997--1998, the subsequent adjustments of its income
threshold only affected a small group of people whose incomes were near the
threshold.
We end this section with some remarks. First, our approach is more suitable
for data with a relatively large proportion of ``switchers'' which make up
the effective sample for the estimator. Second, to implement our method, the
model should have a continuous covariate with large support and ideally weak
dependence on other included covariates. Third, in the
absence of knowledge (or at least a good estimate) of the free-varying
variable's tail probabilities, the asymptotics of our estimator derived in
Section \ref{Sec:asymp_dist} cannot provide a ``rule of thumb'' for choosing
optimal tuning parameter $\sigma_{n}$. Perhaps a practical way is to try
different $\sigma_{n}$'s, use \citeauthor{LeePun2006}'s (
\citeyear{LeePun2006}) proposed method (or other similar methods) to
estimate the convergence rates, and pick the $\sigma_{n}$ that gives the
fastest (estimated) rate. The last remark is for the $m$-out-of-$n$
bootstrap inference. The choice of the bootstrap sample size $m$ is the key
issue. Remark 1 of \cite{LeePun2006} provides some existing data-driven
methods. However, none of them can confirm an (asymptotically) optimal
choice of $m$ in nonstandard M-estimation like ours. Theoretical research on
this topic is necessary, but this is beyond the scope of the current paper.
\begin{comment}
We restrict our attention to those estimates where both inference procedures
do agree. First, the estimated coefficients on the lagged labor force
participation are significantly positive for all three samples. This
indicates that the labor force participation decision is sticky over time, so
ignoring the state dependence may lead to mis-specification of the model. The
estimated coefficients on `` Health Shock''
are significantly negative for all samples. This implies that temporary
deterioration in health status does have a negative effect on labor force
participation, which is not surprising at all. The estimated coefficients
before `` log(Income)'' are only significant
(and positive) for the male sample. The result is also in line with our
economic intuition. Family income is more likely to be a crucial determinant
of labor force participation for males, rather than females. The estimated
coefficients before the long-run health shock (``Activity
Limiting Condition'') are not significant, at 5\% level for
all three samples. This suggests that permanent changes in health status may
not have much effect on the decision to switch from working to not-working, or
the other way around. The estimates on `` Unemployment
Rate'' are not significant for all samples. This basically
means that an individual's working decision is somewhat independent of the
unemployment rate, a macro-level variable indicating general conditions in the
labor market.
Our results agree with the fixed effects estimates in
\cite{DamrongplasitEtal2018}, in terms of the signs of the preference
coefficients. However, the magnitudes of these coefficients are very
different. For example, their results indicate that the effect of
`` Health Shock''\ on an individual's labor
force participation is much smaller in magnitude than that of the state
dependence (less than 1/3), while our results show that their effects are
quite comparable in magnitude. We note that their estimates were obtained
using HK's parametric (conditional Logit) estimator, which
might suffer from mis-specification. Our results can be thought of as a robust
check of their results. \cite{DamrongplasitEtal2018} also obtained the random
effects estimates, and demonstrated that the fixed effects estimates are more
reasonable. We conjecture that one could arrive at similar random effects
estimates using our samples.
As a final note, recall that the preference coefficients are only identified up to scale. The magnitudes of these coefficients are difficult to interpret.
Our semiparametric estimates are most useful if several coefficients are
included in the regression and coefficient estimates are computed to compare
the relative effects of changes in regressors on the choice.
As a final note, it is worth noting that the procedure (\ref{app_Qn}) does not
make full use of the data. In fact, our identification strategy works
with ``switchers'' who make different choices in any two non-initial,
non-adjacent periods (waves). For a general sequence of observed choices
$\{y_{i0},y_{i1},...,y_{iT}\}$, any subsequence $\{y_{it-1},y_{it},...,y_{i\tau}\}$
with $t\geq1$ and $t+1<\tau\le T$ satisfying $y_{it}\neq y_{i\tau}$
can be used to construct objective function of the form similar to
(\ref{app_Qn}). In total, there should be $(T-1)(T-2)/2$ such subsequences to
use. In our application, this number is 78, but (\ref{app_Qn}) only uses 12 of
them, i.e., all consecutive subsequences of length 4. Using longer subsequences can increase the effective sample (number of ``switchers'') and hence improve the finite-sample efficiency of the estimation procedure, especially for applications showing strong state dependence (like here). However, using
longer subsequences involves choosing tuning parameter $\sigma_{n}$
for each $I_{i\iota}$ with $t<\iota<\tau$, which may be cumbersome in practice. In an
ideal case where the distribution of $I_{it}$ is time stationary,
one can apply the same $\sigma_{n}$ to all waves.
\end{comment}
\section{Conclusions}\label{Sec:Conlcusions}
This paper proposes new identification and estimation methods for a class of distribution-free dynamic panel data binary choice models that is first studied in \cite{honore-k}. We show that in the presence of a free-varying continuous covariate with unbounded support, an ``identification at infinity'' strategy in the spirit of \cite{chamberlain1986asymptotic} enables the point identification of the model coefficients without the need of element-by-element matching of covariates over time, in contrast to the method proposed in \cite{honore-k}. This property makes our methods more practical for models with many covariates or important covariates whose support may not overlap over time. Our identification arguments motivate a conditional maximum score estimator that is proven to be consistent and with the convergence rate independent of the model dimension. However, the asymptotic distribution of the proposed estimator is non-Gaussian, in line with well-established literature on cube-root asymptotics. We suggest valid bootstrap methods for conducting statistical inference. The results of a Monte Carlo study demonstrate that our estimator performs adequately in finite samples. Lastly, we use the HILDA data to investigate the demand for private hospital insurance in Australia.
This paper leaves some open questions for future research. For instance, although we suggest several theoretically feasible bootstrap inference methods in Section \ref{Sec:Asymptotics}, their asymptotic validity, finite-sample performance, and implementability (e.g., choice of tuning parameters) are not examined. Alternatively, one can also investigate whether it is possible to achieve a faster rate of convergence and obtain an asymptotically normal distribution by combining \citeauthor{horowitz}'s (\citeyear{horowitz}) and \citeauthor{andrews}'s (\citeyear{andrews}) methods to smooth the sample objective function.
\nocite{HILDA,HILDA2}
\bibliography{references}