\section{Introduction}

In this paper, we propose new identification and estimation methods for
panel data binary choice models with fixed effects and \textquotedblleft
dynamics\textquotedblright\ (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. In panel data literature, $%
y_{it-1}$ is often called \textquotedblleft state
dependence\textquotedblright , and $\alpha _{i}$ is referred to as
\textquotedblleft unobserved heterogeneity\textquotedblright\ or
\textquotedblleft suprious\textquotedblright\ state dependence (see \cite%
{Heckman1981b, Heckman1981a}). In model (\ref{model}), $%
(y_{it},x_{it},z_{it})$ along with the \textquotedblleft initial
status\textquotedblright\ $y_{i0}$ are observed in the data, while $\alpha
_{i}$ and $\epsilon _{it}$ are unobservable to the econometrician.\footnote{%
Though $y_{i0}$ is assumed to be observed, 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 \textquotedblleft short\textquotedblright\ panel
settings, i.e., $n\rightarrow \infty $ an $T<\infty $. Since one can only
identify $\theta $ up to scale$,$ we normalize $\left\Vert \theta
\right\Vert =1,$ where $\left\Vert \cdot \right\Vert $ is the Euclidean norm$%
.$

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 i.i.d. Logistic $\epsilon_{it}$'s, \cite{rasch} and \cite{andersen} show that, in the absence
of $y_{it-1}$, this model is identified and can be estimated by conditional
maximum likelihood method. \cite{Chamberlain2010} proves that the Logit
assumption is crucial for such models to have 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} take one step in this direction to develop a conditional maximum
likelihood estimator for model (\ref{model}) with covariates, and \cite{hahn}
examines the semiparametric efficiency of this estimator. Many recent
advances have been made for variations of the Logit model. See, e.g., \cite%
{johnson}, \cite{BartolucciNigro2010,BartolucciNigro2012}, \cite{bonhomme},
\cite{honore-w}, \cite{dobronyi-gu-kim}, \cite{kiazawa}, and \cite{muris},
to name a few.

Without making distributional assumptions on $\epsilon _{it}$, %
\citet{manski-87} establishes the semiparametric identification of model (%
\ref{model}) with covariates, but no $y_{it-1}$. \cite{honore-k} extend this
approach to a model with both $y_{it-1}$ and covariates, 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%
{OuyangYang2020binary} show that this curse of dimensionality can be
overcome if we impose certain serial dependence conditions on the covariates
and observe an extra time period.

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} show that the model (with or without $y_{it-1}$%
) can be point identified if $z_{it}$ satisfies certain exclusion
restrictions (excluded regressor), even if covariates are endogenous. More
recently, \cite{ChenEtal2019} revisit this method and discuss 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{KhanEtal2021} provide a distinct
partial identification result under even weaker restrictions and prove that
point identification exists in many interesting situations.

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 that is first studied 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 with no 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 show that our estimator is consistent with a slower than
cube-root-$n$ rate of convergence that is independent of the number of
observed covariates and has a non-standard limiting distribution. The
finite-sample performance and implementability of the proposed estimator are
examined using simulated and real-world data.

The rest of this paper is organized as follows. Section \ref%
{sec:identification} establishes the identification of $\theta$, based on
which, Section \ref{Sec:Estimation} provides an MS estimator. 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 in
Section \ref{sec_application} using the HILDA data. Lastly, Section \ref{Sec:Conlcusions} concludes the paper and
briefly discusses some future research directions. Proofs and additional
simulation results are collected in the \hyperref[appendix]{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_{2},...,x_{t}\right) $ and  $%
y^{t}:=\left(y_{1},y_{2},...,y_{t}\right) $.
For two random vectors, $u$ and $v$, the notation $u\overset{d}{=}v|\cdot $
means that $u$ and $u$ 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$), and $a_{n}\asymp b_{n}$ means that
there exist two constants $0<c_{1}\leq c_{2}<\infty $ such that $%
c_{1}b_{n}\leq a_{n}\leq c_{2}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)   \notag  \label{p1}
\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) \notag
\label{likelihood} \\
 =&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)   \notag \\
 =&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)   \notag \\
=& \cdots   =\prod_{t=1}^{T}P\left( y_{it}|\alpha
_{i},y_{i}^{t-1},x_{i}^{T},z_{i}^{T}\right)   \notag \\
 =&\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}}  
\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 latter requires that we know $z_{it}$ is in the model and has
positive effect on the choice probability of $y_{it}=1$. It is straightforward to apply
our method to longer panels and the case with $\varpi <0$ is symmetric. In addition, we will suppress the subscript $i$
to simplify notation whenever it is clear from the context that all
variables are for each individual.\footnote{%
Here we also implicitly assume balanced panel to simplify exposition. The
identification and estimation methods presented in this paper can be applied
to models with unbalanced panels as long as the unbalancedness is not the
result of endogenous attrition.}

