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.
81,582 characters · 9 sections · 138 citation commands
Revisiting Panel Data Discrete Choice Models with Lagged Dependent Variables
\affil[1]{Google} \affil[2]{School of Economics, University of Queensland} \affil[3]{Research School of Economics, Australian National University}
{\it Keywords:} Dynamic Binary Choice Model, Fixed Effects, Identification at Infinity, Maximum Score Estimation
In this paper, we propose new identification and estimation methods for panel data binary choice models with fixed effects and “dynamics” (lagged dependent variables). Specifically, suppose that there are $n$ individuals and $T+1$ time periods, $\{0,1,...,T\} $. In each time period $t\in \{1,...,T\}$, each individual $i$ makes a choice $y_{it}\in \{0,1\}$ according to the following latent utility model:
where $\alpha _{i}$ is an entity fixed effect absorbing all relevant time-invariant factors, $y_{it-1}$ is the lagged dependent variable, $ (x_{it},z_{it})$ is a $(p+1)$-vector of time-varying covariates, and $ \epsilon _{it}$ is an idiosyncratic error term. We separate $z_{it}$ from other covariates because, as it will be clear in the next section, we assign it a crucial role in the identification at infinity. In panel data literature, $ y_{it-1}$ is often called “state dependence”, and $\alpha _{i}$ is referred to as “unobserved heterogeneity” or “spurious” state dependence (see Heckman1981b, Heckman1981a). In model ((ref)), $ (y_{it},x_{it},z_{it})$ along with the “initial status” $y_{i0}$ are observed in the data, whereas $\alpha _{i}$ and $\epsilon _{it}$ are not observable to the econometrician. Note that we do not specify model ((ref)) in the initial period 0. This paper studies the identification and estimation of the preference parameter $\theta :=(\gamma ,\beta ,\varpi )\in \mathbb{R}^{p+2}$ in “short” panel settings, i.e., $n\rightarrow \infty $ and $T<\infty $.
In line with the vast literature on panel data models with entity fixed effects, we do not impose any parametric restrictions on the distribution of $\alpha_i$ conditional on the initial choice $y_{i0}$ and observed covariates in model ((ref)). The prevalent methods for such models assume that $\epsilon_{it}$ are independently and identically distributed (i.i.d.) with a logistic distribution. ArellanoHonore2001, honore2021identification, and Hsiaobook review various conditional likelihood approaches based on these parametric assumptions on $\epsilon _{it}$. Recent advances in the literature focus on constructing moment conditions for variants of dynamic Logit models. Representative works include honore-w, dobronyi-gu-kim, kitazawa2022transformations, and dano2023transition, among others.
Without making distributional assumptions on $\epsilon _{it}$, manski-87 establishes the semiparametric identification of model ((ref)) that includes covariates, but not $y_{it-1}$. honore-k extend this approach to include both $y_{it-1}$ and covariates in the model, showing that model ((ref)) with $T\geq 3$ can be identified under exogeneity and serial dependence assumptions stronger than those in 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. OuyangYang2024binary demonstrates that this curse of dimensionality can be mitigated by imposing certain serial dependence conditions on the covariates and by observing an extra time period. To highlight the novelty and contributions of this paper, we present a thorough comparison of our method with honore-k and OuyangYang2024binary in Appendices (ref) and (ref), respectively.
There are alternative semiparametric and nonparametric approaches to model ((ref)). hl demonstrate that model ((ref)) can be point identified if $z_{it}$ satisfies certain exclusion restrictions. More recently, ChenEtal2019 revisit this method and discuss the sufficient conditions for such exclusion restrictions. Williams2019 studies the nonparametric identification of dynamic binary choice models that satisfy certain exclusion restrictions. In the absence of excluded regressors, arist establishes informative partial identification of model ((ref)) under weak conditions. KhanEtal2020 offer a partial identification result under even milder restrictions and prove that point identification is attainable in many interesting scenarios.
This paper revisits the distribution-free identification and estimation of model ((ref)). We show that the overlapping support restrictions required by honore-k can be removed if $z_{it}$ is a free-varying covariate with full support. Our identification employs an “identification at infinity” strategy, first introduced in chamberlain1986asymptotic and heckman-90, and then applied in more recent work such as tamer2003incomplete, bajari2010identification, wan2014semiparametric, and ouyang2020semiparametric, among others. The combination of this strategy and manski-87's (manski-87) insight yields an estimator in the spirit of honore-k's (honore-k) conditional maximum score (MS) estimator, but without the need to match observed covariates over time. As a result, our estimator can accommodate flexible time effects and escape from the curse of dimensionality, in contrast to honore-k. Through extending KimPollard1990 and SeoOtsu2018, we demonstrate that our estimator converges at a rate slower than cube-root-$n$, is independent of the number of observed covariates, and has a non-standard limiting distribution. The asymptotics share similarities with those in honore-k and OuyangYang2024binary, with an important difference: the rate of convergence for our estimator depends on unknown factors, while the convergence rates of their estimators are known. We evaluate the finite-sample performance and implementability of our proposed estimator using both simulated and real-world data.
The rest of this paper is organized as follows. Section (ref) establishes the identification of $\theta$, which serves as the basis for the MS estimator presented in Section (ref). We then derive asymptotic properties of the proposed estimator in Section (ref). Results of Monte Carlo experiments are reported in Section (ref). We present an empirical illustration using the HILDA data in Section (ref). Finally, Section (ref) concludes the paper with a brief discussion on possible future research directions. All proofs, supplementary discussions, and additional simulation results are included in the Supplementary Appendix.
For ease of reference, we list the notations maintained throughout this paper here.
Suppose in model ((ref)) $\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:
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:
for each individual $i=1,\dots,n$.
In what follows, we will restrict the illustration of our approach to model ((ref)) with $T=3$ and $\varpi>0$ to ease the exposition. The condition $\varpi>0$ implies that we must know that the covariate $z_{it}$ is included in the model and that it has a positive effect on the choice probability of $y_{it}=1$. Applying our method to longer panels is straightforward, and the case with $\varpi<0$ is symmetric. In addition, we will omit the subscript $i$ in our notation whenever the context makes clear that all variables pertain to each individual. Finally, we assume a balanced panel for simplicity. Our methods are applicable to models with unbalanced panels, provided the unbalancedness is not due to endogenous attrition.
Consider two choice histories \[ C =\{y_{0}=d_{0},y_{1}=0,y_{2}=d_{2},y_{3}=1\} \textrm{ and } D =\{y_{0}=d_{0},y_{1}=1,y_{2}=d_{2},y_{3}=0\}, \] where $d_{0},d_{2}\in\{0,1\}$. The conditional probability of the choice history $C$ is equal to
where $p_{0}(\alpha,x^{T},z^{T})$ denotes the conditional probability of $y_{0}=1$. In a similar fashion,
Here, we take $d_{2}=1$ to illustrate, and the case with $d_{2}=0$ is symmetric. Suppose the support of $z_{2}$ is unbounded above. Then, for $z_2>\sigma$, where $\sigma$ is a sufficiently large positive number, these probabilities satisfy
The idea of the above is to make negligible the effect of $y_{1}$ on $y_{2}$, i.e., $F_{\epsilon|\alpha}(\alpha+\gamma+x_{2}'\beta+\varpi z_{2})\approx F_{\epsilon|\alpha}(\alpha+x_{2}'\beta+\varpi z_{2})$ ($\approx1$), by letting $z_2$ be sufficiently large.
Suppose $F_{\epsilon|\alpha}(\cdot)$ is strictly increasing. Then, when $d_{2}=1$, equation ((ref)) implies that
holds for $z_2>\sigma$ as $\sigma\rightarrow +\infty$, where $\text{sgn}\{\cdot\}$ is the sign function, which is equal to 1 if the expression inside the brackets is strictly positive, to 0 if the expression inside the brackets is zero, and to $-1$ if the expression inside the brackets is strictly negative.
Equation ((ref)) reveals that when $z_{2}$ is sufficiently large and $d_{2}=1$, the likelihood of observing event $C$ exceeds that of observing event $D$ if and only if $\gamma(d_{2}-d_{0})+(x_{3}-x_{1})^{\prime}\beta+\varpi(z_{3}-z_{1})>0$. In other words, the sign of $\gamma(d_{2}-d_{0})+(x_{3}-x_{1})^{\prime}\beta+\varpi(z_{3}-z_{1})$ determines the rank order of the conditional probabilities of events $C$ and $D$. Our focus on the subsample with $y_{3}\neq y_{1}$ aligns with manski-87 in forming his MS estimator. The distinction lies in our additional conditioning event of $z_{2}$ being large. It is worth noting that the same identification equation holds true when $-z_{2}$ is sufficiently large and $d_{2}=0$.
A natural way to build a population objective function based on equation ((ref)) is to define
with $w>0$\ for $d_{2}=1$. By a symmetric argument, define
with $w>0$ for $d_{2}=0$. Note that equation ((ref)) implies that $\bar{Q}_{1}(\gamma,\beta,\varpi)\geq\bar{Q}_{1}(r,b,w)$ and $\bar{Q}_{2}(\gamma,\beta,\varpi)\geq\bar{Q}_{2}(r,b,w)$ for all $(r,b,w)\neq(\gamma,\beta,\varpi)$. Establishing that $\theta:=(\gamma,\beta,\varpi)$ is the unique maximum of either objective function ((ref)) or ((ref)) would confirm the point identification of these coefficients. The following conditions are sufficient for this.
Assumptions A1(i)--(iii) place the same restrictions on the joint distribution of $(\alpha,\epsilon^{T},x^{T},z^{T})$ as honore-k, which implies that the unobserved heterogeneity (entity fixed effects) $\alpha$ picks up both the autocorrelation in the unobservables and the dependence between explanatory variables and unobservables. As a result, $\epsilon_{t}$ is independent of $(x^{T},z^{T},y^{t-1})$ conditional on $\alpha$ for all $t\in\mathcal{T}$. Assumption A1(iv) is a regularity condition to guarantee that any possible sequence of $y^{T}$ has a positive probability to occur.
Assumption A2 is a pivotal assumption that enables the “identification at infinity” approach, and when combined with Assumption A1, it establishes the identification equation ((ref)). It is clear from the derivation of equation ((ref)) that relaxing this assumption may require additional restrictions on the parameter space $\Theta$, the support of $x_{2}$, and the distribution of $(\epsilon,\alpha)$.
The support and continuity restrictions on $\xi_{31}$ imposed by Assumption A3 are common for the family of MS-type estimators, which are required to achieve the point identification instead of a set identification. See, e.g., manski-75,manski-85,manski-87, horowitz, honore-k, Fox2007, ShiEtal2018, YanYoo2019, and khan2021inference, among others. Given the significance of Assumptions A2 and A3 in both our theoretical results and empirical application, we provide further discussion on them in Appendix (ref).
Assumption A4(i) is a familiar full-rank condition. Note that Assumptions A3 and A4(ii) require $(x_{3}-x_{1},z_{3}-z_{1})$ to have sufficient variation conditional on $\alpha$ and event $\{z_{2}>\sigma\}$ or $\{z_{2}<-\sigma\}$. Assumption A5 applies the scale normalization and restricts the search of $\theta$ in a compact set, which also facilitates the asymptotic analysis of our estimator proposed in the next section.
Additionally, we compare the key identification assumptions imposed in honore-k and OuyangYang2024binary, along with other aspects, with our method in Appendices (ref) and (ref), respectively.
Our identification results are stated in the following theorem and the proof of which is provided in Appendix (ref).
Applying the analogy principle, the population objective functions ((ref)) and ((ref)) translate into MS estimation procedures (manski-75,manski-85,manski-87).
Assume a random sample of $n$ observations is drawn from model ((ref)) that satisfies Assumption A. Let $\vartheta :=(r,b,w)\in \mathbb{R}^{p+2}$ and $ \sigma _{n}\rightarrow \infty $ as $n\rightarrow \infty $. When the support of $z_{i2}$ is unbounded above, we propose the MS estimator $\hat{\theta}_{n} $ of $\theta $ maximizing the following objective function over the parameter space $\Theta $:
When the support of $z_{i2}$ is unbounded below, one can instead define $ \hat{\theta}_{n}$ with objective function
If the support of $z_{i2}$ is unbounded both above and below, the objective function can be a combination of ((ref)) and ((ref)) such as
Note that (3.3) puts the same weight on $Q_{n1}(\vartheta)$ and $Q_{n2}(\vartheta)$, which is a generic choice and probably not optimal in specific applications. In some cases, it might be preferable to put more weight on one side if additional information, such as restrictions on the error distribution, suggests that the identification at infinity is more effective on that side, especially if $z_{2}$ has a relatively heavier tail against the error term. Since $\text{sgn}(u)=2\cdot \mathds{1}\{u>0\}-1$ almost surely for any continuous variable $u$, objective functions ((ref)) and ((ref)) are sample analogues to monotone transformations of population functions ((ref)) and ((ref)), respectively.
It is clear from expressions ((ref)) and ((ref)) that the effective sample size for the estimator $\hat{\theta}_{n}$ is controlled by the tuning parameter $\sigma_{n}$, and being similar to manski-87's (manski-87) and honore-k's ( honore-k) estimators, only “switchers” who change choices in periods 1 and 3 are used in the estimation. Besides, estimating (identifying) $\gamma$ relies on the variation in $y_{i2} - y_{i0}$, which means that we need some observations with $y_{i2}\neq y_{i0}$ and some with $y_{i2}=y_{i0}$.
Our proposed estimator $\hat{\theta}_n$ has two advantages, compared with honore-k's (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.
This section establishes the asymptotic properties of the MS estimator proposed in Section (ref). Given that objective functions ((ref)) and ((ref)) are symmetric, it suffices to only investigate the estimator $\hat{\theta}_{n}$ obtained from maximizing objective function ((ref)) requiring the support of $z_{2}$ to be unbounded above. The derivation for $\hat{\theta}_{n}$ associated with objective functions ((ref)) or ((ref)) is analogous. Additionally, for the sake of simplicity, we focus on the case where $\xi_{31}=z_{31}$ in Assumption A3.
To ensure the consistency of $\hat{\theta}_n$, we need additional technical conditions.
Assumptions B2 imposes mild restrictions on the tuning parameter $\sigma _{n} $. It is worth noting that Assumption B2(ii) indicates that the choice of $ \sigma_{n}$ depends on the tail behavior of the distribution of $z_{2}$. For example, if $z_{2}$ has a sub-exponential right tail with $P(z_{2} > \sigma_{n}) \asymp e^{-\sigma_{n}}$, then any $\sigma_n$ satisfying $1 \ll \sigma_n \leq (1-\varepsilon) \log(n)$, e.g., $\sigma_n = \log\log(n/\log n)$, meets Assumption B2(ii), for some $\varepsilon\in(1/4,1)$. However, when the distribution of $z_{2}$ has a (too) thin right tail $P(z_{2}>\sigma _{n})\asymp e^{-e^{\sigma _{n}}}$, $\sigma _{n}=\log\log(n/\log n)$ gives $nP(z_{2}>\sigma _{n})/\log n=O(1)$, violating Assumption B2(ii). Notably, $nP(z_{2}>\sigma _{n})$ essentially controls the “effective sample size” for our proposed procedure. As demonstrated in Theorem (ref), 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)) to its population analogue. We provide a set of more concrete sufficient conditions for it in Appendix (ref).
The theorem below states that the proposed procedure described in ((ref))--((ref)) gives a consistent estimator of $\theta$, whose proof is left to Appendix (ref).
We proceed to study the asymptotic distribution of the estimator $\hat{\theta}_n$. Before presenting additional technical conditions and the main results, we introduce some new notation to facilitate exposition:
Assumption C1 is standard in the literature (see, e.g., KimPollard1990 and SeoOtsu2018). This assumption implies that the maximization of $Q_{n1}(\vartheta)$ need not be exact, and any approximate maximizer close enough to the exact one will be enough for the asymptotic analysis. Assumption C2 is an implication of Assumptions A1 and A2. We list it as a separate condition here mainly because it is more directly related to our proof process presented in Appendix (ref). Assumption C3 strengthens Assumption A4(ii). Assumption C4 requires the two conditional probabilities $\kappa_{n}(\bar{\chi})$ and $\kappa ^{+}(\bar{\chi})$ to be smooth enough, which, together with Assumption C3, is important for calculating the expected value of the limiting distribution of the estimator $\hat{\theta}_{n}$.
The smoothness conditions imposed in Assumption C5(i) are standard in the literature as well. Assumption C5(ii) is made to simplify the proof process and can be relaxed to allow for unbounded support, albeit with more tedious discussions. The essential requirement here is to exclude the scenario in which $\alpha + x_2^{\prime} \beta \rightarrow -\infty$ as $z_2 \rightarrow +\infty$. Assumption C5(iii) essentially places a restriction on the relative tail behavior of the observed regressor $z_t$ and unobserved error $\epsilon_t$. As shown in the proof of Theorem (ref), this assumption ensures that the bias of the estimator $\hat{ \theta}_{n}$ shrinks sufficiently fast. If this condition is violated, the bias term dominates the distribution, and inferences are not possible. It is worth noting that such condition plays a crucial role in determining the rate of convergence of estimators based on “irregular identification” strategies including the “identification at infinity” as a special case. See khan2010irregular for an in-depth investigation on this issue.
Assumption C6, together with Assumption C5(iii), guides the selection of the tuning parameter $\sigma_{n}$. These two conditions are in the same spirit of Assumptions 8 and 8* in andrews. On one hand, since $ h_{n}=P(z_{2}>\sigma_{n}|y_{2}=1)$ controls the effective sample size of the estimation procedure, Assumption C6 implies that $\sigma_{n}$ should not increase too rapidly as $n\rightarrow\infty$, ensuring enough effective observations to control the variance of $\hat{\theta}_{n}$. On the other hand, Assumption C5(iii) suggests that $\sigma_{n}$ should grow sufficiently fast as $n\rightarrow\infty$ to lower the bias of $\hat{\theta} _{n}$.
However, there is no way to determine the optimal $ \sigma_{n}$ since this requires the knowledge of relative (unknown) tail behavior of $z_{t}$ and $\epsilon_{t}$. This feature is well known to the “identification at infinity” type of estimators, see, e.g., andrews. We suggest choices of $\sigma_n$ that satisfy both Assumptions C5(iii) and C6 in some special cases in Table (ref). From the table, there are no valid $\sigma_n$ in case (I) when $\lambda^{\prime}>\lambda$, and in case (III). The valid choices of $\sigma_n$, if exists, differ from case to case. As expected, we prefer the cases where $z_2$ possesses heavier tails than $\epsilon$. andrews share similar results, for details, see their discussions after Assumption 8*.
For practice, we propose to take $\sigma_{n}=\sqrt{\log n/2.95}$. This choice of $\sigma_n$ is valid for case (I) with $\lambda^{\prime}=2$ and $\lambda'<\lambda$, and case (II) with $\lambda=2$. Moreover, with this choice of $\sigma_n$, Assumption B2(ii) is satisfied for $z_{2}$ with $P(z_{2} > \sigma_{n}) \gtrsim e^{-\left(1-\varepsilon\right) 2.95\sigma_{n}^{2}}$ for some $\varepsilon\in (1/4,1)$. We show the finite sample properties of our estimator with this choice of $\sigma_n$ by means of simulations in Section (ref). This chosen $\sigma_n$ works well (the bias does not dominate the distribution) even in the situation that belongs to case (I) with $\lambda^{\prime}>\lambda$, where no valid $\sigma_n$ exists.
The above conditions are sufficient to characterize the asymptotic distribution of the estimator obtained by maximizing ((ref))--((ref)), as presented in the following theorem, along with its proof in Appendix (ref).
Note that Theorem (ref) does not determine the exact rate of convergence of $\hat{\theta}_n$, as $h_n$ depends on the unknown tail probabilities of $z_2$. However, the lack of this knowledge does not render statistical inference infeasible. In Section 6, we will apply the $m$-out-of-$n$ bootstrap to conduct the inference in an empirical application. We choose this method for two reasons: it is comparatively easier to implement, and it provides an estimate of the convergence rate for our estimator.
In Remark (ref), we discuss several sampling-based methods with the potential to enable statistical inference in the absence of knowledge of the exact convergence rate of the estimator.
In this section, we investigate the finite-sample performance of the proposed estimators by means of Monte Carlo experiments. We examine two designs, each with a less favorable scenario for our estimator. In these scenarios, $z_{t}$ has a thinner tail than $\epsilon_{t}$, and there is no theoretically valid $\sigma_{n}$. These are the first scenarios in both Designs 1 and 2 presented below. Despite these challenges, our estimator performs reasonably well, as the bias term does not appear to dominate the distribution.
We start by considering a benchmark design similar to that used in honore-k, but we add an additional covariate $z_{it}$ and a time trend that honore-k cannot handle. Specifically, this design (referred to as Design 1) is specified as follows:
where we set $\gamma=\beta_{1}=1$ and $\delta=1/2$. Following the discussion in Remark (ref), we normalize the coefficient on $z_{it}$ to 1 for all designs investigated in this section and Appendix (ref). We consider two scenarios. For each, we let $x_{it,1}\overset{d}{\sim}N\left(0,1\right),$ $\epsilon_{it}\overset{d}{\sim}(\pi^{2}/3)^{-1/2}\cdot$Logistic$\left(0,1\right)$ (the variance of $\epsilon_{it}$ is 1)$,$\ and $\alpha_{i}=\left(x_{i0,1}+x_{i1,1}+x_{i2,1}+x_{i3,1}\right)/4,$\ but we consider $z_{it}$ with different tail behaviors. $x_{\cdot,1},z_{\cdot},$ and $\epsilon_{\cdot}$ are independent of each other, and all covariates are i.i.d. across $i$ and $t.$ In the first scenario, we set $z_{it}\overset{d}{\sim}N(0,1),$ and denote it as “Norm”. In the second scenario, we set $z_{it}\overset{d}{\sim}\text{Laplace}(0,\sqrt{2}/2)$ (with zero mean and unit variance), and denote it as “Lap”. Note that the density function of the Laplace distribution decays like $e^{-\left\vert x\right\vert /c}$ for some constant $c$ at its tail$,$ which is heavier than the tail of the normal density.
In the second design (referred to as Design 2), the setup is the same as that in Design 1, except that we add one more covariate to examine how our estimators perform in a higher dimensional design. Specifically,
where we set $\gamma=\beta_{1}=\beta_{2}=1$ and $\delta=1/2$. Random covariates are generated as$\ x_{it,1},x_{it,2}\overset{d}{\sim}N\left(0,\sqrt{2}/2\right)$, $\epsilon_{it}\overset{d}{\sim}(\pi^{2}/3)^{-1/2}\cdot\text{Logistic}\left(0,1\right)$, and $\alpha_{i}=\sum_{t=0}^{3}(x_{it,1}+x_{it,2})/4$. Similarly, we consider two scenarios with the same $z_{it}$ as in design 1$.$ Again, $x_{\cdot,1}$, $x_{\cdot,2}$, $z_{\cdot}$, and $\epsilon_{\cdot}$ are independent of each other. To investigate only the impact of higher dimension, we set the variance of $x_{it,1}+x_{it,2}$ in Design 2 to be the same as that of $x_{it,1}$ in Design 1.
As discussed in Section (ref), we set $\sigma_{n}$ as \[ \sigma_{n}=\widehat{\text{std}\left(z_{i2}\right)}\sqrt{\log n^{*}/2.95}, \] where $\widehat{\text{std}\left(z_{i2}\right)}$ is the sample standard deviation of $z_{2}$, and $n^{*}$ is the number of “switchers”, that is, observations with $y_{3}\neq y_{1}$. The usage of $n^{*}$ is intended to provide better control over the tuning parameters, based on the features of the data. In practice, one may normalize $z_{it}$ to mean 0 and variance 1 and set $\sigma_{n}=\sqrt{\log n^{*}/2.95}$. We consider sample sizes of $n=5000,10000$, and $20000$. All the simulation results presented in this section are based on 1000 replications of each sample size. We implement MS estimations in R, using the differential evolution (DE) algorithm to attain a global optimum of the objective function. The DE algorithm, developed by storn1997differential, is capable of searching for the global optimum of a real-valued function with real-valued parameters, even if the function lacks continuity or differentiability. This algorithm has been effectively employed in calculating MS-type estimators in the literature, including Fox2007 and YanYoo2019. mullen2011deoptim provides a comprehensive introduction to the R package $\texttt{DEoptim}$, which implements the DE algorithm. We report the mean bias (MBIAS) and the root mean square errors (RMSE) of the estimates for Designs 1 and 2 in Tables (ref) and (ref), respectively.
We summarize the findings in Tables (ref) and (ref). First, the RMSEs of all parameters decrease as the sample size increases, but they converge to zero slower than the parametric rate. Second, the convergence rate is faster with a thicker-tailed $z_{i2}$, as evidenced by comparing the RMSEs from Norm to Lap. Third, the RMSE does not appear to increase for $\gamma$ and $\delta$ as we have one more covariate from Design 1 to Design 2. This confirms our theoretical findings. Note that the RMSE increases a bit for $\beta_{1}$, but this is probably due to the lower variance of $x_{\cdot,1}$ in Design 2. To investigate the sensitivity of the results to $\sigma_{n},$ we consider $\sigma_{n}=0.9\cdot\widehat{\text{std}\left(z_{i2}\right)}\sqrt{\log n^{*}/2.95}$ and $\sigma_{n}=1.1\cdot\widehat{\text{std}\left(z_{i2}\right)}\sqrt{\log n^{*}/2.95}$ (we need larger $\sigma_{n}$ to be in line with the discussion in Section (ref)), and report the corresponding results in Tables (ref) and (ref) in Appendix (ref). We note that the results are not sensitive to the choices of the tuning parameters.
In Appendix (ref), we report additional results from supplementary simulation studies. Firstly, we investigate the impact of auto-correlations of the regressors on the performance of our estimator. Additionally, we compare the performance of our estimator with those proposed by honore-k and OuyangYang2024binary in designs without the time trend term. We direct interested readers to Appendix (ref) for a more detailed discussion. Here, we provide a brief summary of these results. Our estimator still performs reasonably well with certain degrees of auto-correlations, but as expected, not as well as in Designs 1 and 2, where regressors are serially independent. Our estimator's performance is comparable to that of the semiparametric estimators proposed by honore-k and OuyangYang2024binary. It is essential to highlight that these alternative methods are not applicable in scenarios involving time trends or dummies, which are common in empirical applications. In such contexts, our approach offers a valuable alternative.
A final note is that when using observational data, the choice of $\sigma_n$ depends on the unknown tail behavior of the variable $z_2$. As there are no formal methods to determine the appropriateness of a specific $\sigma_n$, we suggest practitioners try different $\sigma_n$'s in estimation and check if the results are sensitive to different choices.
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 cheng2014measuring. We refer interested readers to cheng2014measuring and references therein for more detailed information.
In this section, we apply our MS estimator to analyze the state dependence and the impacts of government incentives on the choice to purchase PHC, using 10 waves (waves 11--20 corresponding to years 2011--2020) of the Household, Income and Labor Dynamics in Australia (\href{https://melbourneinstitute.unimelb.edu.au/hilda}{HILDA}) Survey data. Since 2011, the HILDA survey has begun recording information about respondents' enrollment in PHC.
We denote the dependent variable, $y_{it}$, as whether individual $i$ has PHC in year $t$. We are interested in the effects of “Lifetime Health Cover” (LHC) policy, “Medicare Levy Surcharge” (MLS), and the state persistence $(y_{it-1})$ on one's purchasing PHC.
The age dummy variable $\text{Above30}_{it}$ indicates if individual $i$ is 30 years old or above in year $t$, namely, $\text{Above30}_{it}:=\mathds{1} \{\text{Age}_{it}\geq30\}$. Following the insight of the (sharp) “regression discontinuity” design, its coefficient captures the effects of Australia's LHC policy introduced in 2000 to encourage the uptake of PHC. Loosely speaking, the LHC states that if an individual has not taken out and maintained PHC from the year she turns 31, she will pay a 2% LHC loading on top of her premium for every year she is aged over 30 if she decides to take out PHC later in life. If LHC is a strong incentive, we would expect a significant “jump" in the PHC enrollment rate at this age.
The MLS is a levy paid by Australian taxpayers who do not have PHC and earn above a certain income threshold. In the sample years of our data, MLS rates remain unchanged, while the thresholds have been raised yearly until 2014. It is worth noting that the 2014 rise in MLS thresholds was only 50% of previous years, and the thresholds have remained at the same level until 2022. The time dummy $D_{2014,t}$ is included in model ((ref) ) to examine whether this change in MLS policy would affect people's willingness to purchase PHC. Note that honore-k's ( honore-k) estimators do not allow either age ($\text{Age}_{it}$) or fixed time effects ($D_{2014,t}$) since they do not have overlapping supports across time.
$I_{it}$ represents standardized annual household disposable income using the entire sample in the survey, which serves as the continuous regressor with rich enough support required for point identification (by Assumptions A2 and A3). The standardization is performed before dropping missing data by subtracting the sample mean from each individual value and then dividing the difference by the standard deviation.
We also include a (location) dummy variable $\text{GCC}_{it}$ that indicates whether individual $i$ lives in a major city/greater capital city in year $t$. This variable is included to control the accessibility to private hospital services. Tables (ref) and (ref) provide definitions and summary statistics of all aforementioned variables, respectively. Note that some observations are excluded due to missing information in other variables, so in Table (ref), $I_{it}$ does not have an exact zero mean and unit standard deviation.
With all these covariates, we specify our empirical model as follows:
where $\epsilon _{it}$ and $\alpha _{i}$ are, respectively, the usual idiosyncratic error and unobserved heterogeneity in fixed effects panel data models.
In our analysis, we restrict the coefficient on $I_{it}$ to be 1, following the same convention for scale normalization as in Section (ref). This choice warrants justification; that is, household income enters the model with a significant positive coefficient, as implicitly required by Assumption A5'. We provide the following rationale for this based on common sense and evidence from exploratory regression. Practitioners seeking to justify normalizing the coefficient on $z_{it}$ to 1 can adopt similar argumentation method.
Firstly, in Australia, Medicare provides free access to public hospitals, and Medicare patients in public hospitals receive free treatment and accommodations. However, people can purchase private hospital insurance to cover faster and more premium services. Taking up or maintaining private insurance coverage requires a household to have sufficient disposable income. Besides, as income increases, the marginal utility of saving or other consumption may eventually become lower than that of enhanced private health care. In addition, Australia's tax system also gives considerable financial incentives for high-income households to buy private insurance. Therefore, common sense suggests that income should play a positive and significant role in private insurance purchases.
Secondly, we conduct a simple probit regression using one wave of the data and included income as the only regressor. The estimate is positive and significant, with $p$-value smaller than $10^{-15}$. This result holds true across all data waves, confirming our argument. This finding aligns with the results of more in-depth structural analyses in the health economics literature, such as cheng2014measuring.
Note that Assumption A3 can be demanding. To address this concern, we relax Assumption A3 to Assumption A3' for a model closely resembling the current application and demonstrate identification under this relaxed condition, as detailed in Appendix (ref). Additionally, we justify our use of \( I_{it} \) as \( z_{it} \) under this modified condition in Appendix (ref), specifically by showing the kernel density and summary statistics of \( I_{it+1} - I_{it-1} \). For a more detailed discussion on Assumptions A2 and A3 and their roles in this empirical application, we refer interested readers to Appendix (ref).
Let $x_{it}:=(\text{Above30}_{it},\text{Age}_{it},\text{GCC}_{it})$ and $ \beta :=(\beta _{1},\beta _{2},\beta _{3})$. We estimate $\theta :=(\delta ,\gamma ,\beta)$ through maximizing the objective function
where $u_{it}(\vartheta ):=r(y_{it}-y_{it-2})+d(D_{2014,t+1}-D_{2014,t-1})+(x_{it+1}-x_{it-1})^{ \prime }b+(I_{it+1}-I_{it-1})$ and $\vartheta :=\left(d, r,b\right)$. Objective function ((ref)) extends ((ref)) for longer and unbalanced panels in which the number of waves being observed varies across individuals $i$ ($=:T_{i}$). We select the tuning parameter $\sigma_{n}$ using the same approach as described in Section (ref). It is important to note that the distribution of $I_{it}$ exhibits a significantly longer right tail compared to its left tail (skewness=3.78). Consequently, for sufficiently large $\sigma_{n}$, the objective function ((ref)) has a much larger number of observations to use than the objective function
which extends ((ref)) for left tail observations. In fact, in this application, we set $\sigma_{n}=1.478$, which exceeds the absolute value of the lower bound of $I_{it}$ ($=1.401$ as shown in Table (ref)), thereby effectively using only objective function ((ref)) and observations satisfying $\{I_{it}>\sigma_{n}\}$. Previous versions of this paper explored smaller values of $\sigma_n$ that allowed for the inclusion of left-tail observations (i.e., $\{I_{it}<-\sigma_{n}\}$), yielding similar results.
By construction, procedure ((ref)) 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 honore-k's (honore-k) estimators. As in Section (ref), we use $n^{*}$ to denote the number of “switchers”.
We choose $\sigma _{n}=c\cdot \widehat{\text{std}(I_{it})}\sqrt{\log n^{*}/2.95}$ with $ c=1.0$ and $1.1$ to implement our MS estimation and report the results in Table (ref). We provide summary statistics for the sub-sample of switchers with $I_{it}>\sigma_n$ in Table (ref) of Appendix (ref). We also conducted estimations using $\sigma_n$ with $c=0.5, 0.7,$ and $0.9$. While these results show patterns similar to Table (ref), they highlight the bias-variance trade-off inherent in choosing the tuning parameter, a common challenge in semiparametric methods. These additional results and their discussion are included in Appendix (ref).
In addition to the estimates of $\theta$, we also try calculating the 90% and 95% confidence intervals (CIs) for $\theta $ using the $m$-out-of-$n$ bootstrapping. Here we sample $n$ individuals (clusters) to create the bootstrap sample. The main difficulty in implementing this (or alternative sampling-based) method is that Theorem (ref) 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 LeePun2006 to solve this problem; that is, assume $\hat{\theta}_n$ has convergence rate of $ n^{\lambda }$ and obtain an estimate $\hat{\lambda}$ of $\lambda $ using a double $m$-out-of-$n$ bootstrapping procedure with two bootstrap sample sizes $m_{1}=n^{\rho _{1}}$ and $m_{2}=n^{\rho _{2}}$ for $\rho _{1},\rho _{2}\in (0,1)$. The 90% and 95% CI reported in Table (ref) are calculated with $B=500$ bootstrap replications, $m=n^{7/8}$, and $\hat{ \lambda}=0.309$ (obtained with $\rho _{1}=6/7$ and $\rho _{2}=7/8$).
We can see from Table (ref) that the estimation results are similar for the two tuning parameters. Therefore, the following discussion on the empirical results will be mainly based on the estimates obtained with $c=1.0$. The insignificant coefficient indicates that living in GCC may not affect people's willingness to buy PHC. The significant positive coefficient on $y_{it-1}$ demonstrates the strong state persistence of PHC, which explains why we can only observe a small percentage of switchers in the data. Surprisingly, people's decision to buy PHC is hardly influenced by age. For the two policy variables, the large positive coefficient on $\text{Above30} _{it}$ confirms that the LHC policy is a strong incentive for people to buy PHC, while the change in MLS income threshold does not exhibit a strong impact represented by the coefficient on $D_{2014,t}$. An intuitive explanation for the latter is that although MLS promoted PHC purchases when it was introduced in 1997--1998, the subsequent adjustments of its income threshold only affected a small group of people whose incomes were near the threshold.
We end this section with some remarks. First, our approach is more suitable for data with a relatively large proportion of “switchers” which make up the effective sample for the estimator. Second, to implement our method, the model should have a continuous covariate with large support and ideally weak dependence on other included covariates. Third, in the absence of knowledge (or at least a good estimate) of the free-varying variable's tail probabilities, the asymptotics of our estimator derived in Section (ref) cannot provide a “rule of thumb” for choosing optimal tuning parameter $\sigma_{n}$. Perhaps a practical way is to try different $\sigma_{n}$'s, use LeePun2006's ( 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 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.
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 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 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 honore-k. This property makes our methods more practical for models with many covariates or important covariates whose support may not overlap over time. Our identification arguments motivate a conditional maximum score estimator that is proven to be consistent and with the convergence rate independent of the model dimension. However, the asymptotic distribution of the proposed estimator is non-Gaussian, in line with well-established literature on cube-root asymptotics. We suggest valid bootstrap methods for conducting statistical inference. The results of a Monte Carlo study demonstrate that our estimator performs adequately in finite samples. Lastly, we use the HILDA data to investigate the demand for private hospital insurance in Australia.
This paper leaves some open questions for future research. For instance, although we suggest several theoretically feasible bootstrap inference methods in Section (ref), 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 horowitz's (horowitz) and andrews's (andrews) methods to smooth the sample objective function.
\nocite{HILDA,HILDA2}