EconBase
← Back to paper

Semiparametric Estimation of Dynamic Binary Choice Panel Data Models

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

89,183 characters · 16 sections · 93 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

SEMIPARAMETRIC ESTIMATION OF DYNAMIC BINARY CHOICE PANEL DATA MODELS

\affil[1]{School of Economics, University of Queensland} \affil[2]{Research School of Economics, Australian National University}

abstractWe propose a new approach to the semiparametric analysis of panel data binary choice models with fixed effects and dynamics (lagged dependent variables). The model under consideration has the same random utility framework as in HonoreKyriazidou2000. We demonstrate that, with additional serial dependence conditions on the process of deterministic utility and tail restrictions on the error distribution, the (point) identification of the model can proceed in two steps, and requires matching only the value of an index function of explanatory variables over time, rather than the value of each explanatory variable. Our identification method motivates an easily implementable, two-step maximum score (2SMS) procedure -- producing estimators whose rates of convergence, in contrast to HonoreKyriazidou2000's (HonoreKyriazidou2000) methods, are independent of the model dimension. We then analyze the asymptotic properties of the 2SMS procedure and propose bootstrap-based distributional approximations for inference. Evidence from Monte Carlo simulations indicates that our procedure performs satisfactorily in finite samples. {{\em JEL classification Codes\/}: C14, C23, C35.} {{\em Keywords\/}: Semiparametric estimation; Binary choice model; Panel data; Fixed effects; Dynamics; Maximum score; Bootstrap}

Introduction

In this paper, we propose a two-step estimation method for panel data binary choice models with fixed effects and dynamics. Specifically, we consider binary choice models of the form

equation[equation omitted — 346 chars of source]

where $T$ is small and $n$ is large, $x_{it}$ is a $K\times1$ vector of (time-varying) explanatory variables,\footnote{Any time-invariant covariates can be thought of as being part of the fixed effect $\alpha_{i}$.} $y_{it-1}$ is the lagged dependent variable, $\alpha_{i}$ represents a time-invariant, individual-specific (fixed) effect, and $\epsilon_{it}$ is an idiosyncratic error term. Both $\alpha_{i}$ and $\epsilon_{it}$ are unobservable to the econometrician. Following HonoreKyriazidou2000 (referred to as HK henceforth), we assume the strong exogeneity that $(x_{i1},...,x_{iT})\perp(\epsilon_{i1},...\epsilon_{iT})|\alpha_i$ and that $\epsilon_{it}$ are independent and identically distributed (i.i.d.) across $t$ conditional on $\alpha_i$. Interest centers on estimating the preference parameter $\theta\equiv(\beta^{\prime},\gamma)^{\prime}$. $y_{i0}$ is assumed to be observed, although the model is not specified in the initial period 0. In the literature, the lagged terms $y_{it-1}$ and the fixed effect $\alpha_{i}$ are referred to as the “state dependence” (see Heckman1981b, Heckman1981a) and the \textquotedblleft unobservable heterogeneity\textquotedblright, respectively. The co-existence of these two terms complicates the identification and estimation of $\theta$, owing to the multiple sources of persistence in $y_{it}$.

This paper resembles other panel data discrete response literature using fixed effects methods, in that we impose no restrictions on the distribution of $\alpha_{i}$, conditional on the observed explanatory variables. ArellanoHonore2001 review early works on estimating $\beta$ in model ((ref)) with no state dependence ($y_{it-1}$). Chamberlain2010 shows that, outside of the logistic case, these static binary choice models have a zero information bound, and the identification requires that at least one of the observed covariates have unbounded support. In the presence of lagged dependent variables, various conditional maximum likelihood methods have been developed for variants of model ((ref)) with logistic errors and at least four observations ($T\geq3$) per individual.\footnote{Throughout this paper, this means the data contain $y_{i0}$ and $(y_{i1},y_{i2},y_{i3},x_{i1},x_{i2},x_{i3})$ for each individual $i$.} honore2021identification provide a comprehensive review of this literature.\footnote{Works such as BartolucciNigro2010, BartolucciNigro2012 and al2017exponential study the estimation of model ((ref)) under alternative specifications.} Several new methods have been proposed for dynamic Logit models based on moment conditions. Leading examples include HonoreWeidner2020, dobronyi-gu-kim, kitazawa2022transformations, and dano2023transition, among others.

HK is the first to consider the semiparametric identification and estimation of model ((ref)). They demonstrate that $\theta$ can be identified if, in addition to assumptions analogous to those in Manski1987, all explanatory variables are strictly exogenous, $\epsilon_{it}$'s are serially independent, and $T\geq3$. However, the rate of convergence of their estimator decreases as the number of continuous regressors increases, and is slower than the standard MS rate derived by KimPollard1990.

There are several alternative fixed effects approaches to the semi- and nonparametric analysis of dynamic binary choice models. HonoreLewbel2002 propose an identification strategy that requires an exclusion restriction (excluded regressor). ChenEtal2019 show that the exclusion restriction in HonoreLewbel2002 implicitly assume (conditional) serial independence of the excluded regressor. Williams2019 studies the nonparametric identification of dynamic binary choice models that satisfy certain exclusion restrictions. In the absence of excluded regressors, some recent works, such as KhanEtal2020 and Aristodemou2021, characterize the (sharp) identified set for $\theta$ under mild conditions. We refer interested readers to the survey article by honore2021identification and Chapter 7 in Hsiaobook for a detailed review of this literature.

This paper takes one step in the direction of HK's semiparametric estimator, in the sense that we provide sufficient conditions under which model ((ref)) can be identified and estimated, without needing to match each of the explanatory variables over time, provided we have at least five observations per individual are observed (i.e., $T\geq4$).\footnote{That is, at least $y_{i0}$ and $(y_{i1},y_{i2},y_{i3},y_{i4},x_{i1},x_{i2},x_{i3},x_{i4})$ are observed for each individual $i$. This is a restriction on the minimum panel length, which is satisfied for many longitudinal panel data sets.} The key insight here is that the identification of $\theta$ can proceed in two steps. First, $\beta$ can be identified based on sequences of $\left\{y_{it}\right\} $, for which $y_{is-1}=y_{it-1}$ and $y_{is+1}=y_{it+1}$ for some $1\leq s< t\leq T-1$ with $t\geq s+2$ (e.g., in the simplest case where $T=4$, $\beta$ is identified based on observations with $y_0=y_2=y_4$), if the distribution of explanatory variables $x_{it}$ satisfies certain serial dependence and stochastic dominance restrictions. Then, using the identified $\beta$, $\gamma$ can be identified by simply matching $x_{it}^{\prime}\beta$ over time. We propose an estimation procedure for $\beta$ and $\gamma$, establish the asymptotics for our estimators, and provide an inference method that uses the bootstrap. We investigate their finite-sample properties using Monte Carlo experiments.

As demonstrated by HonoreTamer2006, matching exogenous utilities over time seems to be essential for the point identification in “distribution-free” dynamic discrete choice models.\footnote{More precisely, HonoreTamer2006 provided examples of point identification failure when it is impossible to match $x_{it}$ over time.} However, the approach proposed here involves matching an identified linear combination of $x_{it}$, rather than HK's matching each component of $x_{it}$. Consequently, in contrast to the results presented by HK, the rates of convergence of our proposed estimators are independent of the dimension of the regressor space, making our approach particularly useful for models with a higher dimensional design.

It is known that panel data binary choice models with unobserved heterogeneity and dynamics can be estimated using the random effects or correlated random coefficients approach. Examples include ArellanoCarrasco2003, Wooldridge2005, and HonoreTamer2006. In addition to preference parameters, these approaches often allow the econometrician to calculate other quantities of interest, such as choice probabilities and marginal effects. However, these approaches require the specification of the statistical relation between the explanatory variables and $\alpha_{i}$. Further, they require one to specify the distribution of $y_{i0}$, conditional on the observed explanatory variables and $\alpha _{i}$, which raises the so-called initial condition problem. Conversely, the fixed effects approaches attempt to estimate preference parameters without making these subtle specifications. Finally, there is also literature exploring the identification and estimation of various partial effects in panel data models; see, for example, AltonjiMatzkin2005, ChernozhukovEtal2013, and more recent advancements by Torgovitsky2019, aguirregabiria2021identification, dobronyi-gu-kim, davezies2021identification, and liu2021identification, among others.