Consider two choice histories
\begin{align*}
C& =\{y_{0}=d_{0},y_{1}=0,y_{2}=d_{2},y_{3}=1\}, \\
D& =\{y_{0}=d_{0},y_{1}=1,y_{2}=d_{2},y_{3}=0\},
\end{align*}%
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},x^{T},z^{T})  \notag \\
=& 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}))  \notag \\
& \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}),  \notag
\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},x^{T},z^{T})  \notag \\
=& 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})  \notag \\
& \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})).  \notag
\end{align}

Suppose the support of $z_{2}$ is unbounded either above or below. If $%
d_{2}=1$, these probabilities satisfy
\begin{equation}
\lim_{z_{2}\rightarrow +\infty }\frac{P(C|\alpha ,y_{0},x^{T},z^{T})}{%
P(D|\alpha ,y_{0},x^{T},z^{T})}=\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{equation}%
Else, if $d_{2}=0$, they satisfy
\begin{equation}
\lim_{z_{2}\rightarrow -\infty }\frac{P(C|\alpha ,y_{0},x^{T},z^{T})}{%
P(D|\alpha ,y_{0},x^{T},z^{T})}=\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_2}
\end{equation}%
The idea of the above is to make negligible the effect of $y_{1}$ on $y_{2}$ by
letting $\left\vert z_{2}\right\vert \rightarrow \infty .$

Suppose $F_{\epsilon |\alpha }(\cdot )$ is strictly increasing. Then, when $%
d_{2}=1$, equation (\ref{eq:prob_ratio_1}) implies
\begin{equation}
\text{sgn}\left\{ P(C|\alpha ,y_{0},x^{T},z^{T})-P(D|\alpha
,y_{0},x^{T},z^{T})\right\} =\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{equation}%
when $z_{2}\rightarrow +\infty $, where $\text{sgn}\{\cdot \}$ is the \emph{%
sign function} which is equal to 1 if the expression inside the brackets is
strictly positive, 0 if the expression inside the brackets is zero, and $-1$
if the expression inside the brackets is strictly negative. Similarly, when $%
d_{2}=0$, equation (\ref{eq:prob_ratio_2}) implies
\begin{equation}
\text{sgn}\left\{ P(C|\alpha ,y_{0},x^{T},z^{T})-P(D|\alpha
,y_{0},x^{T},z^{T})\right\} =\text{sgn}\left\{ \gamma
(d_{2}-d_{0})+(x_{3}-x_{1})^{\prime }\beta +\varpi (z_{3}-z_{1})\right\}
\label{eq:iden_ineq2}
\end{equation}%
when $z_{2}\rightarrow -\infty $.

Equations (\ref{eq:iden_ineq1}) and (\ref{eq:iden_ineq2}), therefore,
provide a way to predict $y_{3}-y_{1}$ (whenever $y_{3}\neq y_{1}$) in the
spirit of \cite{manski-87}: In particular, if the likelihood of observing
event $C$ is larger than the likelihood of observing event $D$, then it is
reasonable to expect $y_{3}-y_{1}=+1$; otherwise, it is reasonable to expect
$y_{3}-y_{1}=-1$. Since, when $z_{2}$ is large enough in its absolute value,
the likelihood of observing event $C$ is larger than the likelihood 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$, we can
use the sign of this expression to predict $y_{3}-y_{1}$.

\begin{remark}
It is clear that the derivation of (\ref{eq:prob_ratio_1})--(\ref%
{eq:iden_ineq2}) relies on the assumption that the marginal distribution of $%
z_{2}$, conditional on $\alpha $, $y_{0}$, and all other covariates, has an
unbounded support. \cite{honore-k} identify model (\ref{model}) under the
restriction that $(x_{3}-x_{2},z_{3}-z_{2})$ has a support in some open
neighborhood of zero, conditional on $\alpha $, $y_{0}$, and all other
covariates.\footnote{%
The $x_{t}$ in \cite{honore-k} corresponds to our $(x_{t},z_{t})$.} Specifically, for the model with $T=3$, \cite{honore-k} consider the
following choice histories:
\begin{align*}
A& =\{y_{0}=d_{0},y_{1}=0,y_{2}=1,y_{3}=d_{3}\}, \\
B& =\{y_{0}=d_{0},y_{1}=1,y_{2}=0,y_{3}=d_{3}\},
\end{align*}%
where $d_{0},d_{3}\in \{0,1\}$. Applying argument similar to (\ref%
{eq:prob_ratio_1})--(\ref{eq:prob_ratio_2}) yields
\begin{equation}
\frac{P(A|\alpha ,y_{0},x^{T},z^{T},x_{2}=x_{3},z_{2}=z_{3})}{P(B|\alpha
,y_{0},x^{T},z^{T},x_{2}=x_{3},z_{2}=z_{3})}=\frac{1-F_{\epsilon |\alpha
}(\alpha +\gamma d_{0}+x_{1}^{\prime }\beta +\varpi z_{1})}{1-F_{\epsilon
|\alpha }(\alpha +\gamma d_{3}+x_{2}^{\prime }\beta +\varpi z_{2})}\times
\frac{F_{\epsilon |\alpha }(\alpha +\gamma d_{3}+x_{2}^{\prime }\beta
+\varpi z_{2})}{F_{\epsilon |\alpha }(\alpha +\gamma d_{0}+x_{1}^{\prime
}\beta +\varpi z_{1})}.  \notag  \label{eq:ratio}
\end{equation}%
Then this expression implies
\begin{align}
& \text{sgn}\left\{ P(A|\alpha
,y_{0},x^{T},z^{T},x_{2}=x_{3},z_{2}=z_{3})-P(B|\alpha
,y_{0},x^{T},z^{T},x_{2}=x_{3},z_{2}=z_{3})\right\}  \notag \\
=& \text{sgn}\left\{ \gamma (d_{3}-d_{0})+(x_{2}-x_{1})^{\prime }\beta
+\varpi (z_{2}-z_{1})\right\} ,  \label{eq:HK}
\end{align}%
based on which the point identification can be established. Similar to (\ref%
{eq:iden_ineq1})--(\ref{eq:iden_ineq2}), the relationship above provides a
way to predict $y_{2}-y_{1}$, whenever $y_{2}\neq y_{1}$. The argument to
equation (\ref{eq:HK}) imposes the restriction on the support of $%
(x_{3}-x_{2},z_{3}-z_{2})$, which essentially excludes covariates with
non-overlapping support over time (e.g., time trend or time dummies for
controlling fixed time effects) and may encounter curse of dimensionality
when there are many relevant covariates. Our method does not suffer from
these two issues because we did not match $(x_{t},z_{t})
$ overt time.
\end{remark}

A natural way to build a population objective function based on equation (%
\ref{eq:iden_ineq1}) is to define
\begin{align}
\bar{Q}_{1}(b,r,w):=\lim_{\sigma \rightarrow +\infty }\mathbb{E}& \left[
\left( P(C|\alpha ,y_{0},x^{T},z^{T})-P(D|\alpha ,y_{0},x^{T},z^{T})\right)
\right.  \notag \\
& \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$, and similarly, based on equation (\ref%
{eq:iden_ineq2}), define
\begin{align}
\bar{Q}_{2}(b,r,w):=\lim_{\sigma \rightarrow +\infty }\mathbb{E}& \left[
\left( P(C|\alpha ,y_{0},x^{T},z^{T})-P(D|\alpha ,y_{0},x^{T},z^{T})\right)
\right.  \notag \\
& \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$. A moment's thought reveals that $\bar{Q}_{1}(\beta
,\gamma ,\varpi )\geq \bar{Q}_{1}(b,r,w)$ and $\bar{Q}_{2}(\beta ,\gamma
,\varpi )\geq \bar{Q}_{2}(b,r,w)$ for all
\begin{equation*}
(b,r,w)\neq (\beta ,\gamma ,\varpi )=:\theta .
\end{equation*}
Further proving $\theta $ is the unique maximum of objective functions (\ref%
{eq:pop_obj}) or (\ref{eq:pop_obj_2}) establishes the point identification
of these coefficients, for which the following conditions suffice.

\begin{assumptionp}{A}
For all $\alpha$ and $s,t\in\mathcal{T}:=\{1,2,3\}$ ($T= 3$), the following
conditions hold:

\begin{enumerate}
\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) the CDF of $\epsilon_t$ is absolutely
continuous with support $\mathbb{R}$ conditional on $\alpha$.

\item[A2] $z_{2}$ has unbounded support conditional on $(\alpha
,y_{1},x^{T},z_{1},z_{3})$.

\item[A3] $z_{3}-z_{1}$ has bounded from above and almost everywhere (a.e.)
positive Lebesgue density on $\mathbb{R}$ conditional on $(\alpha
,x_{3}-x_{1})$ and $\{z_{2}>\sigma \}\cup \{z_{2}<-\sigma \}$ as $\sigma
\rightarrow +\infty $.

\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.