Dynamic binary choice models have a wide range of applications, including the study of labor force participation (DamrongplasitEtal2018), poverty dynamics (Biewen2009), health status (Halliday2008), educational attainment (CameronHeckman1998, CameronHeckman2001), stock market participation (AlessieEtal2004), brand loyalty (ChintaguntaEtal2001), welfare participation (ChayEtal1999), and firm behavior (KerrEtal2014), among others. Most applications typically employ parametric forms of the model ((ref)), such as Logit and Probit, or random effects assumptions. The robustness from the distribution-free and fixed effects specification makes the approach proposed here a competitive alternative to existing parametric and random effects methods. Note that the theoretical validity of our approach relies on certain restrictions on the serial dependence of explanatory variables. Thus, before applying our method, we suggest that applied researchers detrend and seasonally adjust $x_{it}$ in a way that makes them resemble “white noise” conditional on $\alpha_i$.\footnote{Monte Carlo results show that relaxing these restrictions does not significantly affect our estimators' finite-sample performances; see Section (ref) for a more detailed discussion.}

The remainder of this paper is organized as follows. Section (ref) establishes the identification of $\theta$ under different sets of sufficient conditions, based on which, a two-step maximum score (2SMS) procedure is proposed in Section (ref). Sections (ref) and (ref) derive the asymptotic properties of the 2SMS estimator and propose bootstrap-based inference methods, respectively. We present the results of Monte Carlo experiments in Section (ref) that examine the finite-sample performance of our proposed method. Section (ref) concludes the paper. We prove the main theorems and present the main simulation results in the \hyperref[appendix]{Appendices}. The Online Supplement to this paper includes proofs of all technical lemmas, technical details for the bootstrap inference, and results for supplementary simulation studies.

For ease of reference, we next describe the notation maintained throughout this paper.

notationAll vectors are column vectors. $\mathbb{R}^{p}$ is a $p$-dimensional Euclidean space equipped with the Euclidean norm $\Vert\cdot\Vert_{2}$. We reserve the letter $i\in\mathcal{N}\equiv\{1,...,n\}$ for indexing individuals, and the letters $s,t\in\mathcal{T}\equiv\{1,...,T\}$ for indexing time periods. An observation is indexed by $(i,t)$. Vector $x_{its}$ denotes $x_{it}-x_{is}$. The first element of $x_{its}$ is denoted by $x_{its,1}$ and the sub-vector comprising its remaining elements is denoted by $\tilde{x}_{its}$. As is common in the panel data literature, we use the notation $\xi^{t}$ to denote $\left( \xi_{1}^{\prime},...,\xi_{t}^{\prime}\right) ^{\prime}$. For example, suppressing the subscript $i$, $y^{t}\equiv(y_{1},...,y_{t})^{\prime},$ a $t\times1$ vector. $F_{\zeta|\cdot}$ and $f_{\zeta|\cdot}$ denote, respectively, the conditional cumulative distribution function (CDF) and probability density function (PDF) of a random vector $\zeta$ conditional on $\cdot$. For two random vectors, $u$ and $v$, the notation $u\overset{d}{=}v|\cdot$ means that $u$ and $v$ have identical distributions, conditional on $\cdot$, and $u\perp v|\cdot$ means that $u$ and $v$ are independent, conditional on $\cdot$. We use $P(\cdot)$ and $\mathbb{E}[\cdot]$ to denote probability and expectation, respectively. Function $1[\cdot]$ is an indicator function, equal to one when the event in the brackets is true, and zero otherwise. Function $sgn(\cdot)$ denotes the sign function, equal to 1 when $\cdot$ is positive, 0 when $\cdot$ is 0, and $-1$ when $\cdot$ is negative$.$ Symbols $\setminus$, $^{\prime}$, $\propto$, $\Leftrightarrow$, $\overset{d}{\rightarrow}$, and $\overset{P}{\rightarrow}$ represent set difference, matrix transposition, proportionality, “if and only if”, convergence in distribution, and convergence in probability, 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 (bounded in probability), and $a_{n}=o(b_{n})$ ($o_{P}(b_{n})$) means that $a_{n}/b_{n}\rightarrow0$ ($a_{n}/b_{n}\overset{P}{\rightarrow}0$).

Identification

This section provides sufficient conditions for identifying the parameter $\theta$ with no need to match observed covariates $x_{it}$ over time. Under these assumptions, we derive a set of identification inequalities that can be taken to data for (point) estimation and inference on the parameter $\theta$.

To simplify the notation, we suppress the subscript $i$ in the rest of this paper whenever it is clear from the context that all variables relate to each individual. Suppose that a random sample from a population of independent individuals\footnote{Here, the term “independent individuals” refers to the assumption that $(\alpha_i, y_{i0}, x_{i1},...,x_{iT}, \epsilon_{i1},...,\epsilon_{iT})$ is independently distributed across $i$.} is observed for $T+1$ ($=|\mathcal{T}\cup\{0\}|$) periods. Recall that, for all $t\in\mathcal{T}$,

equation[equation omitted — 111 chars of source]

Note that the model is incomplete, in the sense that it does not specify the relationship between $y_{0}$ and $(x^{T},\alpha,\epsilon^{T})$. This is known as the initial condition problem in panel data literature. This paper uses a fixed effects approach, in which we attempt to estimate $\theta = (\beta',\gamma)'$ without making any assumptions on the distribution of $\alpha$, conditional on explanatory variables. This helps us to avoid explicitly specifying the functional form of $p_{0}(x^{T},\alpha)\equiv P(y_{0}=1|x^{T},\alpha)$, and thus circumvents the initial condition problem.

As mentioned, we impose no restriction on $F_{\alpha|x^{T}}$, but place the following restrictions on observed covariates $x^{T}$ and unobserved idiosyncratic errors $\epsilon^{T}$:

assumptionHKFor all $\alpha$ and $s,t\in\mathcal{T}$, \begin{enumerate} • (i) $\epsilon^{T}\perp(x^{T},y_{0})|\alpha$, (ii) $\epsilon _{s}\perp\epsilon_{t}|\alpha$, and (iii) $\epsilon_{s}\overset{d}{=} \epsilon_{t}|\alpha$. • $F_{\epsilon_{t}|\alpha}$ is absolutely continuous with PDF $f_{\epsilon_{t}|\alpha}$ and support $\mathbb{R}$. • (i) One of the regressors, without loss of generality (w.l.o.g.) $x_{ts,1},$ has almost everywhere (a.e.) positive probability density on $\mathbb{R}$, conditional on $\tilde{x}_{ts}$ and $\alpha$, and (ii) the coefficient $\beta_{1}$ on $x_{ts,1}$ is nonzero. • The support $\mathcal{X}_{ts}$ of $F_{x_{ts}|\alpha}$ is not contained in any proper linear subspace of $\mathbb{R}^{K}$. • $\theta=(\beta^{\prime},\gamma)^{\prime}\in\mathcal{B} \times\text{int}(\mathcal{R})$, where $\mathcal{B}\equiv\{b=(b_{1} ,...,b_{K})^{\prime}\in\mathbb{R}^{K}|\Vert b\Vert_{2}=1\}$ and $\mathcal{R}$ is a compact subset of $\mathbb{R}$ with a non-empty interior. \end{enumerate}

Assumption \hyperref[Assumption:HK]{A} places the same set of restrictions on the joint distribution of $(x^{T},\alpha,\epsilon^{T})$ as in HK. While not stated explicitly in their Theorem 4, HK use Assumption \hyperref[Assumption:HK]{A}(a), the exogeneity of $(x^{T},y_{0})$ and serial independence of $\{\epsilon_{t}\}$, conditional on $\alpha$, to derive the moment inequalities for the identification. Note that Assumption \hyperref[Assumption:HK]{A}(a) implies that the fixed effects $\alpha$ pick up two types of dependence in the model: the dependence over time in the unobservables, and the dependence between explanatory variables and unobservables. As a result, in model ((ref)), $\epsilon_{t}$ is independent of $\left( x^{T},y^{t-1}\right) $, conditional on $\alpha$. Furthermore, Assumption \hyperref[Assumption:HK]{A}(a) is a special case of the group homogeneity restriction, $\epsilon_{s}\overset{d}{=}\epsilon_{t}|(x_{s},x_{t},\alpha) $, imposed in Manski1987, PakesPorter2016, and ShiEtal2018 for identifying static discrete choice models (without controlling the lagged term $y_{t-1}$ in the model). Thus, we can suppress the time subscript $t$ in $F_{\epsilon_{t}|\alpha}$ and $f_{\epsilon _{t}|\alpha}$ in the rest of this paper without ambiguity. Assumption \hyperref[Assumption:HK]{A}(b) is a regularity condition that ensures that both $y_{s}\neq y_{t}$ and $y_{s}=y_{t}$ occur with positive probabilities for all $\alpha$ and $s,t\in\mathcal{T}$.

It is known and documented in the relevant literature (see, e.g., Lemma 1 of Manski1985) that to establish the point identification of the parameter $\theta$ in a distribution-free setting, $x_{t}$ also needs to satisfy certain regularity conditions. Assumption \hyperref[Assumption:HK]{A}(c) requires the existence of a relevant, continuous regressor, with large support, which is a standard restriction imposed in MS-type estimators. Assumption \hyperref[Assumption:HK]{A}(d) is the familiar full-rank condition. Assumptions \hyperref[Assumption:HK]{A}(c) and \hyperref[Assumption:HK]{A}(d) are identical to Assumption 2 of Manski1987.

Assumption \hyperref[Assumption:HK]{A}(e) is for scale normalization and parameter space. This is a typical practice for discrete choice models, because the identification of $\theta$ is only up to scale. In the semiparametric framework, where no parametric form of $F_{\epsilon|\alpha}$ is specified, identification is often achieved by normalizing the magnitude of the regression coefficients. Assumption \hyperref[Assumption:HK]{A}(e) assumes that $\beta$ is on the unit circle and has a nonzero first element $\beta_{1}$.\footnote{Our procedure identifies $\beta$ and $\gamma$ sequentially, so it is more convenient to normalize the scale of $\beta$, rather than that of $\theta$, as in HK.}

HK demonstrate that, if $T\geq3$ and $x_{t}$ has time-varying overlap support, $\theta$ can be identified under Assumption \hyperref[Assumption:HK]{A}.\footnote{As is stated in HK, Assumption \hyperref[Assumption:HK]{A} is not sufficient for point identifying $\theta$ if $T<3$.} Their proposed approach requires matching all exogenous covariates over time, and results in an estimator with a rate that declines as the number of exogenous covariates increases. The main contribution of this paper is the provision of a set of supplementary conditions, under which the identification of $\theta$ can escape from the necessity of element-by-element matching. Specifically, our approach is based on the following monotonic relationship between a conditional choice probability and an index of the exogenous covariates: For some $s,t\in \mathcal{T}$ such that $t-s\geq2$,

align[align omitted — 255 chars of source]

Note that ((ref)) requires that there be at least five ($T\geq4$) observations per individual observed by the econometrician (i.e., $s=1$, $t=3$, and $s+1=t-1=2$). In the simplest case with $T=4$, ((ref)) reduces to $P(y_{3}=1|x_{1},x_{3},y_{0}=y_{2}=y_{4},\alpha) \geq P(y_{1}=1|x_{1},x_{3},y_{0}=y_{2}=y_{4},\alpha)$ if and only if $x_{3}^{\prime}\beta \geq x_{1}^{\prime}\beta$.

Result ((ref)) states that the indices $x_{s}'\beta$ and $x_{t}'\beta$ rank order the (conditional) probabilities of choosing 1 in periods $s$ and $t$. To ensure this, conditioning on $y_{s-1}=y_{t-1}$ is obviously necessary. However, as $t-1>s$, $y_{s}$ affects $y_{t-1}$ through the dynamics of the model (specifically, through the chain $y_{s}\rightarrow y_{s+1}\rightarrow\cdots\rightarrow y_{t-1}$). Accordingly, we need to include $y_{s+1}$ in the conditioning set to cut off such state dependence, and further impose the restriction $y_{s+1}=y_{t+1}$ to make the events $\{y_{s}=1\}$ and $\{y_{t}=1\}$ have symmetric conditioning sets. Only with this symmetry can we invoke time stationarity restrictions on $x_{t}$ to establish the equivalence in ((ref)).

In particular, to reach ((ref)), we also need to address the following two concerns. First, $x_{s}$ ($x_{t}$) may affect the value of $y_{t+1}$ ($y_{s+1}$) via its serial dependence on $x_{t+1}$ ($x_{s+1}$). Second, the dependence between $x_{t}$ and $y_{t+1}$ (via $x_{t+1}$) may change dramatically over time. Both require additional restrictions to be placed on the serial dependence of the stochastic process of $x_{t}$. Otherwise, as shown in Appendix (ref), $x_t'\beta$ is not the unique factor that can rank order the choice probabilities in ((ref)).

The following condition, together with Assumption \hyperref[Assumption:HK]{A}, is sufficient to establish ((ref)), as shown in Appendix (ref).

assumptionSIFor all $s,t\in\mathcal{T}$ such that $s\neq t$, (a) $x_{s}\perp x_{t}| \alpha$, and (b) $x_{s}\overset{d}{=}x_{t}|\alpha$.

In Appendix (ref), we first prove ((ref)) under Assumptions \hyperref[Assumption:HK]{A} and \hyperref[Assumption:SI]{SI} for a special case of model ((ref)) with $T=4$ and $\gamma<0$. This serves as a roadmap to help readers understand the main ideas. The same arguments can be applied analogously to prove the most general case; see Lemma (ref).

Assumption \hyperref[Assumption:SI]{SI} imposes a strong restriction on the dynamic process of the covariate sequence, which requires the process $\{x_{t}\}$ to be serially independent and strictly stationary, conditional on the individual-specific effects $\alpha$. In a dynamic fixed effects model, $\alpha$ collects all time-invariant covariates, as well as unobserved individual preferences, abilities, or character traits. In such models, if $x_{t}$ includes only observed individual characteristics naturally correlated with $\alpha$, it may be reasonable to further assume that the serial dependence in the process $\{x_{t} \}$ is also derived from $\alpha$. Assumption \hyperref[Assumption:SI]{SI} implies that we cannot accommodate time trends. We do allow random time effects $\lambda_t$ that satisfy $(\lambda_{s},x_{s})\perp(\lambda _{t},x_{t})|\alpha$, and $(\lambda_{s},x_{s})\overset{d}{=}(\lambda_{t} ,x_{t})|\alpha$.\footnote{We can apply essentially the same arguments used in Appendix (ref) to establish a monotonic relationship analogous to ((ref)) for index $\lambda+x_t'\beta $, and identify $\beta$ and $\gamma$ similarly.}

If $x_{t}$ contains covariates related to some institutional factors that lead to exogenous variation in, for example, costs of participation, across individuals, Assumption \hyperref[Assumption:SI]{SI} may be approximately satisfied by using the differencing, demeaning, or de-trending transformation of these variables. This applies to cases where $\{x_{t}\}$ exhibits some long-run equilibrium (trend, deterministic or stochastic). The transformed regressor then measures the deviation of $x_{t}$ from its long-run equilibrium (trend), which, in some cases, might be assumed to be a white noise process affecting the short-run dynamics of the model.\footnote{For example, consider a case where $\{x_{t}\}$ is a random-walk-plus-drift process (i.e., $x_{t}=x_{0}+a_{0}t+\sum_{\tau=1}^{t}e_{\tau}$). Although $\{x_{t}\}$ violates Assumption SI, its first differencing $\Delta x_{t}=x_{t}-x_{t-1}=a_{0}+e_{t}$ is i.i.d. over time.} Note that such variable transformation may involve model reparameterization. For example, when model ((ref)) is the reduced form derived from some structural model, one would need to first reparameterize the structural model accordingly, so that the estimates of the coefficient vector in model ((ref)) with transformed covariates can be interpreted in a meaningful way.

HK assume that the support of $x_{t}$ is overlapping over time, so the differences in the regressors across different time periods have a positive density in a neighborhood of zero. However, evidence presented in HonoreTamer2006 implies that some additional assumption is needed to achieve point identification without performing an element-by-element match, as in HK. Indeed, Assumption \hyperref[Assumption:SI]{SI} is the extra condition needed for our approach, compared with the semiparametric estimator in HK.

Under Assumptions \hyperref[Assumption:HK]{A} and \hyperref[Assumption:SI]{SI}, the identification of $\theta$ proceeds in two steps. Proposition (ref) demonstrates that $\beta$ can be identified based on moment inequality ((ref)), and Proposition (ref) establishes the identification of $\gamma$ by matching the value of the index function $x_{t}^{\prime}\beta$ in different periods.

proposition[Identification of $\beta$] Assume $T \geq 4$. For all $s,t\in\mathcal{T}$ such that $t\geq s+2$, define \begin{align*} Q_{1}(b) = & \mathbb{E}\left\{ [P(y_{t}=1|x_{s},x_{t},y_{s-1} =y_{t-1},y_{s+1}=y_{t+1})- P(y_{s}=1|x_{s},x_{t},y_{s-1}=y_{t-1} ,y_{s+1}=y_{t+1})]\right. \\ & \left. \times sgn(x_{ts}^{\prime}b)\right\} . \end{align*} Suppose Assumptions \hyperref[Assumption:HK]{A} and \hyperref[Assumption:SI]{SI} hold. Then, $Q_{1}(\beta)>Q_{1}(b)$, for all $b\in\mathcal{B} \setminus\{\beta\}$.

The proof of Proposition (ref) can be found in Appendix (ref). Note that our identification strategy for $\beta$ requires $T=4$, as a minimum. In this case, $t=s+2$ must hold with $s=1$ and $t=3$, and thus $Q_{1}(b)$ is

equation[equation omitted — 173 chars of source]

Proposition (ref) establishes the identification of $\beta$, which enables us to identify $\gamma$, with $\beta$ being treated as a known, constant vector. Then, the following proposition shows that $\gamma$ can be identified by matching the deterministic utility $w_{t}\equiv x_{t}^{\prime }\beta$ in different periods; the proof is presented in Appendix (ref). Because the key idea for identifying $\gamma$ uses the insight of Section 4.1 in HK, in what follows, we keep the notation as close to that of HK as possible.

We define the event \[ A=\{y_{0}=d_{0},...,y_{s-1}=d_{s-1},y_{s}=0,y_{s+1}=d_{s+1},...,y_{t-1} =d_{t-1},y_{t}=1,y_{t+1}=d_{t+1},...,y_{T}=d_{T}\}, \] and its counterpart \[ B=\{y_{0}=d_{0},...,y_{s-1}=d_{s-1},y_{s}=1,y_{s+1}=d_{s+1},...,y_{t-1} =d_{t-1},y_{t}=0,y_{t+1}=d_{t+1},...,y_{T}=d_{T}\}, \] where $d_{\tau}\in\{0,1\}$, for $0\leq\tau\leq T$. Note that $y$ takes the same values other than at time periods $(s,t)$ for $A$ and $B$: $y$ switches from $0$ to $1$ at time $s$ and $t$, respectively, for $A$, and $y$ switches from $1$ to $0$ at time $s$ and $t$, respectively, for $B$.

We have two cases, based on whether $s$ and $t$ are adjacent. When $s$ and $t$ are adjacent ($t=s+1$), we define the objective function \[ Q_{2}(r;\beta)=\mathbb{E}\left\{ \left[ P(A|x^{T},w_{t}=w_{t+1} )-P(B|x^{T},w_{t}=w_{t+1})\right] sgn\left( (w_{t}-w_{t-1})+r(d_{t+1} -d_{t-2})\right) \right\} . \] For the case where $s$ and $t$ are not adjacent ($t>s+1$), we define the objective function

align*[align* omitted — 254 chars of source]

In the following proposition, we establish the identification of $\gamma$ by showing that $\gamma$ uniquely maximizes both $Q_{2}(r;\beta)$ and $\tilde{Q}_{2}(r;\beta)$.

From the definitions of $A,B,$ and $Q_{2}(r;\beta),$ we require $T\geq3$. In the simplest case when $T=3$, we have \[ A=\{y_{0}=d_{0},y_{1}=0,y_{2}=1,y_{3}=d_{3}\}\text{ and }B=\{y_{0}=d_{0} ,y_{1}=1,y_{2}=0,y_{3}=d_{3}\}, \] and \[ Q_{2}(r;\beta)=\mathbb{E}\left\{ \left[ P(A|x^{T},w_{2}=w_{3})-P(B|x^{T} ,w_{2}=w_{3})\right] sgn\left( (w_{2}-w_{1})+r(d_{3}-d_{0})\right) \right\} . \] $\tilde{Q}_{2}(r;\beta)$ is not applicable for this case.

proposition[Identification of $\gamma$] Suppose Assumption \hyperref[Assumption:HK]{A} holds. We have \begin{enumerate} • $Q_{2}(\gamma;\beta)>Q_{2}(r;\beta)$ for all $r\in\mathcal{R} \setminus\{\gamma\}$, and • $\tilde{Q}_{2}(\gamma;\beta)>\tilde{Q}_{2}(r;\beta)$ for all $r\in\mathcal{R}\setminus\{\gamma\}$. \end{enumerate}
remarkWhen $T\geq4$, any combination $(s,t)$ of the elements of $\{1,...,T-1\}$ taken two at a time can be used to construct the population objective function to identify $\gamma$. For example, in the simplest case $T=4$, feasible choices of $(s,t)$ include $(1,2)$, $(1,3)$, and $(2,3)$. One can use any of these pairs to define the population objective function, either $Q_{2} (\cdot;\beta)$ or $\tilde{Q}_{2}(\cdot;\beta)$. Clearly, any one (or combination, e.g., by simply summing them up) of these objective functions can be used to identify $\gamma$.

Propositions (ref) and (ref) outline a two-step procedure for identifying the preference parameters $\beta$ and $\gamma$, of which Proposition (ref) uses HK's insight. Note that, as Proposition (ref) suggests, an additional assumption, \hyperref[Assumption:SI]{SI}, enables us to establish the identification of $\beta$ independently to that of $\gamma$ in the first step. As a result, it suffices to match the index $x_{t}^{\prime }\beta$, rather than each component of $x_{t}$ over time, as in HK, when identifying $\gamma$ in the second step. The benefit of doing so is that the two-step procedure avoids the curse of dimensionality caused by matching many explanatory variables (see Theorem (ref) in Section (ref)). Our method is particularly competitive when handling high-dimensional models.

The following theorem is an immediate result of Propositions (ref) and (ref).

theorem[Identification of $\theta$] Suppose $T\geq 4$ and Assumptions \hyperref[Assumption:HK]{A} and \hyperref[Assumption:SI]{SI} hold. Then, $\beta$ is identified based on population objective function $Q_{1}(\cdot)$, and $\gamma$ is identified based on either population objective function $Q_{2}(\cdot;\beta)$ or $\tilde{Q}_{2}(\cdot;\beta)$.

Alternative Sufficient Conditions for Identification

Here, we provide an alternative sufficient condition that permits limited dependence of the covariates for the identification. We show in Lemma (ref) that Assumptions \hyperref[Assumption:HK]{A} and \hyperref[Assumption:SD]{SD} (in below) are sufficient for the inequality in ((ref)), which, in turn, shows the identification. In this section, we present and discuss this assumption.

assumptionSDFor all $\alpha$ and $s,t\in\mathcal{T}$, \begin{enumerate} • $f_{\epsilon|\alpha}(\cdot)/F_{\epsilon|\alpha}(\cdot)$ is a non-increasing function, or equivalently, $f_{\epsilon|\alpha}(\cdot)/[1-F_{\epsilon|\alpha}(\cdot)]$ is a non-decreasing function. • Let $w_t\equiv x_t'\beta$. The joint PDF of $w^{T}$ conditional on $\alpha$ is exchangeable, i.e., \[ f_{w^{T}|\alpha}(\omega_{1},...,\omega_{T})=f_{w^{T}|\alpha}(\omega_{\pi(1)},...,\omega_{\pi(T)}) \] for all permutations $\{\pi(1),...,\pi(T)\}$ defined on the set $\mathcal{T}$. \end{enumerate}

Assumption \hyperref[Assumption:SD]{SD}(a) states that $F_{\epsilon|\alpha}$ has a decreasing inverse Mills ratio, which, together with Assumption \hyperref[Assumption:SD]{SD}(b), guarantees the monotonic relation in ((ref)), as proved in Appendix (ref). Assumption \hyperref[Assumption:SD]{SD}(a) is satisfied by many common continuous distributions, such as the Gaussian, logistic, Laplace, uniform, gamma, log-normal, Gumbel, and Weibull distributions.\footnote{In a mixture model, e.g., \[ f_{\epsilon|\alpha}(e)=\sum_{m=1}^{M}\pi_{m}f_{\epsilon|\alpha}(e;\vartheta _{m}) \] with mixing proportions $\pi_{m}$, $\sum_{m=1}^{M}\pi_{m}=1$, where each component density has a different parameter vector $\vartheta_{m}$, Assumption \hyperref[Assumption:SD]{SD}(a) holds for $F_{\epsilon_{t}|\alpha}(\cdot)$ if it is satisfied by all component distributions $F_{\epsilon|\alpha} (\cdot;\vartheta_{m})$.} However, this property fails if $F_{\epsilon |\alpha}$ has heavy tails (e.g., Student's $t$-distribution and Cauchy distribution).\footnote{More precisely, Assumption \hyperref[Assumption:SD] {SD}(a) does not hold globally for these distributions. For example, it is not difficult to find that this assumption holds for Student's $t$ and Cauchy on $\lbrack-L,\infty)$ for some positive $L$.} Note that Assumption \hyperref[Assumption:SD]{SD}(a) is a key condition imposed in McFadden1976 and Silvapulle1981 for both $-\log F_{\epsilon|\alpha}(\cdot)$ and $-\log(1-F_{\epsilon|\alpha}(\cdot) )$ being convex, which guarantees a unique solution for the MLE in cross-sectional models with errors that follow a general distribution.

In model ((ref)), the exogenous utility $w_{t}$ affects the value of $y_{t+1}$ via $y_{t}$ and its serial dependence with $w_{t+1}$, conditioning on $\alpha$. The former is explicitly captured by the coefficient $\gamma$. For the latter, Assumption \hyperref[Assumption:SD]{SD}(b) restricts the serial dependence of $\{w_{t}\}$. Assumption \hyperref[Assumption:SD]{SD}(b) is weaker than Assumption \hyperref[Assumption:SI]{SI}. Under Assumption \hyperref[Assumption:SI]{SI}, $w_{t}$'s are i.i.d. over time, conditional on $\alpha$. Then, $w_{t}$'s have an exchangeable joint PDF, as defined in Assumption \hyperref[Assumption:SD]{SD}(b). However, the other direction is not always true. As noted in Fox2007, a common example of an exchangeable PDF for non-independent $w_{t}$'s is a multivariate normal density with $\mathbb{E}[w_{t}]=\mu$, $\textrm{Var}(w_{t})=\sigma^{2}$, and $\textrm{Corr}(w_{s},w_{t})=\rho$, for all $s,t\in\{1,...,T\}$. More generally, if each $w_{t}$ can be expressed as $w_{t}=\varphi(u_{t},v)$, where $u_t$'s are i.i.d. random variables, $v$ is a random variable independent of all $u_t$'s, and $\varphi(\cdot,\cdot)$ is some measurable function, then $w^{T}$ satisfies Assumption \hyperref[Assumption:SD]{SD}(b), but not Assumption \hyperref[Assumption:SI]{SI}. Similar exchangeability assumptions are imposed in AltonjiMatzkin2005 and ChenEtal2018.

A few remarks are in order about how our identification conditions are related to the existing literature.

remarkCompared with HK, our approach relies on additional assumptions restricting the serial dependence of strictly exogenous regressors $x_{t}$ and requires $T\geq4$. These conditions make identification without element-by-element matching of $x_{t}$ possible. HK construct identifying inequalities similar to our ((ref)). To obtain point identification, HK use probabilities of specific sequences of $y^{T}$ conditional on event $\{x_{s}=x_{t}\}$, for some $s,t\in\{1,...,T\}$. Instead, our approach matches $y_{t}$ in different time periods to construct identifying inequalities, allowing for point identification without element-by-element matching of $x_{t}$. However, as discussed after the introduction of our key identifying condition ((ref)), $x_{t}$ can affect the choice probabilities in ((ref)) through the utility index $x_{t}'\beta$ and its serial dependence with $x$ in other time periods. Assumptions \hyperref[Assumption:SI]{SI} and \hyperref[Assumption:SD]{SD} restrict this dependence, ensuring that choice probabilities are solely rank-ordered by $x_{t}'\beta$. Additionally, our approach requires comparing $x_{t}'\beta$ in non-adjacent time periods (i.e., $t>s+1$ in ((ref))). As a result, we need observation of one more period for each individual, compared with HK's method, which can achieve point identification using $x_{t}'\beta$ in adjacent time periods ($t=s+1$).
remarkSecond, our identification conditions are non-nested with those in the literature, assuming exclusion restrictions, such as HonoreLewbel2002, ChenEtal2018, ChenEtal2019, and Williams2019. ChenEtal2019 show that HonoreLewbel2002 essentially require the serial independence of the excluded regressor. Williams2019 requires that the other strictly exogenous regressors are conditionally independent of the past values of the excluded regressor. In addition to specific restrictions on the dynamic process for the covariates, the identification results of these studies rely on the existence of at least one \textquotedblleft excluded regressor\textquotedblright \ conditionally independent of the individual fixed effects $\alpha$. Conversely, our approach allows for arbitrary correlation between $x_{t}$ and $\alpha$.

Estimation

Applying the analogy principle, the identification results presented in Section (ref) can be translated into a two-step estimation procedure. In the first step, we obtain an MS estimator (with binary weights) $\hat{\beta}$ of $\beta$. In the second step, $\gamma$ is estimated by a localized MS procedure matching the estimated index $x_{t}^{\prime}\hat{\beta}$ over time. Each of the two steps is described, in turn, below.

In Sections (ref) and (ref), we restrict our discussion to the model with $T=4$ to streamline exposition in subsequent sections. The same method can be applied, with straightforward modification, to models with longer panels. We provide objective functions for general cases with $T\geq4$ in Section (ref).

Estimation of \texorpdfstring{$\beta$}{Lg} with \texorpdfstring{$T=4$}{Lg}

Assuming a random sample of $n$ individuals, we propose the following weighted MS estimator $\hat{\beta}$ of $\beta$, defined as the maximizer over the parameter space $\mathcal{B}$:

equation[equation omitted — 79 chars of source]

where

equation[equation omitted — 133 chars of source]

Because we restrict the search within a compact set $\mathcal{B}$ and the objective function ((ref)) is bounded and continuous, the maximizer $\hat{\beta}$ of the maximization problem ((ref)) does exist. However, $\hat{\beta}$ may not be unique because the objective function ((ref)) is essentially a step function for any finite samples. Nonetheless, as guaranteed by Theorem (ref) in Section (ref), $\hat{\beta}$ is in a small neighborhood of the true parameter $\beta$ with a high probability for a sufficiently large sample size.

It is clear from expression ((ref)) that only observations that satisfy $y_{i1}\neq y_{i3}$, $y_{i0}=y_{i2}$, and $y_{i2}=y_{i4}$ are used in the estimation. That is, the objective function uses only “switchers” with choice changes in periods $1$ and $3$, with the same choices in their previous and subsequent periods, respectively. This feature reduces the “effective” sample size for the estimator of $\beta$. HK's estimator has a similar problem, because it also uses only switchers, and needs to match $x_t$ over time. Our estimator ((ref)) is more applicable when the model has many regressors, especially discrete regressors that must be matched exactly over time when applying HK's procedure.

Estimation of \texorpdfstring{$\gamma$}{Lg} with \texorpdfstring{$T=4$}{Lg}

Proposition (ref) motivates a localized MS estimator $\hat{\gamma }$ of $\gamma$, defined here as the maximizer over the parameter space $\mathcal{R}$ of the objective function\footnote{If one knew $\beta$, the estimation of $\gamma$ requires only $T=3$, that is, using the first line of ((ref)) as the objective function.}

align[align omitted — 341 chars of source]

Expression ((ref)) is the sample analogue of $Q_{2}(r;\beta)$ in Proposition (ref) after taking the union of events $A$ and $B$ for all possible values of $d_0,d_{1},...,d_{4}$. As with objective function ((ref)), ((ref)) also uses only data on switchers (i.e., satisfying $A\cup B$) who make different choices in the two periods compared. In addition, ((ref)) also requires a match in $x_{t}^{\prime}\beta$.

Note that this estimator is not feasible, because $\beta$ is unknown, and it is of probability zero to have exactly matched indices ($x_{is}^{\prime} \beta=x_{it}^{\prime}\beta$) in the presence of continuous regressors. To resolve the first concern, we propose replacing the unknown parameter $\beta$ in expression ((ref)) with the $\hat{\beta}$ obtained from ((ref)), which is shown to be (cube-root $n$) consistent in Section (ref).

For the second concern, we use kernel weights \[ \mathcal{K}_{h_{n}}((x_{it}-x_{is})^{\prime}b),\text{ for all }s,t\in \mathcal{T}\text{ and }b\in\mathcal{B}, \] instead of $1[x_{is}^{\prime}b=x_{it}^{\prime}b].$ $\mathcal{K}_{h_{n}} (\cdot)$ is defined as $h_{n}^{-1}\mathcal{K}(\cdot/h_{n})$, where $\mathcal{K}(\cdot)$ is a kernel density function and $h_{n}$ is a bandwidth sequence that converges to zero as $n\rightarrow\infty$. The idea is to replace the binary weights for $x_{is}^{\prime}\hat{\beta}=x_{it}^{\prime}\hat{\beta}$ with weights that depend inversely on the magnitude of $(x_{it}-x_{is})^{\prime}\hat{\beta}$, giving more weight to observations with $(x_{it}-x_{is})^{\prime}\hat{\beta}$ closer to zero. We discuss the choice of the tuning parameter $h_n$ with illustrating examples in Section (ref).

Then, we propose the following kernel weighted MS estimator $\hat{\gamma}$ of $\gamma$:

equation[equation omitted — 96 chars of source]

where

align[align omitted — 372 chars of source]
remarkNote that objective function ((ref)) is associated with the population objective function $Q_{2}(r;\beta)$ in Proposition (ref), which uses only observations of adjacent time periods. Applying the same idea to {the} population objective function $\tilde{Q}_{2}(r;\beta)$ yields the following objective function, using observations which are not adjacent: \[ \tilde{Q}_{2n}^{K}(r;\hat{\beta})=\frac{1}{n}\sum_{i=1}^{n}1[y_{i2} =y_{i4}]\mathcal{K}_{h_{n}}(x_{i42}^{\prime}\hat{\beta})(y_{i3}-y_{i1})\cdot sgn(x_{i31}^{\prime}\hat{\beta}+r(y_{i2}-y_{i0})). \] In practice, to make full use of all observations, one can consider using $Q_{2n}^{K}(r;\hat{\beta})+\tilde{Q}_{2n}^{K}(r;\hat{\beta})$ as {the} objective function for the estimation of $\gamma$.
remarkCalculating the MS-type of estimators ((ref)) and ((ref)) is challenging, as it is for the semiparametric estimator of HK. Following Fox2007 and YanYoo2019, we suggest using a global optimization method called the differential evolution (DE) algorithm. The DE algorithm, introduced by StornPrice1997, is specifically designed to search the global optimum of a real-valued function with real-valued parameters. Notably, it does not require the objective function to be continuous or differentiable. The DE algorithm has been widely used in engineering applications, and its performance as a global optimization algorithm has been studied extensively (See, e.g., StornPrice2006). To implement the DE algorithm, one can use the “DEoptim” package in R.\footnote{\url{https://cran.r-project.org/web/packages/DEoptim/index.html}.} The following statement is quoted from the R documentation for the “DEoptim” package, which provides a brief introduction. Interested readers are advised to consult this documentation for more information on the implementation and usage of this algorithm. \textquotedblleft Differential Evolution (DE) is a search heuristic introduced by StornPrice1997. Its remarkable performance as a global optimization algorithm on continuous numerical minimization problems has been extensively explored; see StornPrice2006. DE belongs to the class of genetic algorithms which use biology-inspired operations of crossover, mutation, and selection on a population in order to minimize an objective function over the course of successive generations. As with other evolutionary algorithms, DE solves optimization problems by evolving a population of candidate solutions using alteration and selection operators. DE uses floating-point instead of bit-string encoding of population members, and arithmetic operations instead of logical operations in mutation. DE is particularly well-suited to find the global optimum of a real-valued function of real-valued parameters, and does not require that the function be either continuous or differentiable.\textquotedblright

Note that the 2SMS procedure described in ((ref))--((ref) ) and ((ref))--((ref)) does not require matching each covariate in $x_{it}$ over time, as it does in HK. As a result, the rates of convergence of $\hat{\beta}$ and $\hat{\gamma}$ are independent of the number of continuous covariates in $x_{it}$, in contrast to the procedure of HK. In view of existing results on the MS estimators (e.g., Manski1985, Manski1987, KimPollard1990, and SeoOtsu2018), we expect the limiting distributions of $\hat{\beta}$ and $\hat{\gamma}$ to be non-Gaussian and their rates of convergence to be $O_{P}(n^{-1/3})$ and $O_{P}((nh_{n})^{-1/3})$, respectively. Section (ref) states sufficient conditions under which these asymptotic properties can be derived.

Estimation with \texorpdfstring{$T\geq4$}{Lg}

A longer panel allows for more objective functions of similar form. Collectively, these objective functions (by, for example, summing them) can be used to obtain more accurate estimates of $\theta$ for finite samples. For the case with $T\geq 4$, estimators for $\beta$ and $\gamma$ that best use the data can be obtained as follows. For $\beta,$ we find $\hat{\beta}$ by maximizing \[ Q_{1n}(b)=\frac{1}{n}\sum_{i=1}^{n}\sum_{t>s+1}1[y_{is-1}=y_{it-1} ]1[y_{is+1}=y_{it+1}](y_{it}-y_{is})sgn\left( (x_{it}-x_{is})^{\prime }b\right) . \] Once $\hat{\beta}$ is obtained, we estimate $\gamma$ by maximizing \[ Q_{2n}^{K}(r;\hat{\beta})+\tilde{Q}_{2n}^{K}(r;\hat{\beta}) \] with respect to $r$, where \[ Q_{2n}^{K}(r;\hat{\beta})=\frac{1}{n}\sum_{i=1}^{n}\sum_{t=2}^{T-1} \mathcal{K}_{h_{n}}((x_{it+1}-x_{it})^{\prime}\hat{\beta}) (y_{it} -y_{it-1})sgn((x_{it}-x_{it-1})^{\prime}\hat{\beta}+r(y_{it+1}-y_{it-2})) \] is for the case with $t=s+1$, and

align*[align* omitted — 347 chars of source]

is for the case with $t>s+1$.

Asymptotic Properties

The estimators proposed in Section (ref) are of the same structure and differ only in that they each use a different fraction of observations in the sample. We expect that they have similar asymptotic properties. Therefore, it suffices to show the asymptotics for the estimators in Sections (ref) and (ref), for the case $T=4$. The asymptotic properties of the estimators in Section (ref) can be derived in a similar way.

As is standard in the literature, such as KimPollard1990, we start the analysis by introducing modified objective functions for $\hat{\beta}$ and $\hat{\gamma}$. The new objective functions are monotone (linear) transformations of ((ref)) and ((ref)), respectively. As a result, working with them does not change the values of $\hat{\beta}$ and $\hat{\gamma}$, but can facilitate the derivation process.

Because adding terms not related to $b$ will not affect the optimization over $b,$ and $1[a>0]=(sgn(a)+1)/2$ for all $a\in\mathbb{R}$, $\hat{\beta}$ obtained from the following objective function is identical to that from ((ref)), \[ \hat{\beta}=\arg\max_{b\in\mathcal{B}}n^{-1}\sum_{i=1}^{n}\xi_{i}\left( b\right) , \] where

equation[equation omitted — 227 chars of source]

For the same reason,\ $\hat{\gamma}$ can be obtained equivalently from \[ \hat{\gamma}=\arg\max_{r\in\mathcal{R}}n^{-1}\sum_{i=1}^{n}\varsigma _{ni}( r,\hat{\beta}) , \] where

align[align omitted — 437 chars of source]

The following technical assumptions are needed for the asymptotics of $\hat{\beta}$ and $\hat{\gamma}$.

assumptionThe vectors $\left( x_{i}^{T},y_{i}^{T},y_{i0}\right) ^{\prime}$ with $T\geq4$ are i.i.d. across individuals.
assumption$n^{-1}\sum_{i=1}^{n}\xi_{i}( \hat{\beta}) \geq\max _{b\in\mathcal{B}}n^{-1}\sum_{i=1}^{n}\xi_{i}\left( b\right) -o_{P}\left( n^{-2/3}\right) $ and $n^{-1}\sum_{i=1}^{n}\varsigma_{ni}( \hat{\gamma} ,\hat{\beta}) \geq\max_{r\in\mathcal{R}}n^{-1}\sum_{i=1}^{n}\varsigma_{ni}( r,\hat{\beta}) -o_{P}\left( ( nh_{n}\right) ^{-2/3})$.
assumptionThe joint density function for $\alpha$, covariates $x^{T}$, and $\epsilon^{T}$ are continuously differentiable. The density function and its first-order derivatives are uniformly bounded. Further, \[ f\left( \epsilon_{t}|\alpha,x^{T},\epsilon_{1},...,\epsilon_{t-1},\epsilon _{t+1},...,\epsilon_{T}\right) \text{and } f\left( x_{t}|\alpha,x_{1} ,...,x_{t-1},x_{t+1},...,x_{T},\epsilon^{T}\right) \] are continuous differentiable with respect to all arguments. The conditional densities and their derivatives are uniformly bounded.
assumptionThe kernel function $\mathcal{K}\left( u\right) $ is nonnegative, symmetric about zero, continuous differentiable, has compact support, and satisfies $\int_{\mathbb{R}}\mathcal{K}\left( u\right) du=1$.
assumption$h_{n}\rightarrow0,$ $nh_{n}\rightarrow\infty,$ and $nh_{n} ^{4}\rightarrow0$ as $n\rightarrow\infty.$

Assumption (ref) is standard in the literature and precisely defines our estimator. Assumption (ref) is for technical convenience; it ensures certain functions defined in the proof of Theorem (ref) are differentiable, so that $V_1$ and $V_2$ (defined in Theorem (ref)) have simple representations. In the case of discrete explanatory variables which violates Assumption (ref), $V_1$ and $V_2$ can be shown to be well defined, but with more tedious calculations and notation. Assumption (ref) collects some standard restrictions on kernel functions. The symmetry of $\mathcal{K}\left( u\right)$ ensures that the bias term is of the order $h_{n}^{2}$. In Assumption (ref), $nh_{n}\rightarrow\infty$ is standard, and $nh_{n}^{4}\rightarrow0$ ensures the bias term from the kernel estimation is asymptotically negligible.

theoremSuppose $T\geq 4$ and Assumptions \hyperref[Assumption:HK]{A}, \hyperref[Assumption:SI]{SI} (or \hyperref[Assumption:SD]{SD}), and (ref)--(ref)\ hold. Then, \begin{enumerate} • $\hat{\beta}-\beta=O_{P}\left( n^{-1/3}\right) ,$ and \[ n^{1/3}( \hat{\beta}-\beta) \overset{d}{\rightarrow}\arg \max_{\boldsymbol{s}\in\mathbb{R}^{K}}Z_{1}\left( \boldsymbol{s}\right) , \] where $Z_{1}\left( \boldsymbol{s}\right) $ is a Gaussian process with continuous sample paths, expected value $\frac{1}{2}\boldsymbol{s}^{\prime }V_{1}\boldsymbol{s}$, and covariance kernel $H_{1}\left( \boldsymbol{s} ,\boldsymbol{t}\right) .$ $V_{1}$ and $H_{1}$ are defined in expressions ((ref)) and ((ref)), respectively$.$$\hat{\gamma}-\gamma=O_{P}( \left( nh_{n}\right) ^{-1/3}) ,$ and \[ \left( nh_{n}\right) ^{1/3}\left( \hat{\gamma}-\gamma\right) \overset{d}{\rightarrow}\arg\max_{s\in\mathbb{R}}Z_{2}\left( s\right) . \] where $Z_{2}\left( s\right) $ is a Gaussian process with continuous path, expected value $\frac{1}{2}V_{2}s^{2}$, and covariance kernel $H_{2}\left( s,t\right) .$ $V_{2}$ and $H_{2}$ are defined in expressions ((ref)) and ((ref)), respectively. \end{enumerate}

KimPollard1990 and SeoOtsu2018 derive the cube-root asymptotics for a class of estimators by means of empirical processes. For a comprehensive treatment on this technique, see vdVaartWellner2000. Our estimators fall into this category. In particular, they are more closely related to SeoOtsu2018. The main body of the proof for Theorem (ref) verifies the technical conditions in SeoOtsu2018, applies their asymptotics results to our estimators, and calculates the technical terms needed for the asymptotics such as $V_{1},H_{1},V_{2},$ and $H_{2}.$

Note that the asymptotics of $\hat{\gamma}$ are the same as in the case where the true value of $\beta$ is used. Intuitively, $\hat{\beta}$ converges to $\beta$ faster than $\hat{\gamma}$ does to $\gamma,$ and the objective function ((ref)), after proper normalization, uniformly converges to the limit over a compact set of $\left(b',r\right)'$ around $\left(\beta',\gamma\right)'$. The details can be found in the proof of Theorem (ref), which is presented in Appendix (ref).

Inference

The asymptotic distributions of $\hat{\beta}$ and $\hat{\gamma}$\ are complicated and do not have an analytical form. As a result, inference using the asymptotic distribution directly is difficult to implement. One may consider smoothing the objective functions, in the spirit of Horowitz1992, to attain faster rates of convergence and asymptotic normality.\footnote{See also Kyriazidou1997 and charlier1997limited.} However, this requires selecting additional kernel functions and tuning parameters, and then computing consistent estimates for asymptotic variances. As an alternative, we seek to use more direct sampling methods (e.g., bootstrap). Unfortunately, AbrevayaHuang2005 have proved the inconsistency of the classic bootstrap for the MS score estimators. We expect that the classic bootstrap does not work for our estimators either.

For the ordinary MS estimator, valid inference can be conducted using subsampling (DelgadoEtal2001), $m$-out-of-$n$ bootstrap (LeePun2006), the numerical bootstrap (HongLi2020), and a model-based bootstrap procedure that analytically modifies the criterion function (CattaneoEtal2020), among other procedures.\footnote{The case-specific, smooth bootstrap method proposed by PatraEtal2018 is also valid for the MS estimator of Manski1975, Manski1985. However, this method is difficult to generalize to our case.} These methods, with certain modifications, can be justified to be valid for our estimators.

Monte Carlo evidence demonstrated in HongLi2020 suggests that their proposed approach outperforms the subsampling and the $m$-out-of-$n$ bootstrap in finite samples. Based on these results, we focus on the numerical bootstrap. We provide a brief discussion on the classic bootstrap and the $m$-out-of-$n$ bootstrap in Appendix E.\footnote{We show in Appendix E that the classic bootstrap is not consistent for our estimators (Appendix E.2), while the $m$-out-of-$n$ bootstrap is still valid (Appendix E.3). Note that we re-use some notation in this appendix for notational convenience. To avoid confusion, all notation in each subsection is specific to the procedure discussed in that subsection.}

We next introduce some additional notation. Let $( y_{j}^{T\ast\prime},x_{j} ^{T\ast\prime}) ^{\prime},$ $j=1,...,n,$ be a random sample drawn with replacement from the collection of the sample values $\left( y_{1}^{T\prime} ,x_{1}^{T\prime}\right) ^{\prime},$ $\left( y_{2}^{T\prime},x_{2}^{T\prime }\right) ^{\prime},$ $...,$ $\left( y_{n}^{T\prime},x_{n}^{T\prime}\right) ^{\prime}.$ Let $\xi_{j}^{\ast}\left( b\right) $ denote $\xi\left( b\right) $ evaluated at $( y_{j}^{T\ast\prime},x_{j}^{T\ast\prime}) ^{\prime}$, specifically, \[ \xi_{j}^{\ast}\left( b\right) \equiv1\left[ y_{j0}^{\ast}=y_{j2}^{\ast }=y_{j4}^{\ast}\right] \left( y_{j3}^{\ast}-y_{j1}^{\ast}\right) \left( 1\left[ x_{j31}^{\ast\prime}b>0\right] -1\left[ x_{j31}^{\ast\prime} \beta>0\right] \right) . \] Similarly, we define $\varsigma_{nj}^{\ast}\left( r,b\right) $ as

align*[align* omitted — 672 chars of source]

Numerical Bootstrap

HongLi2020 develop a numerical bootstrap procedure for cases in which the classic bootstrap does not work. They demonstrate that their method works for a class of M-estimators that converge at rate $n^{a}$ for some $a\in(1/4,1].$ The estimator $\hat{\beta}$ proposed in Section (ref) fits their framework directly, but $\hat{\gamma}$ does not. With a slight modification of their proof, we show that the numerical bootstrap also works for $\hat{\gamma}$.

The numerically bootstrapped $\hat{\beta}^{\ast}$ and $\hat{\gamma}^{\ast}$ are constructed from

equation[equation omitted — 307 chars of source]

and

equation[equation omitted — 368 chars of source]

where $\varepsilon_{n}\rightarrow0,$ $n\varepsilon_{n}\rightarrow\infty$ ,$\ $and $( y_{j}^{T\ast\prime},x_{j}^{T\ast\prime}) ^{\prime},$ $j=1,...,n,$ are drawn independently from the collection of the sample values $\left( y_{1}^{T\prime},x_{1}^{T\prime}\right) ^{\prime},\left( y_{2}^{T\prime} ,x_{2}^{T\prime}\right) ^{\prime},...,\left( y_{n}^{T\prime},x_{n}^{T\prime }\right) ^{\prime}$, with replacement$.$ $\varepsilon_{n}^{-1}$ plays a similar role as the $m$ in the $m$-out-of-$n$ bootstrap procedure. For $\hat{\gamma}^{\ast}$, we additionally require $\varepsilon_{n}^{-1} h_{n}\rightarrow\infty$ and $\varepsilon_{n}^{-1}h_{n}^{4}\rightarrow0$. Following the same arguments as in the discussion below ((ref)) for $\hat{\beta}$, the maximizer $\hat{\beta}^{\ast}$ exists, but its uniqueness is not guaranteed due to the non-smoothness of its step objective function. The second term in ((ref)) can be shown to be $o_P(1)$. The first term, sharing a similar structure to the objective function ((ref)), dominates in ((ref)). As a result, we anticipate $\hat{\beta}^{\ast}$ to be close to $\beta$ as $n\rightarrow\infty$. Similar arguments apply to $\hat{\gamma}^{\ast}$.

We claim that \[ \varepsilon_{n}^{-1/3}\left( \hat{\beta}^{\ast}-\hat{\beta}\right) \overset{d}{\rightarrow}\arg\max_{\boldsymbol{s}\in\mathbb{R}^{K}}\left( \frac{1}{2}\boldsymbol{s}^{\prime}V_{1}\boldsymbol{s}+W_{1}\left( \boldsymbol{s}\right) \right) \] and \[ \left( \varepsilon_{n}^{-1}h_{n}\right) ^{1/3}\left( \hat{\gamma}^{\ast }-\hat{\gamma}\right) \overset{d}{\rightarrow}\arg\max_{s\in\mathbb{R} }\left( \frac{1}{2}V_{2}s^{2}+W_{2}\left( s\right) \right) , \] where $W_1$ and $W_2$ are mean zero Gaussian processes with covariance kernels $H_1$ and $H_2$, respectively. An outline of the proof of why the numerical bootstrap works and the way to modify the proof in HongLi2020 to accommodate $\hat{\gamma}$ is provided in Appendix E.1.

Procedures in Details

We investigate the finite-sample properties of the bootstrap method discussed in Section (ref) using Monte Carlo experiments in Section (ref), and defer the discussion on the choices of their tuning parameters to Section (ref). Here, we provide the algorithm for constructing the 95% confidence intervals (CIs) for $\beta$ and $\gamma$.

The numerical bootstrap proceeds as follows.

enumerate• Draw $( y_{j}^{T\ast\prime},x_{j}^{T\ast\prime}) ^{\prime},$ $j=1,...,n,$ independently, with replacement, from the original sample. • Obtain $\hat{\beta}^{\ast}$ and $\hat{\gamma}^{\ast}$ from equations ((ref)) and ((ref)). • Repeat Steps 1 and 2 for $B$ times independently and obtain a sequence of $(\hat{\beta}^{\ast},\hat{\gamma}^{\ast})$, say, $\{ (\hat{\beta}_{{} }^{\ast\left( b\right) },\hat{\gamma}_{{}}^{\ast\left( b\right) })\} _{b=1}^{B}.$ • Let $Q_{\hat{\beta}^{\ast}}\left( \tau\right) $ denote the $\tau$-th quantile of $\{ \hat{\beta}_{{}}^{\ast\left( b\right) }\} _{b=1}^{B},$ $0\leq\tau\leq1$. Define $Q_{\hat{\gamma}^{\ast}}\left( \tau\right) $ analogously. The 95% CIs for $\beta$ and $\gamma$ are constructed, respectively, as \[ \left[ \hat{\beta}-n^{-1/3}\cdot\varepsilon_{n}^{-1/3}( Q_{\hat{\beta }^{\ast}}\left( 0.975\right) -\hat{\beta}) ,\text{ }\hat{\beta }-n^{-1/3}\cdot\varepsilon_{n}^{-1/3}( Q_{\hat{\beta}^{\ast}}\left( 0.025\right) -\hat{\beta}) \right] \] and \[ \left[ \hat{\gamma}-n^{-1/3}\cdot\varepsilon_{n}^{-1/3}\left( Q_{\hat {\gamma}^{\ast}}\left( 0.975\right) -\hat{\gamma}\right) ,\text{ } \hat{\gamma}-n^{-1/3}\cdot\varepsilon_{n}^{-1/3}\left( Q_{\hat{\gamma}^{\ast }}\left( 0.025\right) -\hat{\gamma}\right) \right] . \]

Monte Carlo Experiments

Simulation Setup

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 HK. Specifically, this design (referred to as Design 1 hereafter) is specified as follows:

align*[align* omitted — 251 chars of source]

where

itemize$\beta\equiv(\beta_{1},\beta_{2})^{\prime}=(1,1)^{\prime}$ and $ \gamma=-1$, • $x_{it,j}=\frac{\sqrt{15}}{4}u_{it,j}+\frac{1}{4}u_{it,3},j=1,2$, $ \left( u_{it,1},u_{it,2},u_{it,3}\right)\overset{d}{\sim }N\left( 0_{3\times 1},I_{3\times 3}\right)$, and $\left( u_{it,1},u_{it,2},u_{it,3}\right) $ are i.i.d. across $i$ and $t$,\footnote{Note that all covariates are correlated with each other and have a standard deviation of one.} • $\alpha _{i}=\left( x_{i0,2}+x_{i1,2}+x_{i2,2}+x_{i3,2}+x_{i4,2}\right) /5$, • $\epsilon _{it}\overset{d}{\sim }\left( \pi ^{2}/3\right) ^{-1/2}\cdot \text{Logistic}\left( 0,1\right) $ and are i.i.d. across $i$ and $t$, and • $\left( u_{\cdot ,1},u_{\cdot ,2},u_{\cdot ,3}\right) $ and $ \epsilon _{\cdot }$\ are independent of each other.

In the second design (hereafter, Design 2), the model and the coefficients are the same as in Design 1, but $x_{\cdot ,1}$ and $x_{\cdot ,2}$ are autocorrelated over time. Specifically, we have

itemize$x_{i0,j}=\frac{\sqrt{15}}{4}u_{i0,j}+\frac{1}{4}u_{i0,3},$ $j=1,2,$ and $x_{it,j}=\frac{1}{2}x_{it-1,j}+\frac{\sqrt{3}}{2}\left( \frac{\sqrt{15}}{4}u_{it,j}+ \frac{1}{4}u_{it,3}\right) ,$ $j=1,2$ for all $t\geq 1$, where $\left( u_{it,1},u_{it,2},u_{it,3}\right) \overset{d}{\sim }N\left( 0_{3\times 1},I_{3\times 3}\right) $ and $\left( u_{it,1},u_{it,2},u_{it,3}\right) $ are i.i.d. across $i$ and $t$, • $\left( u_{\cdot ,1},u_{\cdot ,2},u_{\cdot ,3}\right) $ and $ \epsilon _{\cdot }$\ are independent of each other.

Note that the setup of Design 2 violates both Assumption \hyperref[Assumption:SI] {SI} and the exchangeability condition stated in Assumption \hyperref[Assumption:SD] {SD}. We conduct this Monte Carlo study to develop insight into the practical consequences of the failure of these sufficient (but not necessary) conditions. That is, we examine the extent to which serial dependence in exogenous covariates may affect the identification.

In the third to fifth designs (Designs 3, 4, and 5, respectively), the setup is the same as that in Design 1, except that we add one, two, and three more covariates (in Designs 3, 4, and 5, respectively) to examine how our estimators perform in higher-dimensional designs. Specifically, in Design $k,$ $k=3,4,$ and $5,$

align*[align* omitted — 305 chars of source]

where

itemize$\beta \equiv (\beta _{1},\beta _{2},...,\beta _{k})^{\prime }=(1,1,...,1)^{\prime }$ and $\gamma =-1$, • $x_{it,j}=\frac{\sqrt{15}}{4}u_{it,j}+\frac{1}{4} u_{it,k+1},j=1,2,...,k$, $\left( u_{it,1},u_{it,2},...,u_{it,k+1}\right) \overset{d}{\sim }N\left( 0_{(k+1)\times 1},I_{(k+1)\times (k+1)}\right) $, and $\left( u_{it,1},u_{it,2},...,u_{it,k+1}\right) $ are i.i.d. across $i$ and $t$, • $\alpha _{i}=\left( x_{i0,2}+x_{i1,2}+x_{i2,2}+x_{i3,2}+x_{i4,2}\right) /5$, • $\epsilon _{it}\overset{d}{\sim }\left( \pi ^{2}/3\right) ^{-1/2}\cdot \text{Logistic}\left( 0,1\right) $ and are i.i.d. across $i$ and $t$, and • $\left( u_{\cdot ,1},u_{\cdot ,2},...,u_{\cdot ,k+1}\right) $ and $ \epsilon _{\cdot }$\ are independent of each other.

We also explore the impact of the serial dependence of $x_{it}$ on the estimation and inference for the models in Designs 3 through 5 in Appendix F. We adopt a similar method to Design 2 for this analysis, which offers additional insights into the robustness and performance of our proposed method.

For the estimation of $\beta$, we adopt the objective function ((ref)). To estimate $\gamma$, we use the objective function ((ref)) with the Epanechnikov kernel function. That is,

equation*[equation* omitted — 123 chars of source]

which satisfies Assumption (ref) with a compact support. We discuss the choice of bandwidth sequence $h_{n}$ in Section (ref).

For inference, we investigate the finite-sample performance of the numerical bootstrap (Section (ref)). The 95% CIs are obtained from $B=199$ independent draws and estimations. See Section (ref) for the details of the implementation.

Recall that only the observations with $\left\{ y_{i0}=y_{i2}=y_{i4}\text{ and }y_{i1}\neq y_{i3}\right\} $ are used to estimate $\beta$. In all designs, the effective observations, which are useful for estimating $\beta$, comprise about $14\%$ of the whole sample. Similarly, for $\gamma$, only observations with either $\left\{ y_{i1}\neq y_{i2}\text{ and }y_{i0}\neq y_{i3}\right\} $ or $ \left\{ y_{i2}\neq y_{i3}\text{ and }y_{i1}\neq y_{i4}\right\} $ are useful. In all designs, about $31\%$ to $39\%$ of the observations are effective for $\gamma$. For each design, we consider sample sizes of 2,500, 5,000, 10,000, and 20,000. All the estimation and inference results (based on 199 draws and estimation) presented in this section are based on 1,000 replications of each design and each sample size.

Furthermore, we compare our method with the parametric (Logit) and semiparametric (distribution-free) estimators of HK. We use the objective functions linked to these two estimators, as defined in Section 4.1 of HK (for panel data where $T>3$), to implement their methods. To facilitate the comparison, we apply the same scale normalization, specified in Assumption \hyperref[Assumption:HK]{A}(e), to both HK's two estimators and our own. Specifically, we normalize the vector of $\beta$'s to have a Euclidean norm of one. Note that HK's Logit estimator does not require scale normalization for the preference coefficients, because it assumes that the error terms follow a standard logistic distribution. Therefore, if we choose to apply scale normalization to $\beta$, we should also estimate a scale parameter in the Logit model to regain one degree of freedom in the parameter space. For example, for two adjacent time periods $t$ and $t+1$ in Designs 1 and 2, the log-likelihood function is written as

align*[align* omitted — 565 chars of source]

where we impose restriction $\sqrt{b _{1}^{2}+b _{2}^{2}}=1$ and add the scale parameter of the logistic distribution, $s>0$, to estimate.\footnote{HK normalize $s$ to one, but do not require $(b_1, b_2)$ to be on the unit circle. Both methods of scale normalization are equivalent.}

commentFinally, we compare our method with the parametric (assuming Logit) and semiparametric (without assuming Logit) methods in HK. For Design 1, the parametric estimator in HK is obtained by \begin{equation*} \left( \hat{\beta}_{1},\hat{\beta}_{2},\hat{\gamma},\hat{s}\right) =\arg \max_{\left( b_{1},b_{2},r,s\right) ,\sqrt{b_{1}^{2}+b_{2}^{2}} =1,s>0}\sum_{i=1}^{n}\left( Q_{i1}\left( b_{1},b_{2},r,s\right) +Q_{i2}\left( b_{1},b_{2},r,s\right) \right) \end{equation*} where \begin{eqnarray*} Q_{it}\left( b_{1},b_{2},r,s\right) &\equiv &1\left[ y_{it}+y_{it+1}=1\right] \sigma _{n}^{-2}\mathcal{K}\left( \frac{x_{it+1,1}-x_{it+2,1}}{\sigma _{n}} \right) \mathcal{K}\left( \frac{x_{it+1,2}-x_{it+2,2}}{\sigma _{n}}\right) \\ &&\times \log \left( \frac{\exp \left( \left[ \left( x_{it,1}-x_{it+1,1}\right) b_{1}+\left( x_{it,2}-x_{it+1,2}\right) b_{2}+r\left( y_{it-1}-y_{it+2}\right) \right] /s\right) ^{y_{it}}}{1+\exp \left( \left[ \left( x_{it,1}-x_{it+1,1}\right) b_{1}+\left( x_{it,2}-x_{it+1,2}\right) b_{2}+r\left( y_{it-1}-y_{it+2}\right) \right] /s\right) }\right) . \end{eqnarray*} We normalize $\sqrt{\beta _{1}^{2}+\beta _{2}^{2}}$ to 1, and we add the scale parameter $s,$ $s>0$, into the estimation.\footnote{ In HK, they normalize the scale parameter, $s,$ to 1 instead. The normalization of $\sqrt{\beta _{1}^{2}+\beta _{2}^{2}}$ to 1 is equivalent.} Since $T=4$ here$,$ HK estimators are able to utilize observations from periods 1-- 4, reflected in $Q_{i2}\left( b_{1},b_{2},r,s\right) .$ The objective functions for other designs can be similarly written down. For Design 1, the semiparametric estimator in HK is obtained by \begin{equation*} \left( \hat{\beta}_{1},\hat{\beta}_{2},\hat{\gamma}\right) =\arg \max_{\left( b_{1},b_{2},r\right) ,\sqrt{b_{1}^{2}+b_{2}^{2}} =1}\sum_{i=1}^{n}\left( \tilde{Q}_{i1}\left( b_{1},b_{2},r\right) +\tilde{Q} _{i2}\left( b_{1},b_{2},r\right) \right) , \end{equation*} with \begin{eqnarray*} \tilde{Q}_{it}\left( b_{1},b_{2},r\right) &\equiv &\sigma _{n}^{-2}\mathcal{K }\left( \frac{x_{it+1,1}-x_{it+2,1}}{\sigma _{n}}\right) \mathcal{K}\left( \frac{x_{it+1,2}-x_{it+2,2}}{\sigma _{n}}\right) \left( y_{it+1}-y_{it}\right) \\ &&\times sgn\left( \left( x_{it+1,1}-x_{it,1}\right) b_{1}+\left( x_{it+1,2}-x_{it,2}\right) b_{2}+r\left( y_{it+2}-y_{it-1}\right) \right) . \end{eqnarray*} $\tilde{Q}_{i2}$ plays a similar role to $Q_{i2}$. Again, the objective functions for other designs can be similarly written down.

Tuning Parameters and Computation

There is only one tuning parameter used for estimation, namely, $h_{n},$ in the objective function ((ref)). In Assumption (ref), we restrict $ nh_{n}^{4}\rightarrow0$, so that the bias term (of order $h^{2}_{n}$) is a small order term of $\left( nh_{n}\right) ^{-2/3}$. Because the convergence rate of $\hat{\gamma}$ is $\left( nh_{n}\right) ^{-1/3}$, the condition, $ nh_{n}^{4}\rightarrow0,$ makes the bias term much smaller than the convergence rate. To attain a faster convergence rate, we set $h_{n}$ as large as possible, and thus set $h_{n}=n^{-1/4}\left( \log n\right) ^{-1}$.

For the numerical bootstrap, we have one additional tuning parameter $\varepsilon _{n}.$ As recommended in HongLi2020, we set $\varepsilon_{n}$ proportional to $n^{-2/3}\log n$ for the inferences of $\hat{\beta}$ and $ \hat{\gamma}.$ Apparently, $\varepsilon_{n}$ of this order satisfies the additional requirements for $\hat{\gamma}_{{}}^{\ast}$\ that $ \varepsilon_{n}^{-1}h_{n}\rightarrow\infty$ and $ \varepsilon_{n}^{-1}h_{n}^{4}\rightarrow0.$ To check how sensitive the procedure is to the choice of $\varepsilon_{n}$, we conduct the procedure with $\varepsilon_{n}=c\cdot n^{-2/3}\log n$ and $c=0.8,0.9,1.0,1.1,$ and $ 1.2.$

Following the recommendation in HK, we adopt the bandwidth $\sigma_n = c \cdot n^{-1/(4+k)}$ for HK's estimators, where $k$ is the dimension of $x_{it}$. We conduct experiments with $c=1,2,3,4$, and report the simulation results corresponding to $c=3$. For this value, the HK estimators of $\gamma $ exhibit the smallest bias and have relatively smaller root mean squared errors among all tested values.

For all simulation designs, we employ the DE algorithm (using the DEoptim package in R, see Remark (ref)) to compute our proposed estimators and those of HK. To ensure efficient convergence and robustness, we set the lower and upper bounds for searching each parameter to $[-3, 3]$, the maximum number of iterations to 500, and the relative convergence tolerance to $10^{-8}$. The DE algorithm does not require explicit initial values; instead, it randomly assigns $NP \times (\text{the number of parameters})$ initial values, where we adopt the default value of $NP=10$ in DEoptim. Additionally, we use the default settings for the other algorithm controls. In all simulation runs for each estimator, we consistently observe successful convergence of the algorithm. All estimators can be computed very quickly for all sample sizes considered. Using an Intel\textsuperscript{\tiny\textregistered} Core\textsuperscript{\tiny TM} i7-4790 processor, each replication takes only seconds to complete.

Simulation Results

We normalize the preference coefficients $\beta$ on exogenous covariates to $ 1$ in Euclidean norm. Because of this normalization, we lose one degree of freedom in the parameter space. As a result, we report only the results for $\left( \beta_{2},\gamma\right) $ in Designs 1 and 2 and the results for $\left( \beta_{2}, ..., \beta_{k},\gamma\right) $ in Design $k$, for $k=3,4,$ and $5$.

We report the mean bias (BIAS), the standard deviation (STD), the median absolute deviation (MAD), and the root mean squared error (RMSE) for $\hat{ \beta}$ and $\hat{\gamma}.$ All results are expressed as percentages of the true values of the parameters, so that the results are independent of how we normalize the parameters.\footnote{We thank the co-editor for this suggestion.} For inference, we report the coverage rates (COV) of the true values and lengths (LEN) of the 95% CIs for the numerical inference procedure for our method only. Furthermore, we report the computation time (TIME) for one replication of each estimator for Design 1, and omit this for other designs, owing to their similarity.

Results for Design 1 are reported in tables numbered “1”, and so on for other designs. We report the performance of the estimators and the numerical bootstrap procedure in tables labeled “A” and “B”, respectively. For example, Table 1A reports the performance of the estimators for Design 1. The results for our estimator are denoted as “OY” in the tables. The parametric and semiparametric estimators in HK are denoted as “HK1” and “HK2”, respectively. Due to space limitations, only the results of Designs 1 and 2 are presented in Appendix (ref). The results of Designs 3--5, along with additional simulation studies, are reported in Appendix F of the Online Supplement. In what follows, we briefly summarize our findings.

The RMSEs of $\hat{\beta}$ and $\hat{\gamma}$ become smaller as the sample size increases in all designs, with the RMSE of $\hat{\gamma}$ slightly greater than that of $\hat{\beta}$. This shows the consistency of our estimators, though the rates of convergence are clearly slower than $\sqrt{n}$. The numerical bootstrap inference procedures perform reasonably well in all designs. In general, they yield shrinking CIs with coverage rates approaching 95% as the sample size grows. The coverage rates of these CIs are greater than 90%, but are slightly lower than 95% in most cases. The coverage rates of the CIs for $\gamma$ do not perform as well as those for $ \beta, $ which is not surprising, considering the complication of using two tuning parameters. The inference procedure is not very sensitive to the choice of tuning parameters.

Despite Design 2 not satisfying Assumption \href{assumptionSI}{SI} or \href{assumptionSD}{SD}, our proposed estimators still perform reasonably well in this setting.\footnote{These assumptions guarantee that the utility index $x_t'\beta$ rank orders the (conditional) probabilities in ((ref)). Relaxing them could potentially invalidate the “if and only if” result in equation ((ref)) and cause bias.} Surprisingly, its performance is similar to that of Design 1, where these assumptions are completely fulfilled. This finding indicates that our method exhibits certain robustness, and can function effectively even when the two sufficient identifying assumptions are not met. Furthermore, the results from Designs 3--5 provide evidence supporting our asymptotic analysis. In these designs, we observe that the convergence rates of our estimators remain relatively stable as the dimension of the model increases.

The HK1 estimator (HK's Logit estimator) performs the best for Designs 1 and 2, which is not surprising, because the error terms are scaled logistic. Our proposed estimators exhibit higher RMSEs compared with those of HK1, approximately more than twice for these designs. For Designs 1 and 2, the HK2 estimator demonstrates finite-sample performance similar to our proposed method.\footnote{In particular, our method exhibits smaller RMSEs for $\gamma$, and larger RMSEs for $\beta$ compared to HK2.} However, when the number of regressors increases in Design 3, our estimator outperforms HK2, particularly in estimating the parameter $\gamma$. In this setting, the RMSEs of our estimators become about 50% greater than those of HK1. In Design 4, where there are two more regressors than in Design 1, our estimators' RMSEs are comparable with those of HK1 for all sample sizes considered. Notably, for sample sizes of $n=10,000$ and $20,000$, our RMSEs are about 40% lower than those of HK2, demonstrating the advantage of our method in high-dimensional settings. In Design 5, with three more regressors than Design 1, our RMSEs are slightly lower than those of HK1 and only half of HK2's RMSEs for sample sizes of $n=10,000$ and $20,000$. These findings highlight the favorable properties of our method, particularly its resilience to the curse of dimensionality.

We also conduct additional simulations to assess how our estimators perform on a relatively small sample size of $n=1,000$ for Design 1. The results presented in Table 1C suggest that our estimators do not perform poorly in terms of parameter estimation. The RMSEs are 30%-35% of the true parameter values, suggesting reasonable accuracy in estimation, despite the smaller sample size. For inference, we report only the results for $c = 1$. We can see the CIs have lower coverage of approximately 85%. In conclusion, our estimators perform reasonably well for this sample size, but the CIs may be too short.

A final note is that the serial dependence of $x_{it}$ has limited impact on the estimation and inference, as demonstrated by the results associated with Design 2 and additional simulation studies for higher-dimensional designs presented in Appendix F of the Online Supplement.

Conclusions

This paper presents new identification results for preference parameters in panel data binary choice models that allow for both fixed effects (Heckman's “spurious” state dependence) and lagged dependent variables (“true” state dependence). The same semiparametric random utility framework as in HonoreKyriazidou2000 is considered. A key innovation in this paper is the assertion that, given additional restrictions on the dynamic process of observed covariates and the tail behavior of the error distribution, the point identification no longer needs element-by-element matching of regressors over time, in contrast to the method proposed in HonoreKyriazidou2000. Our approach requires a minimum panel length of five ($T\geq 4$), which fits in most empirical settings. Our identification arguments motivate a two-step estimation procedure, adapting Manski's MS estimator. The proposed estimators are consistent with rates of convergence independent of the model dimension, unlike the estimator proposed in HonoreKyriazidou2000. We further derive the limiting distributions of the proposed estimators, which are non-Gaussian, aligning with existing literature. We justify the use of several bootstrap procedures for conducting statistical inference. A Monte Carlo study indicates that our estimators and inference procedures perform well in finite samples.

This paper leaves some open questions for future research. For example, it might be worthwhile extending the framework in this paper to study the identification with more than one lag of the dependent variable or the identification in panel data multinomial response models.