\item[A5] Let $\Theta:=\{\vartheta:=(r,b,w)\in\mathbb{R}^{p+2}:\Vert\vartheta\Vert=1\cap \{w\geq\iota\}\}$ and $\iota$ is some positive constant. $\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 \cite{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 enables the approach of ``identification at infinity'' we use
and, together with Assumption A1, establishes identifying equations (\ref%
{eq:iden_ineq1}) and (\ref{eq:iden_ineq2}).

The support and continuity restrictions on $z_3-z_1$ imposed in Assumption
A3 are common for the family of maximum-score type estimators, which are
required to achieve the point identification instead of a set
identification. See, e.g., \cite{manski-75, manski-85, manski-87}, \cite%
{horowitz}, \cite{honore-k}, \cite{Fox2007}, \cite{ShiEtal2018}, \cite{YanYoo2019}, and \cite%
{KOT2021} among many others. 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.

Our identification results are stated in the following theorem and proved in
Appendix \ref{appendixA}.

\begin{theorem}
\label{T:identify} Suppose Assumption A hold. 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}) built
respectively on identifying equations (\ref{eq:iden_ineq1}) and (\ref%
{eq:iden_ineq2}) 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})
satisfying 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 either above or 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}%
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 make different
choices in periods 1 and 3 are used in the estimation. Besides, only
``switchers'' whose choice changes in periods 0 and 2 are relevant in
estimating (identifying) $\gamma$.

Our proposed estimator $\hat{\theta}_n$ has two advantages compared to %
\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}. Since 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.

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 $P(z_{2}>\sigma
_{n})\asymp e^{-\sigma _{n}}$, then $\sigma _{n}=O(\log \log (n/\log n))$
satisfies Assumption B2(ii). However, when the distribution of $z_{2}$ has a
(too) thin right tail $P(z_{2}>\sigma _{n})\asymp e^{-e^{\sigma _{n}}}$%
, the same $\sigma _{n}$ gives $nP(z_{2}>\sigma _{n})/\log n=O(1)$ violating
Assumption B2(ii). Note that $nP(z_{2}>\sigma _{n})$ is essentially the
\textquotedblleft effective sample size\textquotedblright\ for our proposed
procedure.\footnote{%
Precisely speaking, the effective sample for (\ref{eq:Qn1}) should be
observations with $z_{2}>\sigma _{n}$, $y_{2}=1$, and $y_{1}\neq y_{3}$.
Here we suppress the dependence of effective sample on $y^{T}$ to ease the
discussion.} For example, consider sub-Gaussian\footnote{%
See, e.g., Chapter 2 of \cite{WainwrightHDS} for a gentle introduction to
sub-Gaussian and sub-exponential distributions.} $z_{2}$ with a right tail $%
P(z_{2}>\sigma _{n})\asymp e^{-\sigma _{n}^{2}}$ and $\sigma _{n}=\sqrt{%
\log \log n}$. Then objective function (\ref{eq:Qn1}) uses $nP(z_{2}>\sigma
_{n})\asymp n/\log n$ observations, while for the same $\sigma _{n}$ but $%
z_{2}$ with heavier tails (e.g., sub-exponential tail $P(z_{2}>\sigma
_{n})\asymp e^{-\sigma _{n}}$ or geometrically-decaying (Pareto) tail $%
P(z_{2}>\sigma _{n})\asymp \sigma _{n}^{-2}$), $nP(z_{2}>\sigma _{n})$ tends
to be larger.\footnote{%
To see this, note that the sample size of $n/\log n$ is associated with $%
\sigma _{n}=\log \log n$ and $\sigma _{n}=\sqrt{\log n}$ for these two
distributions, respectively.} As will be clear in Theroem \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.

The theorem below states that the proposed procedure described in (\ref%
{eq:Qn1})--(\ref{eq:Qn})\footnote{%
Some straightforward adjustments to Assumption B are necessary for
establishing similar results for $\hat{\theta}_n$ associated with (\ref%
{eq:Qn2}) or (\ref{eq:Qn}).} gives 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}

\begin{comment}
\textcolor{red}{Introducing a smooth function to replace the $\mathbf{1}[z_{i2}>\sigma_{n}]$
may be useful in deriving the limiting distribution of $\hat{\theta}$.
An example of such a smooth function was given by formula (4.1) of
\cite{andrews}. Do we really need this? Read \cite{SeoOtsu2018}.}

\emph{[An alternative to this estimator would maximize something like:
\[
\begin{gathered}
	\sum_{i=1}^nK\left(\frac{z_{i2}-q_n}{\sigma_n}\right)y_{i2}(y_{i3}-y_{i1})\; \text{sgn}\big\{g(y_{i2}-y_{i0})+b(x_{i3}-x_{i1})+(z_{i3}-z_{i1})\big\},
\end{gathered}
\]
when $z_{i2}$ is continuous, where $K$ is a kernel, but I find this estimator to be more difficult to work with because it depends on the relative speeds at which the bandwidths $\sigma_n$ and $q_n$ converge.}

\emph{Could alternatively propose a smooth version like Horowitz, but that has more bandwidths.}

\emph{The case when $z_{i2}$ is unbounded both above and below would yield an estimator that maximizes:
\begin{equation}
\begin{gathered}
	\sum_{i=1}^n\big[y_{i2}\mathds{1}\{z_{i2}>\sigma_n^+\} + (1-y_{i2})\mathds{1}\{z_{i2}<\sigma_n^-\}\big]\;(y_{i3}-y_{i1})\\ \mathds{1}\big\{g(y_{i2}-y_{i0})+b(x_{i3}-x_{i1})+(z_{i3}-z_{i1})>0\big\}.
\end{gathered}
\end{equation}
]}
\end{comment}

\subsection{Asymptotic Distribution}

\label{Sec:asymp_dist} We move on to study the asymptotic distribution of
the estimator $\hat{\theta}_n$. Before presenting additional technical
conditions and the main results, we introduce some new notations to
facilitate exposition:

\begin{itemize}
\item[-] $h_{n}:=P(z_{2}>\sigma_{n}|y_{2}=1)$.\footnote{%
Note that under Assumptions A1 and A3, $P(z_{2}>\sigma_{n})\asymp
P(z_{2}>\sigma_{n}|y^{T})$, and so $h_n$ satisfies Assumption B2(ii).}

\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}%
\left((nh_{n})^{-2/3}\right)$.

\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 neighbourhood 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] There exist some constants $c>0$ and $\varrho \in (1/6,1/2)$ such
that $h_{n}\geq cn^{1/2-3\varrho }$.
\end{enumerate}
\end{assumptionp}

Assumption C1 is standard in the literature and defines our MS estimator.
Assumption C2 is an implication of Assumptions A1--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, plays an important role
in calculating the expected value of the limiting distribution of the
estimator $\hat{\theta}_{n}$.

The smoothness conditions imposed in Assumption C5(i) are also standard in
the literature. Assumption C5(ii) is made to simplify the proof process and
can be relaxed.\footnote{%
Here what we really want 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 observed regressor $z_t$ and unobserved error $\epsilon_t$.%
\footnote{%
Assumption C5(iii) implies $P(\epsilon_2\geq\varsigma+%
\sigma_n)=O((nh_{n})^{-1/3})$ for all $\varsigma\in\mathbb{R}$.} As shown in
the proof of Theorem \ref{T:limiting_dist} presented in Appendix \ref%
{appendixB}, this assumption ensures that the bias of the estimator $\hat{%
\theta}_{n}$ shrinks sufficiently fast. 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), provides guidance on
choosing the tuning parameter $\sigma_{n}$. 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 can be understood as $\sigma_{n}$ should
not increase too fast as $n\rightarrow\infty$ so that we have enough
effective observations to control the variance of $\hat{\theta}_{n}$. It is
easy to verify that for any distribution of $z_{2}$ with a right tail
heavier than sub-Gaussian, the choice of $\sigma_{n}=\sqrt{\log\log n}$ can
give an $h_{n}$ satisfying Assumption C6 (with e.g., $\varrho=1/4$). On the
other hand, Assumption C5(iii) indicates that $\sigma_{n}$ should grow
sufficiently fast as $n\rightarrow\infty$ to lower the bias of $\hat{\theta}%
_{n}$. However, there is no ``rule of thumb'' for selecting optimal $%
\sigma_{n}$ since this requires the knowledge of relative (unknown) tail
behavior of $z_{t}$ and $\epsilon_{t}$.\footnote{%
Let $\tau_{z}(\sigma_{n}):=P(z_{2}>\sigma_{n}|y_{2}=1)$ and $%
\tau_{\epsilon}(\sigma_{n}):=P(\epsilon_{t}>\sigma_{n})$. Suppose we have
identified and (nonparametrically) estimated $\tau_{z}(\cdot)$ from the
data. Since $h_{n}=\tau_{z}(\sigma_{n})$, it is easy to find a candidate $%
\sigma_{n}$ that satisfies Assumption C6 but diverges slowly as $%
n\rightarrow\infty$. If we also (magically) knew the functional form of $%
\tau_{\epsilon}(\cdot)$, then we can move on to check if this $\sigma_{n}$
also satisfies Assumption C5(iii), i.e., there is some constant $c>0$ such
that $\tau_{\epsilon}(\sigma_{n})\leq c(nh_{n})^{-1/3}$ for $n$ large
enough. If the result is ``true'', we can choose this $\sigma_{n}$;
otherwise, we need to pick a $\sigma_{n}$ diverging faster to repeat this
process until we get a ``true'' result. But it is clear that this selection
process is not feasible as it is unrealistic to know $\tau_{\epsilon}(\cdot)$%
.}

The above conditions suffice for asymptotic distribution of our proposed
procedure (\ref{eq:Qn1})--(\ref{eq:Qn}) as stated in the following theorem,
proved 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}

\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 to have a $z_{t}$ whose distribution has thicker tails relative
to $\epsilon_{t}$. To see this, consider any eligible $\sigma_{n}$
satisfying both Assumptions C5 and C6. Suppose $h_{n}^{\lambda}%
\asymp(nh_{n})^{-1/3}$ for some $\lambda>0$. Loosely speaking, $%
\lambda\in(0,1)$ ($\lambda>1$) means that $z_{t}$ has thinner (thicker)
tails relative to $\epsilon_{t}$. Then we have $h_{n}\asymp
n^{-1/(3\lambda+1)}$ and hence $(nh_{n})^{-1/3}\asymp
n^{-\lambda/(3\lambda+1)}$. Thus, for any eligible $\sigma_{n}$, the
convergence rate of $\hat{\theta}_{n}$ increases in $\lambda$ and approaches
$n^{-1/3}$ for large $\lambda$.
\end{remark}

\begin{remark}
Theorem \ref{T:limiting_dist} indicates that the proposed estimator $\hat{%
\theta}$ has a slower than cube-root $n$ rate of convergence and its
asymptotic distribution is not normal. As a result, usual 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}), \emph{m}-out-of-\emph{n} bootstrap (%
\cite{LeePun2006}), the numerical bootstrap (\cite{HongLi2020}), and a
model-based bootstrap with modified objective function (\cite%
{cattaneo2020bootstrap}), among other procedures. \cite{OuyangYang2020match,
OuyangYang2020binary} 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 may apply to the estimator of this paper. However, it is
beyond the scope of this paper to rigorously develop a bootstrap inference
algorithm and justify its asymptotic validity. We leave this work to future
research.
\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 start by considering a benchmark design similar to that used in \cite%
{honore-k}, but we add additional covariate $z_{it}$ and a time trend that
\cite{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 (0-2)+\beta
_{1}x_{i0,1}+\varpi z_{i0}\geq \epsilon _{i0}\right\} , \\
y_{it}& =\mathds{1}\left\{ \alpha _{i}+\delta (t-2)+\gamma y_{i,t-1}+\beta
_{1}x_{it,1}+\varpi z_{it}\geq \epsilon _{it}\right\} ,\text{ \ }t\in
\left\{ 1,2,3\right\} ,
\end{align*}%
where we normalize the coefficients to norm 1 and $\gamma =\beta _{1}=\varpi =2/\sqrt{13}\approx 0.555$ and $\delta =1/\sqrt{13}\approx 0.277$. We
consider two cases. For both cases, 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. All covariates are i.i.d. across $i$ and $t.$
$x_{\cdot ,1},z_{\cdot },$ and $\epsilon _{\cdot }$ are independent of each
other. In the first case, we set $z_{it}\overset{d}{\sim }N(0,1),$ and
denote it as \textquotedblleft Norm\textquotedblright . In the second case,
we set $z_{it}\overset{d}{\sim }\text{Laplace}(0,\sqrt{2}/2)$ (with zero
mean and unit variance), and denote it as \textquotedblleft
Lap\textquotedblright . 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 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 (0-2)+\beta _{1}x_{i0,1}+\beta
_{2}x_{i0,2}+\varpi z_{i0}\geq \epsilon _{i0}\right\} , \\
y_{it}& =\mathds{1}\left\{ \alpha _{i}+\delta (t-2)+\gamma y_{it-1}+\beta
_{1}x_{it,1}+\beta _{2}x_{it,2}+\varpi z_{it}\geq \epsilon _{it}\right\} ,%
\text{ \ }t\in \left\{ 1,2,3\right\} ,
\end{align*}%
where we normalize the coefficients to norm 1 and $\gamma =\beta _{1}=\beta _{2}=\varpi =2/\sqrt{17}\approx0.485$ and $\delta =1/%
\sqrt{17}\approx0.243$. 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 cases
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.
Note that the convergence rate of the estimator depends on the relative tail
thickness of $z_{it}$ to other regressors. To investigate the impact of higher
dimension only, we set the variance of $x_{it,1}+x_{it,2}$ in Design 2 the
same as that of $x_{it,1}$ in Design 1.

For the estimation, we restrict the search of the parameters on the unit
sphere. We consider the estimations using $\left\vert z_{i2}\right\vert $
that exceeds
\begin{equation*}
\sigma _{n}=\widehat{\text{std}\left( z_{i2}\right) }\sqrt{\log \log n},
\end{equation*}%
where $\widehat{\text{std}\left( z_{i2}\right) }$ is the sample standard
deviation of $z_{2}$. In practice, one may normalize $z_{it}$ to mean 0 and
variance 1 and set $\sigma _{n}=\sqrt{\log \log n}$. We consider sample
sizes of $n=5000,10000$, and $20000$. All MS estimations are implemented in
R, using the differential evolution algorithm ($\mathtt{DEoptim}$) to attain
a global optimum of the objective function. All the simulation results
presented in this section are based on 500 replications of each sample size.
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 from the simulation results. First, the RMSE 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}$, comparing the
RMSE 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 the RMSE does increase a bit
for $\beta _{1}$, but this is probably due to the lower variance of $x_{\cdot ,1}$ in
Design 2. To investigate the sensitive of the results to $\sigma _{n},$ we
consider $\sigma _{n}=0.9\cdot\widehat{\text{std}\left( z_{i2}\right) }\sqrt{\log
\log n}$ and $\sigma _{n}=1.1\cdot\widehat{\text{std}\left( z_{i2}\right) }\sqrt{%
\log \log n},$ 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 here.

\newpage
\begin{table}[ptb]
\caption{Simulation Results of Design 1}
\label{T:D1}\centering
\begin{tabular}{r|cc|cc|cc|cc}
\hline\hline
& \multicolumn{2}{c|}{$\beta_{1}$} & \multicolumn{2}{c|}{$\gamma$} &
\multicolumn{2}{c|}{$\delta$} & \multicolumn{2}{c}{$\varpi$}\\
& MBIAS & RMSE & MBIAS & RMSE & MBIAS & RMSE & MBIAS & RMSE\\ \hline
$n_1$ & 0.011 & 0.095 & -0.067 & 0.188 & 0.018 & 0.083 & -0.008 & 0.098 \\
Norm  $n_2$ & 0.012 & 0.075 & -0.058 & 0.143 & 0.025 & 0.071 & 0.001 & 0.076 \\
 $n_3$ & 0.015 & 0.061 & -0.055 & 0.116 & 0.024 & 0.058 & 0.007 & 0.060 \\\hline
$n_1$  & 0.005 & 0.082 & -0.052 & 0.157 & 0.014 & 0.073 & 0.001 & 0.075 \\
Lap  $n_2$ & 0.011 & 0.065 & -0.048 & 0.122 & 0.016 & 0.059 & 0.005 & 0.063 \\
 $n_3$ & 0.016 & 0.052 & -0.047 & 0.100 & 0.014 & 0.049 & 0.009 & 0.047 \\\hline\hline
\multicolumn{9}{l}{Note: $n_{1}=5000, n_{2}=10000, n_{3}=20000$.} \\
\end{tabular}%
\end{table}

\begin{table}[ptb]
\caption{Simulation Results of Design 2}\label{T:D2}\centering
\begin{adjustbox}{width=\textwidth}
\begin{tabular}{r|cc|cc|cc|cc|cc}
\hline\hline
& \multicolumn{2}{c|}{$\beta_{1}$} & \multicolumn{2}{c|}{$\beta_{2}$} &
\multicolumn{2}{c|}{$\gamma$} & \multicolumn{2}{c|}{$\delta$} & \multicolumn{2}{c}{$\varpi$}\\
& MBIAS & RMSE & MBIAS & RMSE & MBIAS & RMSE & MBIAS & RMSE & MBIAS & RMSE \\ \hline
$n_1$  & 0.001 & 0.108 & 0.006 & 0.100 & -0.057 & 0.156 & 0.018 & 0.071 & -0.019 & 0.090 \\
Norm $n_2$  & 0.012 & 0.082 & 0.009 & 0.082 & -0.059 & 0.136 & 0.014 & 0.056 & -0.010 & 0.071 \\
$n_3$ & 0.002 & 0.068 & 0.008 & 0.067 & -0.045 & 0.110 & 0.015 & 0.046 & 0.000 & 0.055 \\\hline
$n_1$  & 0.002 & 0.084 & 0.006 & 0.097 & -0.041 & 0.129 & 0.009 & 0.057 & -0.015 & 0.071 \\
Lap $n_2$   & 0.001 & 0.069 & 0.004 & 0.073 & -0.031 & 0.110 & 0.008 & 0.048 & -0.005 & 0.054 \\
$n_3$   & 0.001 & 0.059 & 0.005 & 0.057 & -0.033 & 0.086 & 0.011 & 0.040 & 0.004 & 0.043 \\\hline\hline
\multicolumn{11}{l}{Note: $n_{1}=5000, n_{2}=10000, n_{3}=20000$.}%
\end{tabular}%
\end{adjustbox}
\end{table}

\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)\footnote{%
Information of one's enrollment in PHC is only available from 2011.} of the
Household, Income and Labor Dynamics in Australia (HILDA) Survey data. The
HILDA Survey follows the lives of more than 17,000 Australians each year. As
a rich panel data set, it collects information on many aspects of life in
Australia.\footnote{%
More details about the HILDA data can be found at \href{https://melbourneinstitute.unimelb.edu.au/hilda}%
{https://melbourneinstitute.unimelb.edu.au/hilda}.}

Our empirical model is 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}+\varpi I_{it}\geq \epsilon _{it}\}.  \label{app_model}
\end{equation}%
The dependent variable $y_{it}$ is whether individual $i$ has PHC in year $t$%
. We are interested in the effects of \textquotedblleft Lifetime Health
Cover\textquotedblright\ (LHC) policy, \textquotedblleft Medicare Levy
Surcharge\textquotedblright\ (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$, namlely, $\text{Above30}_{it}:=\mathds{1}%
\{\text{Age}_{it}\geq30\}$. Following the insight of the (sharp)
``regression discontinuity'' design, its coefficient $\beta_1$ 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
now (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.

In addition, we control a (location) dummy variable $\text{GCC}_{it}$
indicating if 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. $I_{it}$ is standardized\footnote{%
Subtract the sample mean from each individual value and then divide the
difference by the standard deviation.} annual household disposable income
using all the sample in the survey, which serves as the continuous regressor
with rich enough support required for point identification. $\epsilon _{it}$ and $\alpha _{i}$ are, respectively, the
usual idiosyncratic error and unobserved heterogeneity in fixed effects
panel data models. Tables \ref{tab_def} and \ref{tab_sum} provide
definitions and summary statistics of all variables, respectively.\footnote{%
Incomes in different years are expressed in the 2011 Australian dollar, and
the standardization is done before dropping missing data.} Note that some
observations are dropped due to missing information in other variables, so the standard deviation of $I_{it}$ reported in Table \ref{tab_sum} 
is not exactly 1. 

\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 \\\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 ,\varpi )$ through maximizing 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 varying 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.

By construction, procedure (\ref{app_Qn}) 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.

We choose $\sigma _{n}=c\cdot \widehat{\text{std}(I_{it})}\sqrt{\log \log n}$ with $%
c=1/2$ and $2/3$ to implement our MS estimation\footnote{%
We conduct the estimation using $\sigma _{n}$ with many different values of $%
c$, ranging from $0.2$ to $1.0$. It turns out that the estimation results
vary little for $c\in \lbrack 0.5,0.7]$. We suspect that the estimation
results may not be reliable for larger $c$ due to too few effective
observations and the dependence between household income and age.} and
report the results in Table \ref{tab_res}. In addition to the estimates of $%
\theta $, we also try calculating the 90\% and 95\% confidence intervals
(CI) for $\theta $ using the $m$-out-of-$n$ bootstrapping. 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 $%
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.302$ (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{6}[2]{*}{$c = 1/2$} & $I_{it}$ & 0.189 & -0.217 & 0.586 & -0.263 & 0.634 \\
          & $y_{it-1}$ & 0.730^{\ast\ast} & 0.691 & 1.554 & 0.651 & 1.636 \\
          & $\textrm{Above30}_{it}$ & 0.655^{\ast\ast} & 0.561 & 1.347 & 0.501 & 1.431 \\
          &  $\textrm{Age}_{it}$ & 0.019 & -0.492 & 0.350 & -0.539 & 0.404 \\
          & $\textrm{GCC}_{it}$   & 0.019 & -0.491 & 0.289 & -0.516 & 0.364 \\
          & $D_{2014,t}$ & 0.039 & -0.458 & 0.357 & -0.498 & 0.419 \\
\hline
    \multirow{6}[2]{*}{$c = 2/3$} & $I_{it}$ & 0.208 & -0.188 & 0.656 & -0.226 & 0.723 \\
          & $y_{it-1}$ & 0.732^{\ast\ast} & 0.697 & 1.503 & 0.630 & 1.624 \\
          & $\textrm{Above30}_{it}$ & 0.624^{\ast\ast} & 0.517 & 1.403 & 0.463 & 1.456 \\
          & $\textrm{Age}_{it}$  & -0.002 & -0.550 & 0.304 & -0.614 & 0.381 \\
          & $\textrm{GCC}_{it}$ & -0.170 & -0.805 & 0.030 & -0.859 & 0.082 \\
          & $D_{2014,t}$ & 0.059 & -0.419 & 0.433 & -0.470 & 0.481 \\
\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. So the following discussion on the empirical
results will be mainly based on the estimates obtained with $c=2/3$. The
estimated results show that household income ($I_{it}$) has a positive but
not significant effect on PHC purchases. However, the existing literature
(e.g., \cite{cheng2014measuring}) has shown a significant positive
correlation between the two. Considering that the confidence interval is
obviously biased to the right of the origin, we doubt that the statistical
insignificance here is due to the large standard error caused by the
relatively small effective sample. There may be a similar situation for the
coefficient on $\text{GCC}_{it}$. The negative (but marginally
insignificant) coefficient indicates that living in GCC may reduce people's
willingness to buy PHC. A possible explanation is that GCC has relatively
sufficient and quality public medical resources, which have a strong
substitution effect on private hospital services. In fact, historical data
show that PHC enrollment experienced a steady decline during the 1990s since
the introduction of Medicare in 1984. The significant positive coefficient
on $y_{it}$ 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 (free-varying). 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
the 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.

\textcolor{purple}{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 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.

