Tsukuba Economics Working Papers
No. 2019-004

Simple Semiparametric Estimation of Ordered Response Models: with an Application
to the Interdependence Duration Models
byRuixuan Liu and Zhengfei Yu
November 2019

UNIVERSITY OF TSUKUBA
Department of Economics
1-1-1 Tennodai
Tsukuba, Ibaraki 305-8571
JAPANSimple Semiparametric Estimation of Ordered
Response Models: With an Application to theInterdependent Durations Model
Ruixuan Liu and Zhengfei Yu
Emory University and University of Tsukuba

Abstract. We propose two simple semiparametric estimation methods for ordered response models with an unknown error distribution. The proposed methods do not require
users to choose any tuning parameter and they automatically incorporate the monotonicity
restriction of the unknown distribution function. Fixing finite dimensional parameters in
the model, we construct nonparametric maximum likelihood estimates (NPMLE) for theerror distribution based on the related binary choice data or the entire ordered response
data. We then obtain estimates for finite dimensional parameters based on moment conditions given the estimated distribution function. Our semiparametric approaches deliverroot-n consistent and asymptotically normal estimators of the regression coefficient andthreshold parameter. We also develop valid bootstrap procedures for inference. We apply our methods to the interdependent durations model in Honoré and de Paula (2010),
where the social interaction effect is directly related to the threshold parameter in thecorresponding ordered response model. The advantages of our methods are borne out insimulation studies and a real data application to the joint retirement decision of married
couples.

Key words: Ordered Response Models, Semiparametric Estimation, Shape
Restriction, NPMLE, Social Interaction, Interdependent Durations
JEL Classification: C14, C25, C35, C41

Corresponding Address: Zhengfei Yu, Faculty of Humanities and Social Sciences, University of Tsukuba,
1-1-1, Tennodai, Tsukuba, Ibaraki, Japan, 305-8571. E-mail: yu.zhengfei.gn@u.tsukuba.ac.jp.
We would like to thank Jason Abrevaya, Bo Honore, Hidehiko Ichimura, Hiroyuki Kasahara, Toru Kitagawa,
Brendan Kline, Zhongjian Lin, Oliver Linton, Essie Maasoumi, Taisuke Otsu, Aureo de Paula, Adam
Rosen, Katsumi Shimotsu, Haiqing Xu, and seminar participants at UT Austin, SUNY Albany, as well as
participants at the 2018 North American Econometric Society Meetings, the 2018 CEME Conference on
Inference in Nonstandard Problems, and Workshop on Advances in Econometrics 2018 for their insightfulcomments and suggestions. Yu gratefully acknowledges the support of JSPS KAKENHI Grant Number
17K13713. The usual disclaimer applies.
12

1. Introduction
We consider the following ordered response model in which the discrete dependent variable Yi is defined by the threshold-crossing rule given covariates Xi , a latent error term εi ,an unknown threshold parameter α0 , and a regression coefficient β0 :0
1 if εi ≤ Xi β0 ,(1.1)
Yi = 2 if Xi0 β0 < εi ≤ Xi0 β0 + α0 ,

3 if εi > Xi0 β0 + α0 ,for i = 1, · · · , n. We maintain the independence assumption between X and ε and let F0 (·)
be the true unknown distribution function of the latent error ε throughout the paper. Given
independent and identically distributed (i.i.d.) observations of (Yi , Xi )ni=1 , the likelihood
function takes the following form:
o
n
∆3i
∆2i
0
n
0
∆1i
0
0
,
(1.2) Ln (α, β, F ) = Πi=1 F (Xi β) [F (Xi β + α) − F (Xi β)] [1 − F (Xi β + α)]
where the indicators ∆1i , ∆2i , and ∆3i are defined by ∆1i = I{Yi = 1}, ∆2i = I{Yi = 2},
and ∆3i = 1 − ∆1i − ∆2i for i = 1, · · · , n.
The ordered response model dates back to Aitchison and Silvey (1957) where the error
distribution F is parameterized and is widely used to characterize ordered categorical outcome in economics, such as consumers’ demand of differentiated products (Prescott andVisscher, 1977; Shaked and Sutton, 1982), symmetric entry games (Bresnahan and Reiss, 1991), schooling choices (Cameron and Heckman, 1998; Cunha, Heckman, and Navarro,2007), credit/liquidity constraints (Attanasio, Koujianou, and Kyriazidou, 2008; Hajivassiliou and Ioannides, 2007) and discrete time (or interval-censored) duration models (Ridder,
1990; Manski and Tamer, 2002). We refer readers to Greene and Hensher (2010) for acomprehensive review. However, the fully parametric procedure leads to an inconsistent estimate and misleading inference if the parametric model of the error distribution is
misspecified. Flexible semiparametric estimation of (1.1) has been studied by Lee (1992),
Melenberg and Van Soest (1996), Klein and Sherman (2002), Lewbel (2000, 2002), Chen
and Khan (2003) and Coppejans (2007), allowing for an arbitrary error distribution. This
literature can be roughly divided into two categories. The first branch employs either kernelor sieve based nonparametric estimation of the functional nuisance component as in Klein
and Sherman (2002), Lewbel (2002), Chen and Khan (2003) and Coppejans (2007). The

3

implication is that the user has to choose a tuning parameter 1, such as the bandwidth inkernel smoothing or the number of sieve basis functions, and there is no clear answer about
the optimal choice in this context. Inevitably, this requires a considerable amount of intervention and judgment on the part of practitioners. The second approach, which estimates
the finite dimensional parameter (α0 , β0 ) (but not F0 ) without tuning parameters, makes
use of the maximum score type estimation (Lee (1992)). 2 Only the consistency result isavailable for Lee’s (1992) estimator and it is expected to have a non-standard limiting distribution with cubic-root convergence rate which also complicates the inference in practice,not to mention that the maximum score type estimator is very hard to compute. Part of
the aforementioned issues can be alleviated by smoothing the sample criterion function, as
done by Melenberg and Van Soest (1996), which again introduces the kernel bandwidth.
Also, the convergence rate of the smoothed maximum score estimator is slower than thestandard root-n rate (Horowitz, 1992, 2009).
In this paper, we propose two simple semiparametric estimation methods for ordered
response models that are fully automatic and free of any tuning parameter. The resulting estimators of the slope coefficient β0 and the threshold α0 are root-n consistent andasymptotically normal. The first approach, which consists of two stages, starts with thelikelihood function for related binary choice data (∆1i , Xi ) for i = 1, · · · , n:
o
n
1−∆1i
0
n
0
∆1i
,
(1.3)
L1n (β, F ) = Πi=1 F (Xi β) [1 − F (Xi β)]to get our estimated distribution function F̂n (·; β) for any given β following Cosslett (1983).
This is the nonparametric maximum likelihood estimator (NPMLE) in the sense of Kiefer
and Wolfowitz (1956) for the binary choice data. In order to differentiate from the NPMLE
in our second approach, we will refer F̂n as the isotonic estimator throughout the paper.
We then estimate the regression coefficient and threshold parameter sequentially by using
certain moment conditions (or solving estimation equations). In fact, the first stage isadapted from the tuning-parameter-free method in Groeneboom and Hendrickx (2018) for
binary choice models, whereas in the second stage, we obtain the estimated threshold froma simple moment condition concerning the binary choice data (∆3i , Xi ) for i = 1, · · · , n.
Our second approach directly maximizes the full likelihood in (1.2) for the ordered response
1In Lewbel (2000), under an additional independence assumption between the special regressor and othercovariates (see Assumption A.5’ [p.157] of Lewbel (2000)), one can apply the ordered data estimator(Lewbel and Schennach, 2007) without tuning parameter. However, in general, Lewbel (2000)’s estimatorneeds a kernel or sieve type estimator of the conditional density in its first stage.
2Another possibility is a two-stage rank estimator that combines Cavanagh and Sherman (1998) and Chen
(2002). Because the original focus of these papers are not on the ordered response model, we discuss theapproach in Section 3.3. We thank an anonymous associate editor for suggesting this method in an earlier
submission of the paper to another journal.

4data to obtain the NPMLE F̃n (·; α, β) for any given (α, β). Thereafter, we estimate theregression coefficient and threshold jointly by using the same moment conditions as in our
two-stage approach. Throughout this paper, we name the first approach the (Isotonic)
two-stage estimation and the second one the (NPMLE-based) joint estimation.
Our estimation approaches have three main appealing features. First, both methods arefree from any tuning parameter, as opposed to Klein and Sherman (2002), Lewbel (2002),
Coppejans (2007), and Melenberg and Van Soest (1996). This is because we estimate theerror distribution F in (1.3) by a well-defined NPMLE (either using the binary choice dataor the ordered response data), which only exploits the shape (monotonicity) restriction. As
a result, the estimator F̂n or F̃n does not rely on any kernel smoothing or sieve penalization.
Second, both estimators of the latent distribution F obtained from our methods are automatically non-decreasing piece-wise constant functions by construction. In contrast, thekernel-based approach in Klein and Sherman (2002) may not yield a monotonic estimateof the latent distribution3 and the sieve estimator in Coppejans (2007) has to incur additional computation costs by restricting spline coefficients to accommodate monotonicity.
Finally, our approach is easy to implement. The isotonic estimator F̂n (·) is computed using
the pool-adjacent-violators algorithm (PAVA) (see Groeneboom and Jongbloed (2014) for
details), whereas for the NPMLE F̃n (·), we adapt the fast hybrid approach in Wellner andZhan (1997) that combines both the E-M algorithm and the iterative greatest convex minorant algorithm. Our two-stage estimation is particularly attractive from the computational
point of view, in the sense that for the given (β̂n , F̂n (·; β̂n )), the estimating equation for thethreshold parameter α is monotonic. In comparison, the second stage of Klein and Sherman(2002) involves repeating a grid search many times in order to obtain the estimator of thethreshold parameter.4
Our interest in (1.1) stems from a recent application of the interdependent durations
model proposed by Honoré and de Paula (2010). The study of social interaction effects hasreceived considerable attention in economics literature and we refer readers to Durlauf andIoannides (2010) for fast expanding results based primarily on discrete choice models. The
seminal papers of de Paula (2009), Honoré and de Paula (2010) bring strategic interactions
explicitly into the optimal timing decisions of economic agents and greatly extend the scopeof prior static discrete choices models by incorporating the time dimension. It is evident inde Paula (2009) and Honoré and de Paula (2010) that the social interaction effect captures
notable externality and explains interesting timing coordination behaviors. Furthermore,

3Figure 3 in Section 5.2 plots the estimated CDF of the latent error in a real data example.
4We refer readers to Stewart (2005) about detailed comparisons of existing semi-nonparametric estimatorsfrom their computational aspect.

5

the scalar representing the social interaction effect in Honoré and de Paula (2010) is directly related to the threshold parameter in the resulting ordered response model, which
motivates us to propose the current flexible semiparametric estimation methods. Both theMonte Carlo simulation and a real data application demonstrate a good performance of our
proposal. We believe that our work makes the Honoré-de Paula model a viable benchmark
for analyzing interdependent durations data. Other applicable examples in economics include Park and Smith (2008), who study late rushes in market entry after a pioneer develops
a new product; Katz and Shapiro (1986), who examine technology adoption; Hong, Kubik,
and Stein (2004), who model stock market participation; de Paula (2009), who observesgroup desertion in the Union Army after the American Civil War; and Honoré and de Paula(2018), who look at the joint retirement decisions of married couples.
We contribute to the literature in several ways. First, the proposed estimation methods complement the existing semiparametric approaches in the sense they free appliedresearchers from choosing any tuning parameter. Meanwhile, unlike Lee (1992) where only
the consistency result is available for the maximum score type estimator, we establish theroot-n consistency and asymptotic normality for our estimators of the finite dimensionalparameters. Second, from a theoretical perspective, our work contributes to the literature of(two-stage) semiparametric estimation and inference that involves shape-restricted components. In the seminal works of Newey (1994) and Chen, Linton, and Van Keilegom (2003),
general theorems are presented for semiparametric estimators involving some first-stage
nonparametric estimation, in which the linear representation for the directional derivativeof the nonparametric estimate is maintained as a high-level assumption and verified forsieve or kernel-type estimators under sufficient smoothness restrictions. A distinction of
our paper is that our nonparametric estimators are only piece-wise constant functions withrandom jump locations determined by the data. In fact, the crux of our theoretical investigation is to prove that certain linear functional (or the directional derivative) of the
shape-restricted estimate is asymptotically normal, for which we adapt the recent breakthrough made by Groeneboom and Hendrickx (2018) for our two-stage estimation. A close
examination reveals that the proof in Groeneboom and Hendrickx (2018) regarding theregression coefficient β is easier, as they can utilize an orthogonal (to the nuisance tangentset) score function to account for the estimation effect of the error distribution implicitly.
This orthogonal direction is well-known for single-index models (Ichimura (1993); Klein andSpady (1993)). In contrast, we have to explicitly characterize the influence of estimating
the distribution through its linear functional in our two-stage estimation. Third, our joint
estimation provides a new example of utilizing the NPMLE in semiparametric models. It
substantially goes beyond Groeneboom and Hendrickx (2018). In our model, the NPMLE making use of information in all three categories lacks an explicit characterization in6contrast to the binary choice case considered by Groeneboom and Hendrickx (2018). As
a result, determining the asymptotic behavior of its linear functional becomes much morechallenging. Aiming at those challenges, our proofs that make novel use of the empirical process theory and the characterization of NPMLE for the “interval censoring, case 2
model”5 (Geskus and Groeneboom (1996, 1997, 1999) and Van de Geer (1997)) are alsoof independent interest. Therefore, this paper also contributes to the literature of shaperestricted estimation and inference for semiparametric models (Cosslett, 1983; Groeneboom
and Wellner, 1992; Matzkin, 1991, 1993; Banerjee, Mukherjee, and Mishra, 2009; Groeneboom, Jongbloed, and Witte, 2010; Groeneboom and Hendrickx, 2017, 2018). Fourth, we
propose bootstrap procedures with general exchangeable weights to facilitate inference, as
the asymptotic variances of our estimators are quite involved. The bootstrap consistency
does not trivially follow from Van Der Vaart and Wellner (1996) or Cheng and Huang(2010), since the bootstrap is known to fail for the point-wise distribution of the isotonic
estimator or NPMLE (Abrevaya and Huang (2005)). Heuristically speaking, the bootstrap
is valid for the finite dimensional parameter in our model, because the influence of the isotonic estimator or the NPMLE only involves certain linear functional (Huang and Wellner(1995),Groeneboom and Hendrickx (2017)) whose asymptotic behavior can be mimickedby the bootstrap. Last but not least, we demonstrate the usefulness of the Honoré-dePaula model by taking it to the real data based on our proposed estimation and inference
methods. Indeed, it shows the interaction effect between couples is large enough to counter
their age difference in explaining the joint retirement behavior. Considering the continued
interest in modeling social interactions among economists and the level of sophistication of
the Honoré-de Paula model, we believe our empirical illustration is worthwhile.
The rest of our paper is organized as follows. We briefly discuss the connection to theinterdependent durations model in Honoré and de Paula (2010) in Section 2. In Section 3,
we introduce the notation and propose two simple semiparametric estimation methods. In
Section 4, we derive the asymptotic properties of our estimators for the finite dimensionalparameters, proving their consistency and asymptotic normality. Given the complicatedinfluence functions, we also develop novel bootstrap procedures so that one can easily
formulate the bootstrap confidence sets. Section 5 conducts simulation studies to evaluate
the finite-sample properties of the estimators and also illustrates the proposed methods
using a real dataset. The final section concludes. Proofs of main theorems are in Appendix
A, whereas proofs of technical lemmas are delegated to the supplemental note.

5To clarify the comparisons with Groeneboom and Hendrickx (2018), the binary choice model there is also

known as the “interval censoring, case 1 model”; see Groeneboom and Wellner (1992).

7

2. The Interdependent Durations Model
Interactions among economic agents play key roles in characterizing simultaneity of multiple durations models. This endogenous interaction effect generates interesting synchronization for numerous economic phenomena. Honoré and de Paula (2010) study a gametheoretic model in which two players respectively decide (T1 , T2 ) as the time of switchingfrom an initial activity to an alternative activity. The utility flow of the alternative activity
for one player depends on whether the other player has switched or not, which causes anendogenous interaction effect. The equilibrium of two duration variables (T1 , T2 ) with theendogenous interaction effect are characterized by(2.1)

T1 ≡ inf {t1 : φ (X1 ) Λ (t1 ) exp [α∗ I {T2 ≤ t1 }] ≥ 1 } ,T2 ≡ inf {t2 : φ (X2 ) Λ (t2 ) exp [α∗ I {T1 ≤ t2 }] ≥ 2 } ,where the unknown scalar α∗ captures the social interaction effect. The common Λ(·) isa deterministic trend function. Each player has his/her own covariate Xj for j = 1, 2 and(1 , 2 ) are random utility flows from the initial state. Throughout the paper, we focus on
the case where an individual’s switching is complementary to each other, so the interaction
effect α∗ is positive.
In our empirical application, we utilize the model (2.1) to analyze the joint retirement decision of married couples so (T1 , T2 ) represents the retirement timing choices of the couple.
Understanding the retirement mechanism helps to guide the optimal design of employerprovided and government benefit programs (Gustman and Steinmeier (2000), Gustman andSteinmeier (2004), An, Christensen, and Gupta (2004), Honoré and de Paula (2018)). Our
focus on the married couple stems from the fact that most people approaching retirement
age are married and a significant portion of them choose to retire at the same time. Empirical studies documenting the joint retirement of couples abound based on different datasets(Honoré and de Paula (2018)). Standard bivariate mixed proportional hazards models do
not apply to this scenario, because the simultaneous failure (when Pr{T1 = T2 } > 0) is
ruled out from the very beginning in those models. In the current context, the complementary interaction effect is associated with the pension accrual profiles or the complementarity
in leisure time; see Gustman and Steinmeier (2000) for a more detailed discussion.
Despite the sophisticated game structure and presence of multiple equilibria in (2.1),
Honoré and de Paula (2010) prove the identifiability of all model primitives. Moreover,they point out a close connection with the corresponding ordered responses model, which
directly suggests a natural semiparametric estimation route for the interaction effect. We

8start to define the auxiliary discrete dependent variable:
1 if T1 < T2 ,

Y = 2 if T1 = T2 ,
3 otherwise.
Then we parameterize the covariate effect by setting φ(xj ) = x0j β0 for j = 1, 2, as in Honoré
and de Paula (2010). Therefore, one gets
(2.2)

P {Y = 1|X1 , X2 } = H ((X1 − X2 )0 β0 − α∗ ) ,P {Y = 2|X1 , X2 } = H ((X1 − X2 )0 β0 + α∗ ) − H ((X1 − X2 )0 β0 − α∗ ) ,P {Y = 3|X1 , X2 } = 1 − H ((X1 − X2 )0 β0 + α∗ ) ,where H(w) = P {log 1 − log 2 ≤ w}. Hence, it is straightforward to observe that (2.2)
corresponds to the semiparametric ordered response model by taking F0 (·) = H(· − α∗ ),X = (X1 − X2 ), and α0 = 2 × α∗ .
Remark 2.1. Note that the ordered response model does not capture all model primitives
in the original interdependent durations model, in particular the common deterministic
trend function Λ(·) is missing. However, the constructive identification in Theorem 3 ofHonoré and de Paula (2010) suggests a straightforward nonparametric estimate for ln Λ(t),combined with the key insight of Horowitz (1996). To elaborate on the proposal, leth(t1 , t2 ; x1 , x2 ) ≡ Pr{T1 ≤ t1 , T2 > t2 |X1 = x1 , X2 = x2 }fort1 < t2 ,then the proof of Theorem 3 in Honoré and de Paula (2010) states6 that
(2.3)

∂ ln Λ(t1 )∂h/∂t1= β0k.
∂t1∂h/∂x1kThus, one can adopt the estimator defined by equation (2.4) in Horowitz (1996) by plugging
in some nonparametric kernel estimator for the partial derivative of h(t1 , t2 ; x1 , x2 ). Given
that the theoretical properties follow largely from Horowitz (1996) and Λ(t) does not directly
enter the resulting ordered response model, we will not focus on the issue.

3. Simple Semiparametric Estimation
Throughout the paper, we work with the i.i.d. data (Yi , Xi ) for i = 1, ..., n. It isconvenient to introduce the indicators ∆1i , ∆2i , and ∆3i defined by ∆1i = I{Yi = 1} and6Note that in our setup, the covariate effects are parameterized by φ(x) = exp(x0 β ).
∂φ/∂x1kφ(x1 )

Hence, the ratio
0
defined in Theorem 3 of Honoré and de Paula (2010) is always equal to β0k in our specification.

9∆2i = I{Yi = 2}, and ∆3i = 1−∆1i −∆2i for i = 1, · · · , n. Let K denote the dimensionalityof covariates X and write β0 ≡ (β01 , β02 , ..., β0K )0 . Note that the regression coefficient β0 is
only identified up to some scale normalization for an unspecified F0 (see Klein and Sherman(2002)). Without loss of generalization, we normalize β01 = 1.
We use the standard empirical process notations as follows. For a function f (·) of aR
random vector Z = (Y, X) that follows distribution P , we let P f = f (z)dP (z), Pn f =P
n−1 ni=1 f (Zi ), and Gn f = n1/2 (Pn − P ) f . Function f can be replaced by a randomRfunction z 7→ fˆn (z; Z1 , · · · , Zn ). Therefore, P fˆn = f (z; Z1 , · · · , Zn )dP (z), Pn fˆn =Pn−1 ni=1 f (Zi ; Z1 , · · · , Zn ) and Gn fˆn = n1/2 (Pn − P ) fˆn ; see Nan and Wellner (2013). Let
η = (θ, F (·; θ)) be the unknown parameter containing both finite dimensional parameterθ ≡ (α, β 0 )0 and the distribution function F . Furthermore, we consider α ∈ A, β ∈ B, andF ∈ F, where A ⊂ R+ , B ⊂ RK , and F is the class of distribution functions. The distance
between two parameter values (η1 , η2 ) is defined in terms of the following metricd (η1 , η2 ) = |θ1 − θ2 | + kF1 (·; θ1 ) − F2 (·; θ2 )kwhere |·| is the standard Euclidean distance, and k·k is some norm for the class of distribution functions. We work with the L∞ −norm related to our two-stage estimation and theL2 − norm regarding our joint estimation for technical convenience.

3.1. Two-stage Semiparametric Estimation
Our two-stage estimation procedure is inspired by Klein and Sherman (2002) and Lewbel(2002) in the sense that we obtain estimated distribution function from the related binary
choice data (∆1i , Xi ) and (∆3i , Xi ). Unlike Klein and Sherman (2002) who resort to thekernel estimator in Klein and Spady (1993) or Lewbel (2002) who requires a preliminary
nonparametric estimation of certain conditional mean function and integration over covariates’ values, we adapt the important breakthrough by Groeneboom and Hendrickx (2018)
to analyze ordered response models. As it will become self-evident in the sequel, our procedure delivers root-n consistent and asymptotically normal estimators of the regression
coefficient β0 and threshold parameter α0 , which does not require tuning parameters.
We now describe the two-stage semiparametric estimation for the ordered response model:
Stage 1(i). For any β, we compute the NPMLE for F (·) based on the binary choice data:(3.1)

F̂n (·; β) = arg max
F

n
X
i=1

[∆1i log F (Xi0 β) + (1 − ∆1i ) log(1 − F (Xi0 β))] .

10Stage 1(ii). Given F̂n (·; β), our estimator β̂n for the regression coefficient is the zerocrossing point of the estimation equation:ni
1X hXi ∆1i − F̂n (Xi0 β̂n ; β̂n ) = 0.
(3.2)
n i=1
Stage 2. Given β̂n and F̂n (·; β̂n ), we estimate α0 by α̂n , which is the zero-crossing pointof the estimation equation Ψn α̂n , β̂n , F̂n (·; β̂n ) = 0, where
(3.3)

Ψn



n
i
1 Xh
α, β̂n , F̂n (·; β̂n ) =1 − ∆3i − F̂n (Xi0 β̂n + α; β̂n ) .
n i=1



Here we provide a heuristic discussion of each step. The NPMLE in Stage 1 and its characterization date back to Ayer, Brunk, Ewing, Reid, and Silverman (1955) in analyzing
the current status data (see Groeneboom and Wellner (1992)). The corresponding optimization problem is well-defined and it generates a piece-wise constant function F̂n (·; β),which can be characterized as follows: fixing the parameter β, we consider the values of
(β)
(β)
(β)
(β)
U1 = X10 β, · · · , Un = Xn0 β. Let U(1) ≤ · · · ≤ U(n) be the order statistics with corre(β)

sponding indicators ∆1,(i) for i = 1, · · · , n. Thereafter, F̂n (·; β) is equal to the left derivativeof the convex minorant of a cumulative sum diagram consisting of the points (0, 0) and
!
i
X
(β)
i,
∆1,(j)
j=1

for i = 1, · · · , n; see Groeneboom and Hendrickx (2018). Within the context of binary
choice models, it is utilized by Cosslett (1983) to define the tuning-parameter-free profile
likelihood estimator. However, only consistency results 7 are available for Cosslett’s estimator given the challenge that the estimated error distribution is neither linear nor smooth.
The key to develop a root-n consistent and asymptotic normal estimator for β0 while maintaining the tuning-parameter-free feature is the Z-estimator adapted from Groeneboom andHendrickx (2018). Modulo the estimated latent distribution function, one makes use of the
population-level moment condition(3.4)

E[X(∆1 − F0 (X 0 β0 ))] = 0,

7The general result in Tanaka (2008) implies the convergence rate of Cosslett’s estimator is cubic root forthe slope coefficient.

11and plug in the NPMLE F̂n (·; β) in the sample analog.8 In the same spirit, the last step ofour procedure is based on a very simple moment condition:(3.5)

E[(1 − ∆3 − F0 (X 0 β0 + α0 ))] = 0.

Given preliminary estimators F̂n (·; β̂n ) and β̂n , the estimating equation Ψn is monotonewith respect to α, which greatly lowers the computational burden. We emphasize that
it is necessary to use both sets of moment conditions for the sake of consistency. The
naive approach where one uses only the binary choice data (∆2i , Xi ) and then directly
applies Groeneboom and Hendrickx (2018) does not work because the intercept α0 andthe distribution function F0 cannot be separately identified in the binary choice data alone(Ichimura (1993)). We focus on the just-identified case to be coherent with Groeneboom
and Hendrickx (2018). In principle, a GMM type estimation with over-identified moment
conditions could be developed.
As F̂n (·; β̂n ) is a piece-wise constant function, the estimating equations in our procedure
may not hold exactly. Therefore, we adopt the following definition from Groeneboom andHendrickx (2018):
Definition 3.1 (Zero-crossing). We say that β∗ is a zero-crossing of a function C : B 7→ R
if each open neighborhood of β∗ contains points β1 , β2 ∈ B such that C(β1 )C(β2 ) ≤ 0. We
say that a function C̃ : B 7→ RK crosses zero at point β∗ if β∗ is a zero-crossing in eachcomponent C̃j for j = 1, · · · , K.

Figure 1 depicts the typical shape of our estimating function Ψn α, β̂n , F̂n (·; β̂n ) . The
data generating process is the interdependent duration model described in Section 5.1. The
true value of α is 1. The joint distributions of (1 , 2 ) are exponential (Panel (a)) and lognormal
 (Panel (b)),respectively. It is obvious that these two depicted estimating functions
Ψn α, β̂n , F̂n (·; β̂n ) are decreasing in α and the zero crossings are very close to the 1; i.e.,they are about 0.980 in Panel (a) and 0.868 in Panel (b).

3.2. Joint Semiparametric Estimation
One might be wondering whether it is possible to develop a similar tuning-parameter-freeestimation approach utilizing the entire ordered response data directly, instead of breaking
it to two sets of binary choice data. Indeed this is feasible and will be termed as the joint
estimation method. The starting point is the following result that guarantees the existence
8The main improvement made by Groeneboom and Hendrickx (2018) over Cosslett (1983) to restore stan-

dard distributional theory for regression coefficient is that one does not need the error’s density functionin the moment condition (3.4). In contrast, one has to handle the error density in the likelihood based
estimation appearing in the score function, whereas the NPMLE F̂n (·; β) itself is not differentiable.

12





Figure 1. The estimating function Ψn α, β̂n , F̂n (·; β̂n ) when the joint
distribution of (1 , 2 ) in the interdependent duration model in Section 2 isexponential or log-normal. α0 = 1, sample size = 500.
(a). Exponential(b). Log-normalof the NPMLE for the ordered response model given any (α, β) in the parameter space.
The proof can be found in Chapter 4 of Groeneboom and Jongbloed (2014).
Lemma 3.1. Given α and β, the NPMLE for F̃n (·; α, β) exists based on the ordered response data; i.e.,F̃n (·; α, β)= arg maxFnX[∆1i log F (Xi0 β) + ∆2i log (F (Xi0 β + α) − F (Xi0 β)) + ∆3i log (1 − F (Xi0 β + α))]i=1

is well defined. Moreover, the NPMLE is a (sub-)distribution function and piece-wise constant with jumps over a subset of {Xi0 β, Xi0 β + α : i = 1, 2, ..., n}.
Now we are ready to present the joint semiparametric method that makes use of information in all three categories to estimate the distribution function and returns the estimates
for the regression coefficient and the threshold parameter simultaneously.
Joint Estimation. For any α and β, we employ the NPMLE for F̃n (·; α, β) based on theordered response data. Given F̃n (·; α, β) from the previous step, we define (α̃n , β̃n ) which
are the zero-crossing points of the estimation equations simultaneously:(3.6)

Φn (α̃n , β̃n ) = 0,

13

where

i 
h
Pn
0
β;
α,
β)
X
∆
−
F̃
(X
i
1i
n
i
i=1
h
i .
Φn (α, β) ≡ Pn1
0
i=1 1 − ∆3i − F̃n (Xi β + α; α, β)
n


1
n

The NPMLE F̃n (·; α, β) can be computed by the iterative convex minorant algorithm in
Groeneboom and Wellner (1992) and Groeneboom and Jongbloed (2014). The number ofmass points is smaller than 2n, because for any i with ∆2i = 0, either Xi0 β +α (if ∆1i = 1) or
Xi0 β (if ∆3i = 1) does not enter the log-likelihood function. Denote the remaining elements(α,β)
in the set {Xi0 β, Xi0 β + α : i = 1, 2, ..., n} as Uj , j = 1, 2, ..., p. Partition the observations
into the following four groups:(α,β)= Xi0 β for some i and ∆1i = 1},(α,β)= Xi0 β for some i and ∆2i = 1},(α,β)= Xi0 β + α for some i and ∆2i = 1},(α,β)= Xi0 β + α for some i and ∆3i = 1}.

I1 = {1 ≤ j ≤ p : Uj

I2l = {1 ≤ j ≤ p : Uj

I2r = {1 ≤ j ≤ p : Uj
I3 = {1 ≤ j ≤ p : Uj

And then set k be a function that maps any index from I2l to I2r for a given observation
(α,β)
(α,β)
= Xi0 β + α, for ∆2i = 1. Let= Xi0 β and Um
i with ∆2i = 1: k(j) = m if Uj(t)(t)
v (t) ≡ (v1 , ..., vp )0 be the output from the t-th iteration, then v (t+1) is the left derivativeof the cumulative sum diagram consisting of the following points:
!
j
j
X
X
(t)
P0 = (0, 0), Pj =Hj (v (t) ),vi Hj (v (t) ) − Gj (v (t) ) , j = 1, ..., p,
i=1

i=1

where

Gj (v) =−vj−1 if j ∈ I1 ,




(v
− v )−1 if j ∈ I ,
k(j)

j

2l

−1



−(vj − vk−1 (j) )if j ∈ I2r ,(1 − vj )−1 if j ∈ I3 ,and
vj−2 if j ∈ I1 ,




(v
−2
if j ∈ I2l ,k(j) − vj )Hj (v) =(vj − vk−1 (j) )−2 if j ∈ I2r ,(1 − vj )−2 if j ∈ I3 .
The initial value can be set as v (0) = (1/p, 2/p, ..., 1)0 , which assigns the same probability
mass on each jump point. The iterative convex minorant algorithm can be implemented

14

using the R package Icens (Gentleman and Vandal, 2018). R package Icens also provides
a function for a faster hybrid algorithm proposed by Wellner and Zhan (1997), which
combines the iterative convex minorant and the E-M algorithm.
Figure 2 plots the estimating functions with respect to α for fixed β = β̃n . The designs
follow the ones presented in Section 5.1, with β0 = (1, β02 , β03 ). Therefore, (3.6) contains
three estimating equations. In Figure 2, the estimate α̃n is the value where the threeestimating functions cross zero, which equals to 0.890 and 1.036 for the exponential and
log-normal design, respectively. Compared with the two-stage estimator, there are mainly
two complications. First, all estimating equations are not guaranteed to be monotone withrespect to α. Second, we need to solve for the joint zero-crossing point for all equations
simultaneously. Nevertheless, as proved in our Theorem 4.3, the zero crossing points exist
with probability 1.

Figure 2. The estimating functions Ψn α, β̂n , F̂n (·; β̂n ) when the joint
distribution of (1 , 2 ) in the interdependent duration model in Section 2 isexponential or log-normal. α0 = 1, sample size =500.
(a). Exponential(b). Log-normal

3.3. Comparisons with Alternative Methods
This section describes some alternative semiparametric estimators for the model (1.1),
and also discusses the general ordered response models with more than three categories.

15

Remark 3.1 (The kernel-based estimator in Klein and Sherman (2002)). In the first stageof Klein and Sherman (2002)’s approach, the regression coefficient β0 is estimated by maximizing the following quasi-likelihood function with respect to β:nnoXτ̂ (Xi ) I{Yi = 1} ln P̂1 (Xi0 β) + I{Yi = 2} ln[P̂2 (Xi0 β) − P̂1 (Xi0 β)] + I{Yi = 3} ln[1 − P̂2 (Xi0 β)] ,i=1

where P̂j (Xi0 β) is the kernel estimator of the conditional probability Pj (Xi0 β) ≡ Pr(Yi ≤0j|Xn i β) forofixed β in Klein and Spady (1993), j = 1, 2, and the trimming function τ̂ (x) =I |x| ≤ ξˆ with ξˆ being certain sample quantile of |Xi |’s. In the second stage, the thresholdparameter α0 is estimated through the shift restriction P2 (Xi0 β0 − α0 ) = P1 (Xi0 β0 ), whichleads to
X
1
α̂ =
(3.7)
(V̂i − Ṽi2 ),I{i ∈ T } i∈Tonwhere V̂i ≡ Xi0 β̂ and Ṽi2 solves P̂2 (Ṽi2 ) = P̂1 (V̂i ) 9 for each i ∈ T , T = V̂i : P̂L ≤ P̂1 (V̂i ) ≤ P̂U
and (P̂L , P̂U ) are determined by the pth and (1−p)th quantile of an a collection of estimatedprobabilities (see Klein and Sherman (2002), p671). Besides the choice of bandwidth forP̂j (Xi0 β), the K-S approach also requires two trimming parameters. Our simulations studies
find that the performance of the threshold parameter estimator α̂ is particularly sensitive tothe trimming parameter p used for constructing the target set T . (See Table 2 in Section
5.1.) Note that the purpose of the trimming is to exclude individual estimators V̂i − Ṽi2 withpoor performance.
Remark 3.2 (The smoothed maximum-score estimator). Horowitz (1992) initially proposed the smoothed maximum-score (SMS) estimator for the binary choice model and Melenberg and Van Soest (1996) have extended it to the ordered response model. The estimatorin Melenberg and Van Soest (1996) can also be viewed as smoothed version of Lee (1992)’s
maximum-score estimator for the ordered response model. Under the median independence
condition, the SMS approach estimates (α0 , β00 )0 by maximizing the smoothed objective function:




n
X
−Xi0 β − α−Xi0 β+ (2I{Yi ≥ 3} − 1) K,
max0
(2I{Yi ≥ 2} − 1) K
α,β
h
h
i=1
where K(v) is an integral kernel function satisfying limv→+∞ K(v) = 1 and limv→−∞ K(v) =0.10 Users need to specify the bandwidth h. The convergence rate of SMS is slower than
9In its actual implementation, Ṽ

i2 is the point for which P̂2 (Ṽi2 ) is closest to P̂1 (V̂i ) over a grid constructedfollowing the procedure given in [p.671-672] of Klein and Sherman (2002).
10Horowitz (1992) adopts K(v) = 0.5 + (105/64)[v − (5/3)v 3 + (7/5)v 5 − (3/7)v 7 ] if |v| ≤ 1; K(v) = 0 if
v < −1 and K(v) = 1 if v > 1.

16

the root-n rate, even with the MSE-optimal bandwidth of order n−1/5 ; see Section 4.3.3 ofHorowitz (2009) for a detailed discussion. Simulation results in Section 5.1 find that SMS
yields larger mean square errors and longer confidence intervals than other methods.

Remark 3.3 (A two-stage rank estimator). An alternative tuning-parameter-free method is
the two-stage rank estimation combining Cavanagh and Sherman (1998) and Chen (2002)11.
The original focus of Chen (2002) is the unknown link function in the transformation model;
however, his method is also applicable to the ordered response model. Specifically, one canestimate β0 following Cavanagh and Sherman (1998) first by
(3.8)

max
β

n
X

Yi Rn (−Xi0 β),

i=1where Rn (−Xi0 β) denotes the rank of −Xi0 β. The second stage adapts the key idea of Chen
(2002) that given P (Yi = 1|Xi ) = F0 (Xi0 β0 ) and P (Yi ≤ 2|Xi ) = F0 (Xi0 β0 + α0 ), one has
E [I{Yi = 1} − I{Yj ≤ 2}|Xi , Xj ] ≥ 0 whenever Xi0 β0 − Xj0 β0 ≥ α0 for i 6= j.
Therefore, a maximum rank correlation estimator for the threshold isn X
X1α̂ = arg max(I{Yi = 1} − I{Yj ≤ 2}) I{Xi0 β − Xj0 β ≥ α}.
α n(n − 1)i=1 j6=i

Our simulation studies demonstrate remarkable performances of this two-stage rank estimator for estimating the finite dimensional parameters, especially when the sample size isrelatively large. However, unlike our two estimators or the Klein-Sherman estimator, therank estimator itself does not provide a nonparametric estimate for the error distribution.
In our empirical application, the error distribution function is still of practical interest,since it captures the (log-)difference of two heterogeneity terms for wives and husbands.

3.4. Models with Four or More Categories
Here we discuss the general ordered response model where the dependent variable cantake more than three values such as in empirical applications of Cameron and Heckman

11This approach credits to an anonymous associate editor in an early submission of the paper to anotherjournal. Given the connection is not obvious, we expand our discussion here.

17(1998) or Klein and Sherman (2002). Formally, the dependent variable is determined by
1 if εi ≤ Xi0 β0 ,






2 if Xi0 β0 < εi ≤ Xi0 β0 + α0,1 ,


(3.9)
Yi = 3 if Xi0 β0 + α0,1 < εi ≤ Xi0 β0 + α0,2 ,


..



.




J + 1 if εi > Xi0 β0 + α0,J−1 ,for i = 1, · · · , n and J ≥ 3, with a set of ordered thresholds (α0,1 , · · · , α0,J−1 ). In thesame spirit of Lewbel (2002), our two-stage estimator in Section 3.1 directly applies tothis setting. In the first stage, the binary data (I(Yi > 1), Xi )ni=1 is employed to obtainβ̂n and F̂n (·; β̂n ) in the same manner as (3.1) and (3.2). And then one utilizes the data
(I(Yi > j + 1), Xi )ni=1 in the second stage to estimate α0,j through the following estimation
equation of αj for j = 1, · · · , J − 1,
(3.10)

Ψj,n



n
i
1 Xh
0
αj , β̂n , F̂n (·; β̂n ) =1 − I(Yi > j + 1) − F̂n (Xi β̂n + αj ; β̂n ) .
n i=1



The large sample properties of resulting estimates (α̂n,1 , · · · , α̂n,J−1 ) are stated in our Corollary (4.2).
Our joint estimator in Section 3.2 is also applicable. From a computational point ofview, if there are more than three categories, only the interval corresponding to the chosen
category and its adjacent ones will be relevant for the computation of the NPMLE; theother intervals can be discarded (Groeneboom (2014), p2093). Therefore, the constructionof the NPMLE is almost the same as the case with three categories. The consistency of
the NPMLE for multiple categories is shown in Schick and Yu (2000). However, the rateof convergence or the asymptotic properties of its linear functionals remain completely unknown. We leave this challenging issue to the future research. Thereafter, we recommend
practitioners use the methods in Klein and Sherman (2002)12 or Coppejans (2007), if efficiency is the main concern. Since our empirical application in Section 5.2 involves three
categories, we will focus on the setup specified by (1.1).

12The K-S estimator is also semiparametrically efficient under additional periodicity restriction on thecovariates; see Section 3.3 of Coppejans (2007).

18

4. Asymptotic Results
This section consists of three subsections. The first two subsections provide root-n consistency and asymptotic normality results for Isotonic two-stage estimator and NPMLEbased joint estimator, respectively. The last subsection provides bootstrap procedures forthe construction of confidence intervals.

4.1. Asymptotic Properties of The Two-stage Estimator
A key challenge in our proof related to the two-stage estimation is to pin down the asymptotic contribution of F̂n (·, β̂n ) to the finite dimensional parameter using the characterizationof the isotonic estimator and empirical process theory. For the slope coefficients β0 , we rely on the recent breakthrough made by Groeneboom and Hendrickx (2018). Building onthat, we need substantial more efforts to determine the influence function for the threshold parameter α0 , as the unknown parameters in the estimating equation are substituted
with the estimates in Groeneboom and Hendrickx (2018). The proof in Groeneboom andHendrickx (2018) regarding the regression coefficient is easier because they can utilize an
orthogonal (to the nuisance tangent set) score function to incorporate the estimation effectof the error distribution implicitly. This orthogonal direction is well-known for single-index
models (Ichimura (1993) and Klein and Spady (1993)), as it involves the conditional meanof covariates X given the true linear index U = X 0 β0 . In contrast, we have to explicitly
characterize the influence of estimating the distribution through its linear functional in our
two-stage estimation.
We introduce additional notations to present our theoretical results. Following Groeneboom and Hendrickx (2018), it is straightforward to observe that the first-stage isotonic
estimator F̂n (·; β) provides an estimate of
n
o Z
(β)
(β)
(4.1)
F0 (u; β) ≡ P ∆1i |Ui = u = F0 (u + x0 (β0 − β))fX|X 0 β (x|X 0 β = u)dx,
In the sequel, we let F0 (u) = F0 (u; β0 ). The density function of the random variable X 0 β
is denoted by g0 (u; β). Denote the true linear index by Ui = Xi0 β0 for i = 1, · · · , n, andlet g0 (u) be the probability density function of the random variable U . The following two
terms appear in the Taylor expansion in our asymptotic analysis:
(4.2)
(4.3)

Vα0 =
Vβ0 =

∂
E[F0 (X 0 β0 + α)] α=α0 ,∂α∂E[F0 (X 0 β + α0 ; β)] β=β0 ,∂β19whereas
(4.4)

h
i
⊗2
Hβ0 = E f0 (X 0 β0 ) {X − E[X|X 0 β0 ]}denotes the Hessian matrix for β̂n in Groeneboom and Hendrickx (2018).
The following regularity conditions are adapted from Ichimura (1993); Klein and Spady(1993); Klein and Sherman (2002); Groeneboom and Hendrickx (2018). The independence
between ε and X in Condition 1 is an standard assumption used by the kernel or sieve based
estimators (Klein and Sherman, 2002; Coppejans, 2007) and the rank estimator describedin Remark 3.3 (Cavanagh and Sherman, 1998; Chen, 2002). In comparison, the (smoothed)maximum score estimator of Lee (1992) or Melenberg and Van Soest (1996) only requires themedian independence, yet it delivers a convergence rate of the finite-dimensional parameterthat is slower than root-n. Root-n consistent estimators that permits general forms ofheteroskedasticity have been proposed by Lewbel (2002) and Chen and Khan (2003) which
again relies on kernel smoothing. The large support of the linear index X 0 β0 in Condition
1 is necessary. The assumption imposed here only requires the support to include thepositive infinity, whereas the support can be bounded away from the negative infinity (sothe linear index could follow exponential or gamma-type distribution). This condition ismade to exclude the non-identifiable example in Manski (1985) and to ensure that thecomponent ψF0 in the influence function of our Theorem 4.2 is well-defined. Considering
the interdependent durations model, the large support condition is also assumed in Honoré
and de Paula (2010) to achieve point identification.
Condition 1. We have i.i.d. data (Yi , Xi ) for i = 1, · · · , n. The covariates X and latenterror ε are independent. The support of X is not contained in any proper linear subspaceKof R . The support of X 0 β0 contains +∞.
Condition 2. The true regression parameter β0 belongs to the interior of B where B is acompact set in RK . The true threshold parameter α0 ∈ A ≡ (αL , αU ), where [αL , αU ] is acompact interval on the positive real line.
Condition 3. The function F0 (·; β) has a strictly positive continuous derivative, whichstays away from zero for all β ∈ B. Moreover, the function F0 (u; β) is twice continuously
differentiable on the interior of the support for β ∈ B.
Condition 4. The density function of random variable X 0 β denoted by g0 (u; β) is continuous and also stays away from zero for all β ∈ B.
Condition 5. The density g0 (u; β) and conditional expectations E[X|X 0 β = u] and E[XX 0 |X 0 β =u] are twice continuously differentiable w.r.t.u. The functions β 7→ g0 (u; β), β 7→ E[X|X 0 β =20u], and β 7→ E[XX 0 |X 0 β = u] are continuous functions for u in the definition domain andall β ∈ B.
R
Condition 6. The matrix Hβ0 is of full rank and the scalar Vα0 6= 0, where Vα0 = f0 (u +α0 )g0 (u)du.
The asymptotic analysis of β̂n and F̂n (·; β) follows directly from Theorem 4.1 on [p.1426]
of Groeneboom and Hendrickx (2018). Specifically, the linear representation for β̂n is asfollows:

 Z
1/2
(4.5)
n Hβ0 β̂n − β0 = {x − E[X|X 0 β0 = x0 β0 ]} {F0 (x0 β0 ) − δ1 } dGn + op (1).
Thus, β̂n is root-n consistent and asymptotically normal with the influence function equalto Hβ−1
φβ0 , where0
φβ0 (Zi ) = (Xi − E[Xi |Ui ])(F0 (Ui ) − ∆1i ).

(4.6)

Regarding the latent error distribution, one gets from Lemma 5.9 on [p.120] of Groeneboom
and Wellner (1992) or Lemma 3.1 on [p.1423] of Groeneboom and Hendrickx (2018) the
following uniform convergence at the cubic root rate (modulo the logarithm factor):(4.7)

k F̂n (u; β̂n ) − F0 (u) k∞ = Op (log n × n−1/3 );The large sample property of α̂n is more complicated and it is our main focus. We
consider the following random map
n
i
1 Xh
0
1 − ∆3i − F̂n (Xi β̂n + α; β̂n )(4.8)
Ψn (α) =n i=1
and its probability limit
Z(4.9)

Ψ(α) =[1 − ∆3 − F0 (X 0 β0 + α)] dP.

First, we prove the uniform convergence by the Glivenko-Cantelli Theorem (Van Der Vaart
and Wellner (1996)), which leads to supα |Ψn (α)−Ψ(α)| → 0. We then verify the population
criterion function has a unique root and the sample analog has a zero-crossing point withprobability tending to 1.
Theorem 4.1 (Consistency of the two-stage estimator). Suppose Conditions (1)-(6) hold.
Then we have: (i). β̂n obtained from equation (3.2) is a consistent estimator of β0 ; (ii).
for all large n, the unique zero-crossing α̂n of Ψn (α) exists with a probability tending to one
and it is a consistent estimator of α0 .

21Because our estimation procedure belongs to the general Z-estimation with bundledparameter and nuisance functional components, we prove the root-n rate and asymptotic
normality of α̂n following the route in Nan and Wellner (2013). Unlike the examples in
Nan and Wellner (2013), which have nuisance nonparametric components either estimablewith root-n rate or subject to certain smoothness restriction, the nonparametric part isestimated utilizing shape restriction in our model. The crux of our proof is to determine
the asymptotic contribution of the estimated latent distribution to the threshold parameter.
Theorem 4.2 (Asymptotic normality of the two-stage estiamtor). Under Conditions (1)(6), the following linear representations hold:√ (4.10)n β̂n − β0 = Gn [ψβ0 (Zi )] + op (1), and√n (α̂n − α0 ) = Vα−1
Gn [(ψ0 + ψF0 + Vβ0 ψβ0 )(Zi )] + op (1),
(4.11)
0where(4.12)

ψβ0 (Zi ) = Hβ−1
(Xi − E[Xi |Ui ])(F0 (Ui ) − ∆1i ),0(4.13)

ψ0 (Zi ) = [1 − F0 (Ui + α0 ) − ∆3i ] ,(4.14)

ψF0 (Zi ) =g0 (Ui − α0 )[∆1 − F0 (Ui )].
g0 (Ui )Therefore,(4.15)(4.16)√ n β̂n − β0 ⇒ N(0, Ωβ0 ), and√n (α̂n − α0 ) ⇒ Vα−1
× N(0, Ωα0 ),0where Ωβ0 = E ψβ0 ψβ0 0 and Ωα0 = E [(ψ0 + ψF0 + Vβ0 ψβ0 )2 ].
Intuitively speaking, the linear representation for the threshold estimator α̂n involves
three parts: the oracle influence function ψ0 given true β0 and F0 , the effect from theestimation of F0 encoded in ψF0 , and the effect from the estimation of β0 collected in ψβ0 .
Despite the fact that we have a closed-form representation here, a brute-force estimation of
the asymptotic variance involves some density function such as g0 (·) or f0 (·). This motivates
us to propose a simple bootstrap approach that circumvents the obstacle in Section 4.3.
Given the linear representation for both α̂n and β̂n , an immediate corollary is the joint
asymptotic normality for θ̂n = (α̂n , β̂n0 )0 as follows. To simplify the presentation, we abuse
the notation somewhat by setting ψα0 ≡ Vα−1
(ψ0 + ψF0 + Vβ0 ψβ0 ).
0Corollary 4.1. Under Conditions (1)-(6), we have√ n θ̂n − θ0 ⇒ N(0, Σ0 ),(4.17)

22

where the asymptotic covariance matrix is Σ0 = E[(ψα0 , ψβ0 0 )0 (ψα0 , ψβ0 0 )].
Remark 4.1. In the seminal works of Newey (1994) and Chen, Linton, and Van Keilegom(2003), general theorems are presented for semiparametric estimators concerning the root-n
consistency and asymptotic normality, which involves some first-stage nonparametric estimation. Maintained as a high-level assumption, Chen, Linton, and Van Keilegom (2003)
assume thati√ h(4.18)n Sn (θ0 ) + Γ(θ0 )[F̂n − F0 ] ⇒ N(0, Σ),for some finite positive definite matrix Σ without specifying how the nonparametric component F0 is estimated in the first stage. In (4.18), Sn (θ0 ) stands for the (normalized) oracle
score function for the parametric part, whereas the direction derivative Γ(θ0 )[F̂n − F0 ] encodes the estimation effect of the nonparametric component. In an earlier celebrated paper,
Newey (1994) also directly assumes the following linear representation holds
n
i
√ h
1 X
√
ψi + op (1),
(4.19)n Γ(θ0 )[F̂n − F0 ] =n i=1
for some zero-mean and square integrable random variables ψi . The illustrative examples
in Newey (1994) or Chen, Linton, and Van Keilegom (2003) are about nonparametriccomponents with sufficient smoothness restrictions and estimated by sieve or kernel-type
estimators. The crux of our theoretical investigation is to show that certain linear functionalof the shape restricted nonparametric estimator is asymptotically normal. The verificationof (4.18) or (4.19) in our context requires considerable more efforts than those examples in
the aforementioned works mainly due to the fact that the isotonic estimator or the NPMLE
is neither smooth nor linear.
Considering models where the ordered responses have more than three categories as
described by the thresholding rules in (3.9), we collect the finite dimensional parameter in
θ0J = (α0,1 , α0,2 ...α0,J−1 , β 0 )0 and denote the two-stage semiparametric estimator in Section
3.4 by θ̂nJ = (α̂n,1 , α̂n,2 ...α̂n,J−1 , β̂n0 )0 with J ≥ 3. Furthermore, we introduce the following
notations:
ψ0,j (Zi ) = [F0 (Ui + α0,j ) − I(Yi ≤ j + 1)] ,∂E[F0 (X 0 β0 + α)] α=α0,j ,Vα0,j =
∂α
ψα0,j ≡ Vα−1
(ψ0,j + ψF0 + Vβ0 ψβ0 ),0,jfor j = 1, 2, ..., J − 1. The next corollary presents the asymptotic normality of θ̂nJ . Its
proof follows from a straightforward modification of our Theorem 4.2 up to some notation
changes, and thus is omitted.

23

Corollary 4.2. Suppose Conditions (1)-(6) hold and Vα0,j 6= 0 for all j. Then we have√  J
(4.20)n θ̂n − θ0J ⇒ N(0, ΣJ0 ),whereΣJ0 = E[(ψα0,1 , ψα0,2 , ..., ψα0,J−1 ψβ0 0 )0 (ψα0,1 , ψα0,2 , ..., ψα0,J−1 , ψβ0 0 )].

4.2. Asymptotic Properties of The Joint Estimator
The asymptotic analysis of our joint estimation method is more involved. There aremainly two differences compared with the results related to the two-stage estimation. First,
both α̃n and β̃n are determined in the second stage estimation together, so the calculationof the Hessian matrix is done with respect to both components; see Lemma S2.11 in the
supplemntal note. More importantly, it is much more challenging to pin down the influencefunction capturing the effect of the NPMLE F̃n (expressed via some linear functional ofF̃n ). Unlike the binary choice case where the isotonic estimator can be characterized as aleft-continuous slope of the greatest convex minorant of certain cusum diagram, such an
explicit interpretation is lacking for NPMLE using entire ordered response data. Thus,
we seek an alternative characterization that builds on a sequence of research by Van deGeer (1997) and Geskus and Groeneboom (1996, 1997, 1999) for the interval censored data(case 2) and combines it with an analytic argument using empirical process theory as in
Groeneboom and Hendrickx (2018).
First, we set up necessary notations. Denote the probability limit of the NPMLE as(4.21)F0 (u; α, β) ≡ E[∆1 |X 0 β + α = u] = E[∆1 |U = u − α].

Consider the following matrix
H(α, β) ≡!
∂E[−X ∂αF (X 0 β; α, β)]E[−X ∂β∂ 0 F (X 0 β; α, β)],∂E[− ∂αF (X 0 β + α; α, β)] E[− ∂β∂ 0 F (X 0 β + α; α, β)]which leads to the Hessian matrix of our estimator when evaluated at the true parametervalues (α0 , β0 ):!
E[(X − E[X|X 0 β0 ])f0 (X 0 β0 )]E[(X − E[X|X 0 β0 ])⊗2 f0 (X 0 β0 )]H0 ≡ H(α0 , β0 ) = −.
E[f0 (X 0 β0 + α0 )]E[(X − E[X|X 0 β0 ])0 f0 (X 0 β0 + α0 )]The fundamental theory related to the linear (or more generally differentiable) functional
is given by Van der Vaart (1991). We also refer readers to Chapter 3 in Groeneboom
and Wellner (1992) and Chapter 25 in Van Der Vaart (1998) for more comprehensive
discussions. We draw on Geskus and Groeneboom (1996, 1997, 1999) where the authors
develop a systematic approach to characterize the linear functional of NPMLE for the24interval censored data (case 2). For that purpose, we define c1 (u) = −E[X|u]F0 (u)f0 (u),c3 (u) = −F0 (u+α0 )f0 (u), and c(u) = (c01 (u), c3 (u))0 . Consider the linear functional κ(F0 ) =Rc(v)dF0 (v) and its canonical (with zero mean) gradientZ(4.22)
κ̃F0 (u) = c(u) − c(v)dF0 (v).
The canonical gradient is also known as the “efficient influence function” and is supposed toRbelong to the space L02 (F ); i.e., the space of square integrable functions satisfying adF =0, Groeneboom and Wellner (1992). A key component in determining the asymptotic
property of (α̃n , β̃n ) is κ(F̃n (·; α0 , β0 )); i.e., the linear functional of the NPMLE when thefinite dimensional parameter is set to be its true value. The influence function of the latter
one crucially depends on whether there is a unique element φF satisfying
L∗ φF = κ̃F ,(4.23)

given the differentiability of κ̃F in the sense of Van der Vaart (1991), where L∗ denotes theadjoint operator of L defined in equation (4.24). We further denote its derivative by κ̃0F .
To present the solution φF0 for the true distribution function, we denote x0 = x0 β0 andthe support of it as [CL , CU ]; see Coppejans (2007). For any function a in the tangent set,
the score operator for the nonparametric component is
Rx
R x +αRC
δ1 CL0 adF
δ2 x00 0 adF(1 − δ1 − δ2 ) x0U+α0 adF(4.24) L[a](x0 , δ1 , δ2 ) =+−F (x0 )
F (x0 + α0 ) − F (x0 )1 − F (x0 + α0 )For any function b(u, δ1 , δ2 ), we also have the adjoint operator L∗ specified as follows:(4.25)
Z CU
Z x0
Z x0 −α∗
L [b](x0 ) =b(u, 1, 0)g0 (u)du +b(u, 0, 1)g0 (u)du +b(u, 0, 0)g0 (u)du.
x0

Let φF (u) ≡(4.26)
∗

Ru
CL

x0

CL

a(v)dF (v) denote the integrated score function, then we have

Z CU

L L[a](x0 ) =

x0 −α

Z x0
Z x0 −α0φ(u + α0 ) − φ(u)φ(u + α0 )φ(u)g0 (u)du+g0 (u)du+g0 (u)du.
F (u)1 − F (u + α0 )x0 −α0 F (u + α0 ) − F (u)CLThen one can solve the equation and obtain
L[a](u, δ1 , δ2 ) = − δ1 ((1 − F0 (u))ω(u) + (1 − F0 (u + α0 )ω(u + α0 ))+ δ2 (F0 (u)ω(u) − (1 − F0 (u + α0 )ω(u + α0 ))+ (1 − δ1 − δ2 ) (F0 (u)ω(u) + F0 (u + α0 )ω(u + α0 )) ,

250

where ω(u) ≡ gc0(u). Moreover, the influence function takes the following form (see Example(u)
4.2 of Van de Geer (1997)):(4.27)

−F (u) ((1 − F (u))ω(u) + (1 − F (u + α )ω(u + α )) , for C ≤ u ≤ α
0
0
0
0
0
L
0
φF0 (u) =(1 − F0 (u)) (F0 (u)ω(u) + F0 (u − α0 )ω(u − α0 )) , for α0 ≤ u ≤ CU .
We need one more set of assumptions to guarantee the asymptotic normality for the
linear functional of NPMLE; see [p.31] in Van de Geer (1997).
Condition 7. We assume ω(u) is uniformly bounded for all u in the support. Moreover,the following ratios are all uniformly bounded:ω 0 (u)ω 0 (u)≤ c, sup≤ c,
f0 (u)f0 (u + α0 )uu
ω 0 (u + α0 )
ω 0 (u + α0 )sup≤ c, sup≤ c,
f0 (u)f0 (u + α0 )uu
f0 (u)f0 (u + α0 )≤ c, sup≤ c,
supf0 (u)f0 (u + α0 )uu
supfor some universal finite constant c. Also, for any α in the parameter space, we assume
F0 (u + α) − F0 (u) is uniformly bounded away from zero for any u in the support.
Denote the stacked estimator for the finite dimensional parameter as θ̃n ≡ (α̃n , β̃n0 )0 andthe true unknown parameter as θ0 ≡ (α0 , β00 )0 .
Theorem 4.3 (Consistency of the joint estimator).
Under our Conditions (1)-(7), we have 
for all large n, the zero-crossing θ̃n for Ψn θ̃n exists with a probability tending to one andis a consistent estimator of θ0 .
Theorem 4.4 (Asymptotic normality of the joint estimator). Under our Conditions (1)(7), we have√ 
(4.28)n θ̃n − θ0 ⇒ N(0, Σ̃0 ),where Σ̃0 = H0−1 E[(φ00 + φ0F0 )0 (φ00 + φ0F0 )]H0−1 , φF0 defined in equation (4.27) and!
[∆1i − F0 (Xi0 β0 )]Xi
φ0 =.
∆3i − 1 + F0 (Xi0 β0 + α0 )Remark 4.2. It is known that the NPMLE is indeed more efficient than the isotonic
estimator only using binary choices data in the sense that both n1/3 (F̃n (t; α0 , β0 ) − F0 (t))
and n1/3 (F̂n (t; β0 ) − F0 (t)) converge to the Chernoff distribution yet with different scaling
constant terms. Specifically, the NPMLE has a smaller asymptotic variance than the one

26

associated with the isotonic estimator; see Exercise 4.27 in Chapter 4 of Groeneboom andJongbloed (2014). Our simulation results also confirm this theoretical claim.
Remark 4.3. Apropos of the asymptotic covariance matrices for finite dimensional parameters, the comparison between our two-stage and joint estimation is not obvious analytically,
as both influence functions are very complicated. As the joint approach simultaneously estimates α0 and β0 , one may naturally expect that it works better. This is supported by our
simulation studies. Overall, performances of both estimators are reliable and not subject tothe sensitivity of tuning parameter selection.
Remark 4.4. Referring to the semiparametric efficient score function in Coppejans (2007),
it is inevitable to use either the kernel or sieve approach, because it involves the density
function of latent error. Embedding shape restriction, we suggest an efficient estimation
procedure making use of the maximum smoothed likelihood estimator (MSLE) in Groeneboom (2014), which can also be viewed as the generalization of Klein and Spady (1993) tothe ordered response data. Consider a kernel density function k(·) and kbn (·) ≡ k(·/bn )/bnPwith the bandwidth equal to bn . Define h̄jn (t; β) = n1 ni=1 kbn (Xi0 β, t)∆ji , for j = 1, 2, 3.
The MSLE F̄nbn (t; α, β) is now defined by maximizingZZZ
h̄1n (t; β) log F (t)dt+ h̄3n (t+α; β) log[1−F (t+α)]dt+ h̄2n (t; β) log[F (t+α)−F (t)]dt.
Thereafter, the estimators (ᾱn , β̄n ) for the finite dimensional parameters are defined by
maximizing the profile (smoothed) likelihood function:(4.29)(ᾱn , β̄n ) = arg max Ln (F̄nbn (·; α, β); α, β).
α,βWe conjecture that this (kernel) smoothed maximum likelihood procedure delivers semiparametric efficient estimators for (α0 , β0 ), based on the results in Coppejans (2007) and
Groeneboom (2014). One notable advantage over the sieve method in Coppejans (2007) isthat the estimated distribution is monotone by construction, whereas Coppejans (2007) has
to impose additional restrictions on spline coefficients to accommodate the monotonicity.
Remark 4.5. For some empirical questions, one might be interested in the marginal effectassociated with the ordered response model. In this case, one needs the joint normalityof the estimated finite dimensional parameter and the distribution function evaluated at a
certain point. However, the point-wise asymptotic distribution of the isotonic estimator F̂nor the NPMLE F̃n is not normal. Instead, they both converge to the Chernoff distribution13
with the cubic root rate. In order to restore the asymptotic normality and improve the rate
13A random variable is said to follow the Chernoff distribution if it is the last time where standard two-sided Brownian motion minus the parabola t2 reaches its maximum; see Section 3.9 in Groeneboom andJongbloed (2014).

27

of convergence for the distribution function, one could resort to the smoothed maximum
likelihood estimator (SMLE) or the maximum smoothed likelihood estimator (MSLE) in
Groeneboom, Jongbloed, and Witte (2010) and Groeneboom (2014).

4.3. Bootstrap Inference
In this part, we propose corresponding novel bootstrap procedures for conducting inference for our semiparametric estimators. In a closely related paper for binary choice data,
Groeneboom and Hendrickx (2017) first prove the bootstrap validity for the regression coefficient using Efron’s multinomial weights (Efron (1979)) based on the shape-restricted
estimation in Groeneboom and Hendrickx (2018). We establish the bootstrap consistency
with more general exchangeable bootstrap weights (Mni )ni=1 . Commonly used exchangeablebootstrap schemes include (i) Nonparametric bootstrap in which the bootstrap weight0
s Mn = (Mn1 , ..., Mnn ) follow the multinomial distribution M ult (n, (n−1 , ..., n−1 )); (ii)P
Bayesian bootstrap in which the bootstrap weights Mni = ωi / ( ni=1 ωi ) for 1 ≤ i ≤ nand ωi as the unit exponential distribution (Rubin, 1981); and (iii) Delete-h jackknives
in which the bootstrap weights are generated from permuting the deterministic weights
wni = n/ (n − h) for 1 ≤ i ≤ n − h and let Mnj = wnRn (j) where Rn is a random permutation uniformly over {1, ..., n} (Wu, 1990). All these bootstrap procedures satisfy thefollowing Condition 8. A detailed verification can be found in Cheng and Huang (2010).
0

Condition 8. (i) The vector Mn = (Mn1 , ..., Mnn ) is exchangeable for all n.
P(ii) Mni ≥ 0 for all n, i and ni=1 Mni = n for all n.
(iii) For some positive constant C < ∞, lim supn→∞ kMn1 k2,1 ≤ C, where kMn1 k2,1 =R∞p
Pr (Mn1 ≥ u)du.
0(iv) limλ→∞ lim supn→∞ supt≥λ t2 Pr (Mn1 > t) = 0.
P(v) (1/n) ni=1 (Mni − 1)2 →p 1.
P
In analog with the notations from the empirical process theory, we let P∗n f = n−1 ni=1 Mni f (Zi ),P
and G∗n f = n−1/2 ni=1 (Mni − 1) f (Zi ). Before stating the theoretical results for the bootstrapped quantities, we should be explicit about the underlying probability space and sourceof randomness. There are two sources of randomness coming from the observed data andthe bootstrap weight Mn . The following set of definitions and notations are adapted fromCheng and Huang (2010). We have the product probability space(Z ∞ × M, A∞ × Ω, PZM )for the joint randomness from observed data and bootstrap weights. Furthermore, thebootstrap weights are independent from the sample observations, i.e., PZM = PZ × PM .

28

Definition 4.1. For a real-valued random variable ∆n , we define (i) ∆n = oPM (1) , inPZ −probability if for any ε, δ > 0

PZ PM |Z (|∆n | > ε) > δ → 0.
(ii) We define ∆n = OPM (1) , in PZ −probability if for any δ > 0 there exists a C s.t.

PZ PM |Z (|∆n | > C) > δ → 0.

4.3.1. Bootstrap Two-stage Estimation. We now describe our bootstrap estimator forthe regression coefficients and threshold parameter. Note that bootstrap for regression
coefficients in the binary choice model can be found in Groeneboom and Hendrickx (2017).
Stage 1(i)*. First of all, the bootstrap MLE F̂n∗ (·, β) is computed using the weighted
cumulative sum diagram formed by the point (0, 0) and
!
i
i
X
X
(β)
(β)
(β)
Mn(j) ,
Mn(j) ∆1,(j) ,
j=1
(β)

j=1
(β)

where Mn(i) corresponds to the weight attached to U(i) .
Stage 1(ii)*. The bootstrap estimator of the regression coefficient β̂n∗ is defined as thezero-crossing point of the following estimating equations:nhi
1X(4.30)
Mni Xi ∆1i − F̂n∗ (Xi0 β̂n∗ ; β̂n∗ ) = 0.
n i=1Stage 2*. Finally, the bootstrap version α̂n∗ is then defined as the zero-crossing point of
the following estimating equation:nhi1X
Mni 1 − ∆3i − F̂n∗ (Xi0 β̂n∗ + α̂n∗ ; β̂n∗ ) = 0.
(4.31)
n i=1
Given Condition 8, the bootstrap estimate F̂n∗ is a step-wise monotone function. Therefore, the estimating equation for α̂n∗ is again monotone so that the computational advantage
of our approach is amplified along the bootstrap replications. Regarding the theoretical underpinning, one could easily prove that α̂n∗ → α0 conditional on observations (Z1 , · · · , Zn )almost surely. We characterize the conditional weak limit for α̂n∗ in the next theorem. For
completeness, we also state the result for β̂n∗ , which has been established by Groeneboom
and Hendrickx (2017)[p3464-3465].

Theorem 4.5 (Bootstrap validity for the two-stage estimator). Suppose Conditions (1)(6) and Condition (8) hold. For the bootstrap estimators β̂n∗ and α̂n∗ with exchange weights

29(Mn1 , ..., Mnn ), we have thatβ̂n∗ − β̂n = OPM n−1/2 , and |α̂n∗ − α̂n | = OPM n−1/2 ,in PZ −probability. Furthermore,√  ∗(4.32)n β̂n − β̂n ⇒ N(0, Ωβ0 ),√(4.33)× N(0, Ωα0 ),n (α̂n∗ − α̂n ) ⇒ Vα−10conditional on observations (Z1 , · · · , Zn ) almost surely.

A direct consequence of the above theorem is the validity of percentile bootstrap confidence interval (see Cheng and Huang (2010)). The lower p-th quantile of the bootstrap∗∗
distribution of α̂n∗ is the quantity τn,α(p) satisfying τn,α(p) ≡ inf{τ : PM |Z (α̂n∗ ≤ τ ) ≥ p}.
∗∗
Similar, we can define τn,β(p) ≡ inf{τ : PM |Z (β̂n,k≤ τ ) ≥ p} for the kth component of β̂n∗ ,kk = 2, ..., K. Then the percentile-type bootstrap confidence intervals can be constructed as ∗
∗d α (p) = τn,α(4.34)
BC(p/2), τn,α(1 − p/2) ,andd β,k (p) = τ ∗ (p/2), τ ∗ (1 − p/2) , for k = 2, · · · , K.
BC
n,βk
n,βk

(4.35)

Corollary 4.3. Suppose Conditions (1)-(8) hold, then we haved(4.36)PZM α0 ∈ BC α (p) → 1 − p,andd β,k (p) → 1 − p,
PZM β0k ∈ BC

(4.37)

for k = 2, · · · , K,

as n → ∞.

4.3.2. Bootstrap Joint Estimation. Now we describe the bootstrap inference for our
joint estimation approach.
Joint Estimation*. For any α and β, we can derive the bootstrap NPMLE for F (·; α, β)based on the ordered response data:F̃n∗ (·; α, β)= arg max
F

n
X
i=1

Mni [∆1i log F (Xi0 β) + ∆2i log (F (Xi0 β + α) − F (Xi0 β)) + ∆3i log (1 − F (Xi0 β + α))] ,30

which exists. Given F̃n∗ (·; α, β) from the previous step, we define (α̃n∗ , β̃n∗ ) which are thezero-crossing points of the estimating equations simultaneously:nih
1X(4.38)
Mni Xi ∆1i − F̃n∗ (Xi0 β̃n ; α̃n , β̃n ) = 0,
n i=1
n
h
i
1X
Mni 1 − ∆3i − F̃n∗ (Xi0 β̃n + α̃n ; α̃n , β̃n ) = 0.
n i=1Theorem 4.6 (Bootstrap validity for the joint estimator). Suppose Conditions (1)-8) hold.
For the bootstrap estimator θ̃n∗ with exchange weights (Mn1 , ..., Mnn ), we have thatθ̃n∗ − θ̃n = OPM n−1/2 ,in PZ −probability. Moreover, we have√  ∗n θ̃n − θ̃n ⇒ N(0, Σ̃0 ),(4.39)
conditional on observations (Z1 , · · · , Zn ), almost surely.
Similarly as the two-stage estimator (part 4.3.1), we can construct the percentile-typeg α and BC
g β,k of the regression coefficients and thresholdbootstrap confidence intervals BC
parameter for the joint estimator. The coverage properties are summarized in the following
corollary.
Corollary 4.4. Suppose Conditions (1)-(8) hold, then we haveg(4.40)PZM α0 ∈ BC α (p) → 1 − p,and(4.41)g β,k (p) → 1 − p,
PZM β0k ∈ BC

for k = 2, · · · , K,

as n → ∞.

5. Numerical Results
In this section, we conduct Monte Carlo simulations to evaluate the finite sample performances of our proposals and other semiparametric approaches described in Section 3.3.
The simulation designs follow the interdependent duration model of Honoré and de Paula(2010). We also illustrate the usefulness of all these methods by applying them to theDanish joint retirement data in An, Christensen, and Gupta (2004).

31

5.1. Monte Carlo Studies
Recall the model setup in Section 2 of the interdependent durations model proposedby Honoré and de Paula (2010): two players respectively decide (T1 , T2 ) as the time between switching from an initial activity to an alternative activity. The utility flows fromthe initial activity for players 1 and 2 are random variables (1 , 2 ) with a joint CDF
G(·, ·). The utility flow from the alternative activity at time t gives player j an amountta exp(Xj0 β0 ) exp[α∗ I {Tk ≤ t}] for j, k ∈ {1, 2} and j 6= k. Here we consider individuals
with covariates denoted by Xj = (1, X1j , X2j , X3j , X4 ), where the last coordinate is assumed to be the same for both players, and write the corresponding regression coefficient
as β0 = (β0o , β01 , β02 , β03 , β04 )0 .
The main complication of the model is the characterization of (T1 , T2 ) as the Nash equilibrium outcomes. Based on the value of 1 /2 , the equilibrium (T1 , T2 ) can be categorized
into three cases: T1 < T2 , T1 = T2 , and T1 > T2 . Specifically, Honoré and de Paula(2010) characterize the equilibrium (T1 , T2 ) as follows: if 1 /2 < exp[(X1 − X2 )0 β0 − α∗ ],then T1 = [1 / exp(X10 β0 )]1/a and T2 = [2 / exp(X20 β0 + α∗ )]1/a ; if exp[(X1 − X2 )0 β0 − α∗ ] <1 /2 < exp[(X1 − X2 )0 β0 + α∗ ], then T1 = T2 and can take any value betweenhhi1/a hi1/a 
i1/a hi1/a 
(X10 β0 +α∗ )(X20 β0 +α∗ )(X10 β0 )(X20 β0 )max 1 /e, 2 /e, min 1 /e, 2 /e;if 1 /2 > exp[(X1 −X2 )0 β0 +α∗ ], then T1 = [1 / exp(X10 β0 + α∗ )]1/a and T2 = [2 / exp(X20 β0 )]1/a .
When exp[(X1 −X2 )0 β0 −α∗ ] < 1 /2 < exp[(X1 −X2 )0 β0 +α∗ ], the model produces multiple
equilibria so that in case of the simultaneous failure (T1 = T2 ) there exists a range of possible values. However, the equilibrium selection rule does not matter here when it comes toestimating the regression coefficient β0 or the interaction effect α∗ , because the conditional
probability P {T1 = T2 |X1 , X2 } remains the same. Therefore, one can transform the above
interdependent durations model into the ordered responses model with Y = 1 if T1 < T2 ,Y = 2 if T1 = T2 and Y = 3 if T1 > T2 . Then the conditional probabilities become:(5.1) P {Y ≤ 1|X1 , X2 } = H (β01 (X11 − X12 ) + β02 (X21 − X22 ) + β03 (X31 − X32 ) − α∗ ) ,(5.2) P {Y ≤ 2|X1 , X2 } = H (β01 (X11 − X12 ) + β02 (X21 − X22 ) + β03 (X31 − X32 ) + α∗ ) ,where H(w) = P {log(1 /2 ) ≤ w}. Furthermore, H(·) is an arbitrary distribution function
if no assumption is made regarding (1 , 2 ). Note that the common effect from covariates’ effect β0o +β04 X4 is differenced out here. Also, we normalize β01 = 1 for the identification purpose. In our simulations, we set (a, βo , β1 , β02 , β03 , β04 , α∗ ) = (1.35, −4.00, 1.00, 1.00, 1.00, 0.50, 1.00).
We consider two cases for G(·, ·): the joint distribution of two independent unit exponential variables or two independent log normal variables. Variables (X11 , X12 , X31 , X32 , X4 )

32

are independent standard normal random variables and (X21 , X22 ) are independent χ2 (1)random variables (normalized to a mean zero and variance one). Our Monte-Carlo study
involves 1,000 simulations of the foregoing models. The sample sizes are n = 250, 500, and750.
Table 1 presents the finite sample bias, root of the mean square error (RMSE), and median absolute deviation (MAE) of our Isotonic two-stage estimator and NPMLE-based joint
estimator for (β02 , β03 , α∗ ). The error terms (1 , 2 ) are two independent unit exponential or
log-normal random variables. The initial values of (β2 , β3 , α∗ )0 are set as (0.2, 0.2, 0.1)0 andthe estimating equations are solved using R package BB. As a comparison, Table 2 presents
the finite sample results for three alternative semiparametric estimators described in Remarks 3.1 to 3.3; i.e., Klein and Sherman (2002)’s kernel-based approach (K-S), Melenberg
and Van Soest (1996)’s smoothed maximum score estimator (SMS), and the rank estimator
combining Cavanagh and Sherman (1998) and Chen (2002). For the K-S approach, thebandwidths (including a bandwidth, a pilot bandwidth, and a smoothing parameter in thedamping function) are chosen following the guidelines of Klein and Sherman (2002).14 The
trimming parameter ξˆ for the quasi-likelihood function is set to be the 0.95th sample quantile of the Euclidean norm of covariates. We consider two values for trimming proportion
for the target set: p = 0.05 and p = 0.20. For the SMS approach, we experiment two
bandwidths n−1/4 and n−1/5 . The former reduces bias by under-smoothing and the later
is chosen according to the MSE-optimal rate. Table 3 examines point-wise estimates forthe error distribution function H(·) evaluated at 5 equally spaced points in the support of(log 1 − log 2 ).
We first evaluate the performances of our two-stage estimator and joint estimator. Table
1 shows that bias, RMSE, and MAE for both estimators decrease when the sample size
increases. NPMLE-based joint estimation yields smaller bias than the two-stage estimator.
The RMSE’s are close between two estimators, with the two-stage estimator slightly betterfor (β02 , β03 )0 and the joint estimator slightly better for α∗ . The lower variability of twostage estimator is likely due to its computational simplicity. As shown in Figures 1 and
2, the estimating functions exhibit smoother variation in the two-stage procedure than in
the joint procedure. Nevertheless, joint estimator still yields smaller RMSE for α∗ due
to the fact that it utilizes all the information across three categories by the simultaneous
estimation. Also, the MAE of the joint estimator is slightly smaller than that of twostage estimator. In terms of the nonparametric estimate of H(·), Table 3 shows that the
14The value δ in Klein and Sherman (2002) p669 is set to 1/6. The rate of bandwidth α is set as the middle

point of the allowed range ((3 + δ)/20, 1/6), see p670. The rate of the pilot bandwidth follows Lemma 5Aof Klein and Sherman (2002). Finally,  in the damping function is set as the middle point of the allowed
range (0, 1/40 − δ/20), see p670.

33

NPMLE-based approach performs better in general than the isotonic estimator, especiallyfor H(1) and H(2). This confirms the efficiency gain from utilizing additional informationthat differentiates the category with Y = 2 and the one with Y = 3 (see Remark 4.2).
We further compare our two-stage and joint estimator with other semiparametric estimators as reported in Table 2. The K-S estimator for α∗ is sensitive to the trimming
parameter p, which decides the subset used for averaging individual estimators (see Remark 3.1). It is clear that p = 0.05 does not sufficiently exclude individual estimators that
perform poorly and thus leads to unreliable estimate for α∗ . The performance for α∗ greatly
improves when p = 0.20; however, the RMSE and MAE are still larger than our two-stage
and joint estimators. The SMS estimator has significantly larger RMSE and MAE than allother methods. This can not be fully explained by the under-smoothing bandwidth n−1/4 .
When bandwidth takes the MSE-optimal rate n−1/5 , the RMSE and MAE of SMS decrease,yet still larger than other methods. The rank estimator described in Remark 3.3 exhibits
remarkable finite sample performances in our simulations. It yields the smallest bias forregression coefficients and threshold parameter among all methods. However, referring to
the RMSE for α∗ , our two-stage and joint estimators still have an edge on the rank estimator. In sum, our simulation studies show that the isotonic two-stage estimator and the
NPMLE-based joint estimator perform well and stable. Their tuning-parameter-free advantage is further enhanced by the robust finite sample performances. Namely, in terms of
the RMSE, they clearly outperform the SMS in all scenarios we considered and are better
than the K-S for most cases. Moreover, our estimators, especially the joint estimator, yield
the smallest RMSE for α∗ , which is the parameter of primary interest in the interdependent
durations model.

34

Table 1. Finite sample performances of the isotonic two-stage estimatorand the NPMLE-based joint estimator for (β02 , β03 , α∗ ).
(1 , 2 )Exponential
Log-normal
MethodsnBias RMSE MAE
Bias RMSE MAE
Two-stage 250
β02-.0658 .1985 .1448-.0482 .1701 .1202
β03-.0692 .1717 .1234-.0547 .1582 .1199∗
α-.0663 .1470 .1092-.0496 .1364 .0982
500
β02-.0470 .1370 .0937-.0345 .1255 .0841
β03-.0485 .1259 .0883-.0333 .1113 .0798∗
α-.0476 .1132 .0810-.0342 .0978 .0715
750
β02-.0425 .1163 .0822-.0297 .0980 .0679
β03-.0423 .1033 .0719-.0294 .0861 .0607∗
α-.0376 .0907 .0643-.0299 .0791 .0556
Joint
250
β02
-.0337 .2079 .1412-.0211 .1786 .1197
β03-.0394 .1774 .1258-.0285 .1647 .1166∗
α-.0516 .1463 .1031-.0378 .1357 .0968
500
β02-.0282 .1372 .0916-.0184 .1272 .0835
β03-.0316 .1276 .0889-.0184 .1120 .0789∗
α-.0389 .1099 .0758-.0281 .0953 .0689
750
β02-.0280 .1190 .0828-0160 .0998 .0694
β03-.0300 .1048 .0717-.0162 .0894 .0623∗
α-.0322 .0889 .0618-.0246 .0776 .0542

35

Table 2. Finite sample performances of other semiparametric estimators,
described in Remarks 3.1 to 3.3. Trimming parameter p in the K-S approachis used to exclude individual estimators for α∗ that performs poorly; h is the
bandwidth used in SMS approach.
(1 , 2 )Exponential
Log-normal
MethodsnBias RMSE MAE
Bias RMSE MAE
K-S
250
β02
.0374 .2430 .1593.0267 .1956 .1256
β03.0373 .1943 .1251.0201 .1585 .1000
p = .05α∗.1568 .2797 .1549.1647 .2681 .1541∗
p = .20α.0697 .2104 .1214.0570 .1791 .1058500
β02.0205 .1435 .0942.0142 .1258 .0849
β03.0164 .1267 .0794.0121 .1075 .0671∗
p = .05α.1682 .2501 .1477.1642 .2232 .1484∗
p = .20α.0446 .1405 .0824.0378 .1160 .0738
750
β02.0194 .1202 .0825.0166 .0958 .0655
β03.0090 .1003 .0622.0080 .0847 .0563∗
p = .05α.1718 .2422 .1452.1688 .2143 .1507∗
p = .20α.0347 .1130 .0680.0356 .0937 .0593
SMS
250
β02
.0776 .3093 .1658.0574 .2683 .1419−1/4(h = n)β03.0728 .2791 .1645.0559 .2326 .1359
α∗.0638 .3046 .1939.0427 .2648 .1658500
β02.0534 .2100 .1310.0379 .1960 .1193
β03.0381 .1882 .1138.0318 .1693 .1094
α∗.0365 .2152 .1468.0250 .1847 .1201750
β02.0428 .1896 .1111.0255 .1564 .1057
β03.0267 .1574 .1001.0236 .1410 .0869∗
α.0215 .1823 .1222.0156 .1638 .1049−1/5(h = n) 750
β02.0468 .1823 .1086.0273 .1419 .0910
β03.0303 .1449 .0877.0243 .1277 .0764∗
α.0273 .1628 .1065.0174 .1470 .0905
Rank
250
β02
.0108 .2082 .1296.0111 .1685 .1108
β03.0112 .1725 .1142.0118 .1478 .0977
α∗.0087 .1639 .1103.0087 .1462 .0905
500
β02.0160 .1311 .0900.0090 .1090 .0711
β03.0137 .1179 .0738.0094 .0988 .0648
α∗.0114 .1171 .0750.0054 .0974 .0651
750
β02.0071 .1041 .0696.0073 .0869 .0581
β03.0042 .0909 .0609.0056 .0790 .0522
α∗.0053 .0921 .0580.0037 .0797 .0513

36

Table 3. Finite sample performances of point-wise estimators for the function H(w), Isotoic estimator uses binary data, NPMLE uses entire ordered
response data.
(1 , 2 )Exponential
Log-normal
Methods n
Bias RMSE MAE
Bias RMSE MAE
Isotonic 250 H(−2) -.0186 .0578 .0417-.0223 .0528 .0430
H(−1) -.0053 .0829 .0580-.0088 .0842 .0602
H(0) .0089 .1043 .0726.0131 .1014 .0720
H(1) .0294 .1018 .0698.0296 .1078 .0728
H(2) .0416 .0834 .0673.0340 .0706 .0785
750 H(−2) -.0091 .0386 .0274-.0089 .0354 .0260
H(−1) .0003 .0558 .0391-.0041 .0551 .0359
H(0) .0083 .0682 .0483.0103 .0714 .0474
H(1) .0154 .0693 .0476.0137 .0679 .0462
H(2) .0172 .0550 .0380.0189 .0482 .0368
NPMLE 250 H(−2) -.0084 .0498 .0342-.0073 .0466 .0328
H(−1) .0053 .0660 .0445.0052 .0734 .0495
H(0) .0183 .0847 .0571.0168 .0873 .0611
H(1) .0298 .0823 .0562.0239 .0872 .0627
H(2) .0216 .0588 .0435.0190 .0540 .0421
750 H(−2) -.0051 .0335 .0243-.0032 .0312 .0221
H(−1) .0044 .0469 .0332.0012 .0449 .0302
H(0) .0094 .0569 .0387.0412 .0618 .0428
H(1) .0148 .0542 .0367.0110 .0535 .0372
H(2) .0117 .0393 .0277.0096 .0337 .0237
K-S
250 H(−2) .0029 .0412 .0279-.0021 . 0346 .0234
H(−1) .0040 .0574 .0404.0019 .0561 .0376
H(0) .0004 .0690 .0479.0075 .0744 .0518
H(1) -.0118 .0688 .0463-.0131 .0703 .0431
H(2) -.0352 .0683 .0426-.0425 .0662 .0399
750 H(−2) -.0003 .0240 .0161-.0014 .0211 .0133
H(−1) .0032 .0337 .0231.0010 .0343 .0230
H(0) .0041 .0412 .0280.0081 .0453 .0298
H(1) -.0029 .0440 .0301-.0042 .0401 .0265
H(2) -.0223 .0411 .0268-.0300 .0407 .0285

37

Table 4 reports the empirical coverage proportion (CP) and the median length (ML) of
the confidence intervals (CIs) based on the nonparametric bootstrap described in Section
4.3, for the our two-stage and the joint estimators. We observe that the coverage proportions of their nonparametric bootstrap based CIs are close to the nominal 95% level. For
the threshold parameter α∗ , the CIs of joint estimator has better coverage proportions andshorter median lengths than the two-stage version. Table 4 also summarizes the nonparametric bootstrap based CIs for other semiparemetric estimators: the K-S, the SMS, and therank estimator. The under-smoothing bandwidth n−1/4 is used for the SMS to reduce thebias in the asymptotic distribution (see Section 4.3.3 of Horowitz (2009)). The trimming
proportion p is set to 0.2 in the K-S estimator. We observe that the SMS exhibits overcoverage in all scenarios and the median length of CIs are substantially longer than othermethods. The K-S leads to under-coverage for α∗ when the n = 750. The rank estimator
performs well when n = 750; however, it suffers from under-coverage for the regression
coefficients (β02 , β03 )0 when n = 250 and 500.

Table 4. Coverage proportions of 95% nonparametric bootstrap based confidence intervals, sample size = n, CP = coverage proportion, ML = medianlength of the confidence interval, number of bootstrap replications = 200.
(1 , 2 ) Exponential(1 , 2 ) log-normaln = 250n = 500n = 750n = 250n = 500n = 750
Methods
CP ML CP ML CP ML
CP ML CP ML CP ML
Two-stage β02 .974 .779 .955 .506 .945 .415
.978 .683 .958 .457 .941 .368
β03 .961 .677 .940 .454 .943 .373
.963 .597 .947 .409 .952 .332∗
α .961 .577 .929 .386 .926 .317
.956 .522 .945 .363 .930 .290
Joint
β02 .931 .737 .945 .520 .948 .431
.952 .681 .929 .470 .944 .386
β03 .945 .664 .942 .470 .934 .386
.949 .609 .948 .424 .951 .347∗
α .950 .514 .920 .372 .934 .307
.952 .484 .930 .344 .944 .279
K-S
β02 .988 .864 .985 .590 .987 .476
.986 .713 .970 .484 .983 .402
β03 .983 .804 .980 .548 .983 .448
.994 .648 .979 .459 .979 .380∗
α .952 .748 .924 .463 .910 .358
.942 .613 .905 .384 .884 .295
SMS
β02 .965 1.317 .977 .849 .976 .709
.980 1.102 .976 .736 .975 .625
β03 .973 1.182 .975 .755 .980 .632
.985 .990 .975 .666 .975 .566
α∗ .978 1.281 .971 .869 .980 .716
.984 1.067 .975 .744 .970 .648
Rank
β02 .831 .512 .907 .431 .936 .373
.843 .452 .915 .374 .946 .323
β03 .869 .499 .935 .410 .939 .343
.884 .434 .929 .348 .947 .297
α∗ .939 .565 .939 .423 .950 .355
.935 .512 .952 .387 .956 .316

38

On the basis of our simulation results, we recommend the NPMLE-based estimatorand the isotonic two-stage estimator in combination with the nonparametric bootstrap
inference. The main convenience of the recommended estimators is that applied researchers
do not need to choose any tuning parameter. Moreover, our simulations show that NPMLEbased joint estimator is more efficient (smaller RMSE and shorter CI median length) thanthe K-S and SMS approaches. The rank estimator performs quite well when the sample size
is relatively large, but is outperformed by our estimators when the sample size is relatively
small.

5.2. An Empirical Application Using Joint Retirement Data
We utilize the interdependent durations model to study the joint retirement decisions ofDanish couples (An, Christensen, and Gupta (2004)) so that T1 and T2 represent the timing
choices of retirement dates for the wives and husbands, respectively. The unique feature inthis type of data is the couples’ synchronization of retirement dates (so T1 = T2 ), despite thedifference in their observable characteristics such as age, education level, or working skill.
This posts new challenges to the econometric analysis, because the probability of simultaneous retirement is assumed to be zero in the traditional (mixed) proportional hazards
models in duration analysis. In sharp contrast, the Honoré-de Paula model accommodates
the presence of simultaneous retirement, and this is generated by strategic interaction ofboth players in the noncooperative game. In the literature, noncooperative games havebeen widely used to model joint retirement behavior of couples in Gustman and Steinmeier(2000) and Gustman and Steinmeier (2004), among others. Also see Lundberg and Pollak(1994) for a general discussion of noncooperative games for within-marriage decisions.
The data of elderly couples are drawn from a random 0.5% sample of the Danish population and are provided by An, Christensen, and Gupta (2004). We restrict the sample tobe couples satisfying the following conditions: (1) Both husband and wife were working in
the base year (defined as the year in which the oldest spouse is 54 in November); (2) At
least one spouse retired by the end of the survey, which allows us to determine the timeorder of retirement for each couple; (3) Couples remained married until either the death of
the spouse or the end of the survey; (4) The base year was 1990 or earlier. The resulting
sample consists of 146 couples. For 29% of the couples in our sample, wives retired at least
one year earlier than husbands; 35% of couples retired in the same year; for the remaining
36% couples, wives retiring later than husbands. As explained in Section 2, the interdependent durations model can be estimated via an ordered response model. We consider twocovariates: age and a skill dummy variable indicating the job category. Table 5 presents
the summary statistics of the covariates for husbands and wives. Note that in the order39

response model, we use the covariate of the wife minus the covariate of the husband. The
coefficient on age difference is normalized to one.
Table 5. Descriptive statistics.
Variables
Mean Std.
Age-wife
50.10 3.82
Age-husband
53.41 2.05
Skill-wife
0.60 0.49
Skill-husband
0.68 0.47
Table 6 presents the estimates of the coefficient β02 on the job category difference and theinteraction parameter α∗ . The minus sign of the point estimate of β02 suggests that if a wife
was categorized as a skilled employee while her husband was not, then the probability that
she retired earlier was smaller, which can be driven by the family-wise income effect. Our
two estimators and the rank estimator yield similar point estimates for β02 . Nevertheless,such an effect of job skill difference on the retirement timing turns out to be insignificant at
the 5% significance level except for the SMS estimator with an under-smoothed bandwidth(n−1/4 ) and the rank estimator. On the other hand, the interaction parameter α∗ is significantly positive in all approaches except for K-S estimator, suggesting that a spouse gained
extra utility from life after retirement together with her/his counterpart. In particular, theIsotonic two-stage estimator indicates that couples to synchronize their retirement timingwith a force large enough to compensate on average 4.76 years of age difference in the SSE
estimate. That number rises to 6.17 years for the NPMLE-based joint estimator. Other
approaches yield comparable point estimates. In terms of the 95% confidence intervals, ourtwo-stage and joint estimators yield results close to the rank estimator. The vast majorityof values in their confidence intervals are large enough to counter the average age difference
3.31 of the couples15. This evidence is also consistent with previous empirical results basedon various structural models (Gustman and Steinmeier (2000, 2004)) where the essence
is to argue the interaction effect is significant enough to counter the age difference whencouples retire simultaneously. It is notable that SMS estimator produces substantially larger upper bound for α∗ and smaller lower bound for β02 . This is not surprising given theover-coverage reported in our Monte-Carlo results. As a matter of fact, in a different empirical application, Bellemare, Melenberg, and Van Soest (2002) also found that the SMS
estimates are far away from other methods including the ordered Probit, the partial linear,and the semi-parametric least square, among others16.
15Recall that the coefficient on the age difference is normalized to one.
16See p194-195 of Bellemare, Melenberg, and Van Soest (2002) for details.

40

Table 6. Estimates of an interdependent duration model applied to joint
retirement in Denmark, n = 146, β01 (coef. on age difference) is normalized
to 1. The 95% confidence intervals are constructed from 200 bootstrap samples.
Two-stage

Coef.
95% CI
Joint
Coef.
95% CI
K-S
Coef.
(p = 0.05) 95% CI
K-S
Coef.
(p = 0.20) 95% CI
SMS
Coef.
(h = n−1/5 ) 95% CI
SMS
Coef.
−1/4
(h = n) 95% CI
Rank
Coef.
95% CI

β02 (job category diff.)-2.35[-12.69, 2.32]-3.21[-9.54, 3.32]-0.14[-12.24, 4.02]-0.14[-12.24, 4.02]-5.20[-23.90, 0.01]-5.22[-24.14, -0.32]-2.99[-10.56, -0.40]

α∗ (interaction)
4.76[3.06, 14.31]
6.17[2.90, 15.96]
5.76[-1.71, 9.55]
7.60[-0.57, 9.34]
4.70[4.16, 23.60]
4.83[3.99, 23.94]
8.20[3.15, 14.40]

Figure 3. The estimated CDF of log(1 /2 ): the distribution of the log
ratio of wife-husband unobservables.

41

Figure 3 plots the estimated CDF of the error term in the ordered response model, whichcorresponds to the distribution of the log ratio of unobservables for wives and husbands,
H(·). By construction, both our isotonic and NPMLE estimates for the CDF are nondecreasing step functions. In contrast, the K-S approach yields an estimated CDF thatfluctuates and has some local maximum and minimum points.

6. Conclusion and Extension
In this paper, we have proposed two simple semiparametric estimation methods for ordered response models with an unknown latent error distribution. We rigorously establish
the asymptotic properties of finite dimensional parameters, tackling the challenging issues related to the nonparametric components based on NPMLE. Our methods are easy
to implement and free of any tuning parameter. Also, the methodology is directly applicable to estimate the social interaction effect in the interdependent durations model by
Honoré and de Paula (2010). Complementing several important contributions as in Lee(1992), Klein and Sherman (2002), Lewbel (2002), and Coppejans (2007), our work makes
the Honoré-de Paula model a viable benchmark in analyzing multiple durations data withstrategically-interacting agents. Both the Monte Carlo simulation and a real data application demonstrate the utility of our approach.
In conclusion, we briefly draw the reader’s attention to several interesting questions thatcan provide future research problems in this area. First of all, a rigorous investigationof an efficient procedure proposed in Section 4.2 need to be delegated to another paper,
considering the amount of details required. A second extension involves modeling the joint
retirement phenomenon by cooperative games as considered by Honoré and de Paula (2018).
The Nash bargaining solution in Honoré and de Paula (2018) can not be written as thestandard ordered response model, unlike the non-cooperative game model in Honoré andde Paula (2010). When it comes to the estimation and inference of the model, Honoré andde Paula (2018) resort to the indirect inference, which is parametric and computationallyintensive. It remains a challenge to develop a flexible semiparametric estimation procedure for the model in Honoré and de Paula (2018). Third, one could apply our estimation
methods in estimating a sample selection model where the selection equation involves making ordered choices; see Section VI.A on [p.147] of Vella (1998). An interesting example
is presented by Attanasio, Koujianou, and Kyriazidou (2008) in which they study credit
constraints faced by consumers in the market for durable goods. Last but note least, it isworthwhile to explore the shape-restricted estimation of the ordered response model withheteroskedasticity (Lewbel, 2000; Chen and Khan, 2003). A thorough treatment of all these
extensions is outside the scope of this paper and will be pursued in the future.

42

7. Appendix A: Proofs of Main Results
This appendix provides proofs for the theorems. The lemmas used for proofs (denoted
as Lemma S1’s and S2’s) are collected in the supplemental note.
We denote some positive constants by c or C whose value might change line by line.

Proof of Theorem 4.1. Part (i) has been shown in Theorem 4.1 [p1426] of Groeneboom andHendrickx (2018). Here we proof part (ii). Note that Remark 3.1 [p1423] of Groeneboom
and Hendrickx (2018) showed that F̂n (u; β̂n ) converges uniformly to F0 (u). We start withthe following estimating equation in Stage 2:
n
i
1 Xh
0
(7.1)
Ψn (α) =1 − ∆3i − F̂n (Xi β̂n + α; β̂n ) ,n i=1
and its probability limit(7.2)

Ψ(α) = E [1 − ∆3 − F0 (X 0 β0 + α)] .

We first show the following uniform convergence result(7.3)

sup |Ψn (α) − Ψ(α)| →p 0.
αObserve that
(7.4)

h
i
sup |Ψn (α) − Ψ(α)| ≤ sup |(Pn − P ) ∆3 + F̂n (X 0 β̂n + α; β̂n ) |
α
α
i
h
0
0
+ sup |P F̂n (X β̂n + α; β̂n ) − F0 (X β0 + α) |.
αHence, the uniform convergence (7.3) hold by Glivenko-Cantelli property in Lemma (S1.4)
and the fact that F̂n (u; β̂n ) converges to F0 (u) in sup-norm.
Next we show the existence of a unique zero-crossing point with probability approaching
to 1. Because α is a scalar, the zero-crossing point of Ψn (α) can be equivalently defined asα̂n such that for any α:(7.5)(α̂n − α)Ψn (α) ≥ 0,

see Lemma 4.1 of Groeneboom and Hendrickx (2018). If the zero-crossing point does notexist, then for all α1 there exists some α2 such that
(7.6)

(α1 − α2 )Ψn (α2 ) ≤ −c < 0,for some finite positive constant c. Such a constant term c exists because the isotonic
estimate F̂n (u; β̂n ) is a piece-wise constant function with finitely many jumps for any n, so

43

is Ψn (α) for all α. In particular, we have
(α0 − α2 )Ψn (α2 ) ≤ −c,(7.7)
By (7.3), we get(α0 − α2 )Ψ(α2 ) ≤ −c/2,(7.8)

with a probability tending to 1. However, this contradicts the fact that α0 is the uniquezero crossing point of Ψ(α), since Ψ(α) is monotone and continuous w.r.t. α, given themonotonicity and absolute continuity of F0 . Thus, the zero-crossing point α̂n exists with aprobability tending to 1.
Finally considering the monotone estimating equation Ψn (α) given the monotonic NPMLE F̂n , the consistency result is a direct consequence of Lemma 5.10 in Van Der Vaart(1998).

Proof of Theorem 4.2. The linear representation of β̂n and its asymptotic normality havebeen shown in Theorem 4.1 of Groeneboom and Hendrickx (2018). Here we focus on α̂n .
Following Groeneboom and Hendrickx (2018), we can define the value of Ψn at α̂n by
setting Ψn (α̂n ) = 0. Note that Ψn (α̂n ) is the convex combination of the left and right limit
at α̂n :(7.9)

Ψn (α̂n ) = λΨn (α̂n −) + (1 − λ)Ψn (α̂n +) = 0,

where we can choose λ ∈ [0, 1] such that (7.9) holds. To prove the root-n consistency andasymptotic normality of α̂n , we start with the estimating equation Ψn (α̂n ) = 0 and we
decompose the l.h.s. asni
1 Xh
1 − F̂n (Xi0 β̂n + α̂n ; β̂n ) − ∆3i = I1n + I2n + I3n ,(7.10)
n i=1
where
n

1X
[1 − F0 (Xi0 β0 + α0 ) − ∆3i ] ,n i=1(7.11)

I1n =(7.12)

n
i
1 Xh
I2n =
F0 (Xi0 β0 + α0 ) − F̂n (Xi0 β0 + α0 ; β0 ) ,n i=1(7.13)

n
i
1 Xh
I3n =
F̂n (Xi0 β0 + α0 ; β0 ) − F̂n (Xi0 β̂n + α̂n ; β̂n ) .
n i=1

Here F̂n (Xi0 β0 + α0 ; β0 ) is the (infeasible) NPMLE computed using the true unknown β0 .
Apparently, the term I1n is of Op (n−1/2 ) with its influence function equal to ψ0 as definedin our Theorem 4.2.

44ab
Referring to I2n , we get I2n = I2n
+ I2n
where
h
i
a
(7.14)I2n= P F0 (U + α0 ) − F̂n (U + α0 ; β0 )
and
h
i
b
I2n = (Pn − P ) F0 (U + α0 ) − F̂n (U + α0 ; β0 ) .
b
We shall utilize P -Donsker property (Van Der Vaart and Wellner (1996)) to show I2n
=
−1/2
a
op (n
) in Lemma (S1.4), and we obtain the linear representation for I2n as follows√ a(7.15)
nI2n = Gn ψF0 + op (1),in our Lemma (S1.9). When it comes to I3n , we decompose it into three terms I3n =
c
b
a
where
+ I3n
+ I3n
I3n
(7.16)hiaI3n= P F0 (X 0 β̂n + α̂n ; β̂n ) − F0 (X 0 β0 + α0 ) ,
h
i
b
0
0
0
0
I3n = P F̂n (X β̂n + α̂n ; β̂n ) − F̂n (X β0 + α0 ; β0 ) − F0 (X β̂n + α̂n ; β̂n ) + F0 (X β0 + α0 ) ,
h
i
c
I3n
= (Pn − P ) F̂n (X 0 β̂n + α̂n ; β̂n ) − F̂n (X 0 β0 + α0 ; β0 ) .
cb== op (n−1/2 ), I3n
In Lemma S1.11 and Lemma S1.4 respectively, we prove that I3n
−1/2
op (n
) using P -Donsker Property of related functional classes. Also we have the following
expansion:(7.17)aI3n= Vα0 (α̂n − α0 ) + Vβ0 (β̂n − β0 ) + op (n−1/2 + α̂n − α0 + |β̂n − β0 |).

In sum, the desired linear representation for α̂n follows after collecting the leading termsaaand substituting the linear representation for β̂n based on (4.5).
and I3nin I1n , I2n

Proof of Theorem 4.3. Given the compactness of the parameter space, the sequence of estimator θ̃n has a subsequence θ̃nk converging to some element θ∗ = (α∗ , β ∗0 )0 . In Appendix
C of the supplemental note, we apply Theorem 7.4 in Van de Geer (2000) to show thefollowing convergence in terms of the Hellinger distance (Van de Geer (1993)):
sup h(p̂n,θ , p0,θ ) = Op (n−1/3 log2 n),θ

where the underlying density function is∆∆

p0,θ ≡ F0∆1i (Xi0 β; θ) × (F0 (Xi0 β + α; θ) − F0 (Xi0 β; θ)) 2i × (1 − F0 (Xi0 β + α; θ)) 3i ,

45and p̂n,θ is the corresponding maximum likelihood estimator given θ.
Since the Hellinger distance (Van de Geer (1993)) is equal to
Z 
2
1/2
2h(p̂n,θ , p0,θ ) =F̃n1/2 (u; θ) − F0 (u; θ) dQ
Z 2(F̃n (u + α; θ) − F̃n (u; θ))1/2 − (F0 (u + α; θ) − F0 (u; θ))1/2 dQ+Z 2(1 − F̃n (u + α; θ))1/2 − (1 − F0 (u + α; θ))1/2 dQ,+
it is obvious thatsup k F̃n (·; θ) − F0 (·; θ) k2 = Op (n−1/3 log2 n),θ2RF̃n (u; θ) − F0 (u; θ) g0 (u)du.
where the L2 norm is defined as k F̃n (u; θ) − F0 (u; θ) k22 ≡As a consequence of the uniform convergence for NPMLE, we haveF̃nk (α̃nk + x0 β̃nk ; α̃nk , β̃nk ) → F0 (α∗ + x0 β ∗ ; α∗ , β ∗ ).
Thereafter, the following uniform convergence is immediate.
(7.18)|Φnk (θ̃nk ) − Φ(θ∗ )| →p 0.

Given the true unknown θ0 is the unique root of the continuous limiting function Φ(·), wemust have θ∗ = θ0 , which leads to the consistency of θ̃n .
The proof of θ̃n ’s existence is more tedious (see Theorem 4.1 in Groeneboom and Hendrickx (2018)) and is relegated to Appendix C in the supplemental note.
Proof of Theorem 4.4. Denote the stacked moment conditions by
ζ(Zi ; α, β, F (; α, β)) =[∆1i − F (Xi0 β; α, β)]Xi∆3i − 1 + F (α + Xi0 β; α, β)!

We set the value of Φn to be zero at the point of zero-crossing θ̃n = (α̃n , β̃n ) asΦn (θ̃n ) ≡ Pn ζ(Z; α̃n , β̃n , F̃n (·; α̃n , β̃n )) = 0.
With this definition, we use the representation of the components as a convex combinationof the left and right limit at θ̃n :(7.19)

Φn,j (θ̃n ) = λj Φn,j (θ̃n −) + (1 − λj )Φn,j (θ̃n +) = 0,

where Φn,j denotes the jth coordinate of Φn and where we can choose λj from the unit
interval in such a way that (7.19) holds since we have a crossing of zero component-wise.
Note that this does not change the location of zero-crossing point either.

46

Then we proceed with
0 =(Pn − P )ζ(Z; α̃n , β̃n , F̃n (·; α̃n , β̃n ))+P [ζ(Z; α0 , β0 , F̃n (·; α0 , β0 )) − ζ(Z; α0 , β0 , F0 (·))](7.20)+P [ζ(Z; α̃n , β̃n , F̃n (·; α̃n , β̃n )) − ζ(Z; α0 , β0 , F̃n (·; α0 , β0 ))]In the technical lemma (see Lemma S2.3) of our Appendix C, we prove that by the consistency of our estimator and the P-Donsker property of the corresponding functional class,one has(Pn − P )ζ(Z; α̃n , β̃n , F̃n (·; α̃n , β̃n )) = Pn ζ(Z; α0 , β0 , F0 ) + op (n−1/2 ).
Also, by applying the P-Donsker property and taking a Taylor expansion, one has
P [ζ(Z; α̃n , β̃n , F̃n (·; α̃n , β̃n )) − ζ(Z; α0 , β0 , F̃n (·; α0 , β0 ))]= P [ζ(Z; α̃n , β̃n , F0 (·; α̃n , β̃n )) − ζ(Z; α0 , β0 , F0 (·))] + op (n−1/2 )!
α̃n − α0+ op (n−1/2 + (α̃n − α0 ) + |β̃n − β0 |).
= H0β̃n − β0
When it comes to the second term in (7.20), we integrate by parts to getZ
P [ζ(Z; α0 , β0 , F̃n (·; α0 , β0 )) − ζ(Z; α0 , β0 , F0 (·))] = c(u)d[(F̃n − F0 )(u)] = κ(F̃n ) − κ(F0 ),which reduces the problem to characterizing the asymptotic property of the linear functional
for the NPMLE. Thus, we proceed in Appendix C to show thatP [ζ(Z; α0 , β0 , F̃n (·; α0 , β0 )) − ζ(Z; α0 , β0 , F0 (·))] = (Pn − P )φF0 + op (n−1/2 ).
In sum, we obtain!
α̃n − α0
H0
= Pn ζ(Z; α0 , β0 , F0 ) + (Pn − P )φF0 + op (n−1/2 + (α̃n − α0 ) + |β̃n − β0 |).
β̃n − β0Hence, the desired conclusion follows given the short-hand notation φ0 ≡ ζ(Z; α0 , β0 , F0 ).

Proof of Theorem 4.5. The bootstrap validity of β̂n∗ has been shown in Groeneboom andHendrickx (2017)[p3465, equation (4.19)]. Here we focus on α̂n∗ . To prove its conditional
weak convergence, we start with the bootstrap estimating equation and we decompose it
into
n
h
i
1X
∗
∗
∗
(7.21)+ I3n,
Mni 1 − F̂n∗ (Xi0 β̂n∗ + α̂n∗ ; β̂n∗ ) − ∆3i = I1n
+ I2n
n i=1

47

where
(7.22)

1∗I1n=(7.23)

1∗=I2n(7.24)

1
∗
=
I3n

n
X

n i=1
n
X

n i=1
n
X

n i=1

Mni [1 − F0 (Xi0 β0 + α0 ) − ∆3i ] ,MnihF0 (Xi0 β0 + α0 ) − F̂n∗ (Xi0 β0 + α0 ; β0 )

i

Mni

h

F̂n∗ (Xi0 β0 + α0 ; β0 ) − F̂n∗ (Xi0 β̂n∗ + α̂n∗ ; β̂n∗ ),i.

The general scheme is in analog with our proof of Theorem 4.2. First of all, note that∗
= OPM (n−1/2 ) in PZ -probability.
I1n
∗b
∗a
∗
where
+ I2n
, we get I2n = I2nReferring to I2n(7.25) h
i
h
i
∗a
∗b
I2n
= P F0 (U + α0 ) − F̂n∗ (U + α0 ; β0 )and I2n= (P∗n −P ) F0 (U + α0 ) − F̂n∗ (U + α0 ; β0 ) .
b
We shall utilize P -Donsker property (Van Der Vaart and Wellner (1996)) to show I2n
=
−1/2
op (n
) in Lemma S1.4 in the supplemental note. We also state the linear representation∗aof I2n in Lemma S1.9.
∗c
∗b
∗a
∗
∗
where
+ I3n
+ I3n
= I3n
, we decompose it into three terms I3n
When it comes to I3n(7.26)ih∗aI3n= P F0 (X 0 β̂n∗ + α̂n∗ ; β̂n∗ ) − F0 (X 0 β0 + α0 ) ,
h
i
∗b
I3n
= P F̂n∗ (X 0 β̂n∗ + α̂n∗ ; β̂n∗ ) − F̂n∗ (X 0 β0 + α0 ; β0 ) − F0 (X 0 β̂n∗ + α̂n∗ ; β̂n∗ ) + F0 (X 0 β0 + α0 ) ,
h
i
∗c
∗
∗
0 ∗
∗
∗
∗
0
I3n = (Pn − P ) F̂n (X β̂n + α̂n ; β̂n ) − F̂n (X β0 + α0 ; β0 ) .
∗b∗b== op (n−1/2 ), I3n
In Lemma S1.11 and Lemma S1.4 respectively, we prove that I3n
op (n−1/2 ) using P -Donsker Property of related functional classes. Also we have the following
expansion:(7.27)∗aI3n= Vα0 (α̂n∗ − α0 ) + Vβ0 (β̂n∗ − β0 ) + op (n−1/2 + α̂n∗ − α0 + |β̂n∗ − β0 |).

By taking the difference of two linear representation for α̂n∗ and α̂n respectively, we get√(7.28)G∗n [ψ0 + ψF0 + ψβ0 ] + op (1),n (α̂n∗ − α̂n ) = Vα−10and the claim follows from Theorem 3.6.13 in Van Der Vaart and Wellner (1996).

Proof of Theorem 4.6. The overall structure of the proof is similar as the one for Theorem
4.4. The only necessary change is that one has to apply the maximal inequality withmultiplier bootstrap weights to relevant functional classes. To avoid repetition, we only
sketch the main steps. We skip the steps leading to the consistency of (α̃n∗ , β̃n∗ ) and directly

48

start withP∗n ζ(Z; α̃n∗ , β̃n∗ , F̃n∗ (·; α̃n∗ , β̃n∗ )) = 0.
Then we proceed with
0 =(P∗n − P )ζ(Z; α̃n∗ , β̃n∗ , F̃n∗ (·; α̃n∗ , β̃n∗ ))+P [ζ(Z; α0 , β0 , F̃n∗ (·; α0 , β0 )) − ζ(Z; α0 , β0 , F0 (·))]+P [ζ(Z; α̃n∗ , β̃n∗ , F̃n∗ (·; α̃n∗ , β̃n∗ )) − ζ(Z; α0 , β0 , F̃n∗ (·; α0 , β0 ))]Then one can show that(P∗n − P )ζ(Z; α̃n∗ , β̃n∗ , F̃n∗ (·; α̃n∗ , β̃n∗ )) = P∗n ζ(Z; α0 , β0 , F0 ) + op (n−1/2 ),P [ζ(Z; α0 , β0 , F̃n∗ (·; α0 , β0 )) − ζ(Z; α0 , β0 , F0 (·))] = (P∗n − P )φF0 + op (n−1/2 ),andP [ζ(Z; α̃n∗ , β̃n∗ , F̃n∗ (·; α̃n∗ , β̃n∗ )) − ζ(Z; α0 , β0 , F̃n∗ (·; α0 , β0 ))]= P [ζ(Z; α̃n∗ , β̃n∗ , F0 (·; α̃n∗ , β̃n∗ )) − ζ(Z; α0 , β0 , F0 (·))] + op (n−1/2 )!
α̃n∗ − α0= H0+ op (n−1/2 + (α̃n∗ − α0 ) + |β̃n∗ − β0 |).
∗β̃n − β0In the end, we get!
α̃n∗ − α0= P∗n ζ(Z; α0 , β0 , F0 ) + (P∗n − P )φF0 + op (n−1/2 + (α̃n∗ − α0 ) + |β̃n∗ − β0 |),H0β̃n∗ − β0
which leads toα̃n∗ − α̃nβ̃n∗ − β̃n!
= H0−1 [(P∗n − Pn )ζ(Z; α0 , β0 , F0 ) + (P∗n − Pn )φF0 ] + op (n−1/2 ).


References
Abrevaya, J., and J. Huang (2005): “On the bootstrap of the maximum score estimator,” Econometrica, 73(4), 1175–1204.
Aitchison, J., and S. Silvey (1957): “The generalization of probit analysis to the caseof multiple responses,” Biometrika, 44, 131–140.
An, M. Y., B. J. Christensen, and N. D. Gupta (2004): “Multivariate mixed proportional hazard modelling of the joint retirement of married couples,” Journal of Applied
Econometrics, 19(6), 687–704.

49

Attanasio, O. P., G. P. Koujianou, and E. Kyriazidou (2008): “Credit constraints
in the market for consumer durables: Evidence from micro data on car loans,” International Economic Review, 49(2), 401–436.
Ayer, M., H. Brunk, G. Ewing, W. Reid, and E. Silverman (1955): “An empirical
distribution function for sampling with incomplete information,” Annals of Mathematical
Statistics, 26, 641–647.
Banerjee, M., D. Mukherjee, and S. Mishra (2009): “Semiparametric binary regression models under shape constraints with an application to Indian schooling data,”
Journal of Econometrics, 149(2), 101–117.
Bellemare, C., B. Melenberg, and A. Van Soest (2002): “Semi-parametric models
for satisfaction with income,” Portuguese Economic Journal, 1, 181–203.
Bresnahan, T. F., and P. C. Reiss (1991): “Entry and competition in concentrated
markets,” Journal of Political Economy, 99, 977–1009.
Cameron, S. V., and J. J. Heckman (1998): “Life cyce schooling and dynamic selection
bias: Models and evidence for five cohorts of American males,” Journal of Political
Economy, 106, 262–333.
Cavanagh, C., and R. P. Sherman (1998): “Rank estimators for monotonic index
models,” Journal of Econometrics, 84, 351–381.
Chen, S. (2002): “Rank estimation of transformation models,” Econometrica, 70, 1683–
1697.
Chen, S., and S. Khan (2003): “Rates of convergence for estimating regression coefficients in heteroskedastic discrete response models,” Journal of Econometrics, 117,
245–278.
Chen, X., O. Linton, and I. Van Keilegom (2003): “Estimation of semiparametricmodels when the criterion function is not smooth,” Econometrica, 71, 1591–1608.
Cheng, G., and J. Huang (2010): “Bootstrap consistency for general semiparametric
M-estimate,” Annals of Statistics, 38, 2884–2915.
Coppejans, M. (2007): “On efficient estimation of the ordered response model,” Journalof Econometrics, 137, 577–614.
Cosslett, S. R. (1983): “Distribution-free maximum likelihood estimator of the binary
choice model,” Econometrica, 51, 765–782.
Cunha, F., J. Heckman, and S. Navarro (2007): “The identification and economic
content of ordered choice models with stochastic thresholds,” International Economic
Review, 48, 1273–1309.
de Paula, Á. (2009): “Inference in a synchronization game with social interactions,”
Journal of Econometrics, 148, 56–71.

50

Durlauf, S. N., and Y. M. Ioannides (2010): “Social interactions,” Annual Review of
Economics, 2, 451–478.
Efron, B. (1979): “Bootstrap methods: Another look at the jackknife,” The Annals of
statistics, 7, 1–26.
Gentleman, R., and A. Vandal (2018): “Icens: NPMLE for censored and truncateddata,” R package version 1.53.0.
Geskus, R., and P. Groeneboom (1996): “Asymptotically optimal estimation of smooth functionals for interval censoring, I,” Statistica Neerlandica, 50, 69–88.
(1997): “Asymptotically optimal estimation of smooth functionals for interval
censoring, II,” Statistica Neerlandica, 51, 201–219.
(1999): “Asymptotically optimal estimation of smooth functionals for interval
censoring, case 2,” The Annals of Statistics, 27, 627–674.
Greene, W. H., and D. A. Hensher (2010): Modeling ordered choices: A primer.
Cambridge University Press.
Groeneboom, P. (2014): “Maximum smoothed likelihood estimators for the interval
censoring model,” The Annals of Statistics, 42, 2092–2137.
Groeneboom, P., and K. Hendrickx (2017): “The nonparametric bootstrap for thecurrent status model,” Electronic Journal of Statistics, 11, 3446–3483.
(2018): “Current status linear regression,” The Annals of Statistics, 46, 1415–
1444.
Groeneboom, P., and G. Jongbloed (2014): Nonparametric estimation under shape
constraints. Cambridge University Press.
Groeneboom, P., G. Jongbloed, and B. I. Witte (2010): “Maximum smoothed
likelihood estimation and smoothed maximum likelihood estimation in the current status
model,” The Annals of Statistics, 38, 352–387.
Groeneboom, P., and J. A. Wellner (1992): Information Bounds and NonparametricMaximum Likelihood Estimation. Birkhauser.
Gustman, A. L., and T. L. Steinmeier (2000): “Retirement in dual-career families: astructural model,” Journal of Labor economics, 18, 503–545.
(2004): “Social security, pensions and retirement behaviour within the family,”
Journal of Applied Econometrics, 19, 723–737.
Hajivassiliou, V. A., and Y. M. Ioannides (2007): “Unemployment and liquidity
constraints,” Journal of Applied Econometrics, 22, 479–510.
Hong, H., J. D. Kubik, and J. C. Stein (2004): “Social interaction and stock-market
participation,” The Journal of Finance, 59(1), 137–163.
Honoré, B. E., and Á. de Paula (2010): “Interdependent Durations,” Review of Economic Studies, 77, 1138–1163.

51(2018): “A new model for interdependent durations with an application to joint
retirement,” Quantitative Economics, 9, 1299–1333.
Horowitz, J. L. (1992): “A smoothed maximum score estimator for the binary responsemodel,” Econometrica, pp. 505–531.
(1996): “Semiparametric Estimation of a Regression Model with an Unknown
Transformation of the Dependent Variable,” Econometrica, 64, 103–137.
(2009): Semiparametric and nonparametric methods in econometrics. Springer.
Huang, J., and J. A. Wellner (1995): “Asymptotic normality of the NPMLE of linear
functionals for interval censored data, Case 1,” Statistica Neerlandica, pp. 153–163.
Ichimura, H. (1993): “Semiparametric least squares (SLS) and weighted SLS estimationof single-index models,” Journal of Econometrics, 58.
Katz, M. L., and C. Shapiro (1986): “Technology adoption in the presence of network
externalities,” Journal of Political Economy, 94(4), 822–841.
Kiefer, J., and J. Wolfowitz (1956): “Consistency of the maximum likelihood estimator in the presence of infinitely many incidental parameters,” The Annals of Mathematical
Statistics, 27(4), 887–906.
Klein, R. W., and R. P. Sherman (2002): “Shift restrictions and semiparametricestimation in ordered response models,” Econometrica, 70(2), 663–691.
Klein, R. W., and R. H. Spady (1993): “An efficient semiparametric estimator forbinary response models,” Econometrica, 61(2), 387–421.
Lee, M.-j. (1992): “Median regression for ordered discrete response,” Journal of Econometrics, 51, 59–77.
Lewbel, A. (2000): “Semiparametric qualitative response model estimation with unknown
heteroscedasticity or instrumental variables,” Journal of Econometrics, 97, 145–177.
(2002): “Ordered response threshold estimation,” Working paper.
Lewbel, A., and S. M. Schennach (2007): “A simple ordered data estimator for inverse
density weighted expectations,” Journal of Econometrics, 136(1), 189–211.
Lundberg, S., and R. A. Pollak (1994): “Noncooperative bargaining models of marriage,” The American Economic Review, 84, 132–137.
Manski, C. F. (1985): “Semiparametric analysis of discrete response: Asymptotic properties of the maximum score estimator,” Journal of Econometrics, 27, 313–333.
Manski, C. F., and E. Tamer (2002): “Inference on regressions with interval data on aregressor or outcome,” Econometrica, 70(2), 519–546.
Matzkin, R. L. (1991): “Semiparametric estimation of monotone and concave utility
functions for polychotomous choice models,” Econometrica, 59, 1315–1327.
(1993): “Nonparametric identification and estimation of polychotomous choicemodels,” Journal of Econometrics, 58, 137–168.

52

Melenberg, B., and A. Van Soest (1996): “Measuring the costs of children: parametric and semiparametric estimators 1,” Statistica Neerlandica, 50, 171–192.
Nan, B., and J. A. Wellner (2013): “A general semiparametric Z-estimation approachfor case-cohort studies,” Statistica Sinica, 23(3), 1155.
Newey, W. K. (1994): “The asymptotic variance of semiparametric estimators,” Econometrica, 62, 1349–1382.
Park, A., and L. Smith (2008): “Caller number five and related timing games,” Theoretical Economics, 3, 231–256.
Prescott, E., and M. Visscher (1977): “Sequential location among firms with foresight,” Bell Journal of Economics, 8, 378–893.
Ridder, G. (1990): “The non-parametric identification of generalized accelerated failuretime models,” Review of Economic Studies, 57, 167–181.
Rubin, D. (1981): “Bayesian bootstrap,” The Annals of statistics, 9, 130–134.
Schick, A., and Q. Yu (2000): “Consistency of the GMLE with mixed case intervalcensored data,” Scandinavian Journal of Statistics, 27, 45–55.
Shaked, A., and J. Sutton (1982): “Relaxing price competition through product differentiation,” Review of Economic Studies, 49, 3–13.
Stewart, M. B. (2005): “A comparison of semiparametric estimators for the ordered
response model,” Computational statistics & data analysis, 49, 555–573.
Tanaka, H. (2008): “Semiparametric least squares estimation of monotone single index
models and its application to the iterative least squares estimation of binary choicemodels,” Working Paper.
Van de Geer, S. (1993): “Hellinger-consistency of certain nonparametric maximum likelihood estimators,” The Annals of Statistics, 21, 14–44.
(1997): “Asymptotic normality in mixture models,” ESAIM: Probability andStatistics, 1, 17–33.
(2000): Empirical processes in M-estimation. Cambridge University Press.
Van der Vaart, A. (1991): “On differentiable functionals,” The Annals of Statistics, 19,
178–204.
Van Der Vaart, A. (1998): Asymptotic statistics. Cambridge University Press.
Van Der Vaart, A., and J. A. Wellner (1996): Weak convergence and empirical
processes. Springer.
Vella, F. (1998): “Estimating models with sample selection bias: a survey,” Journal of
Human Resources, 33, 127–169.
Wellner, J. A., and Y. Zhan (1997): “A hybrid algorithm for computation of the nonparametric maximum likelihood estimator from censored data,” Journal of the American
Statistical Association, 92, 945–959.

53

Wu, C. (1990): “On the asymptotic properties of the jackknife histogram,” The Annalsof statistics, 18, 1438–1452.

Supplemental Note for “Simple Semiparametric
Estimation of Ordered Response Models: With an
Application to the Interdependent Durations
Model”

Ruixuan Liu and Zhengfei Yu
Emory University and University of Tsukuba

S1. Appendix B: Technical Proofs Related to Two-stage
Estimation
We first restate some necessary definitions and Theorem 2.4.1 in Van Der Vaart andWellner (1996) that will be used repeatedly in the sequel. Let F be the class of functionsand L2 (Q) be the L2 -norm defined by a probability measure Q . For any probability
measure Q, let N (ε, F, L2 (Q)) be the minimal number of balls of radius ε needed to cover
the class F. The entropy integral J(δ, F) is defined as
Z δpJ(δ, F) ≡ sup
1 + log N (ε, F, L2 (Q))dε.
Q0

An envelope function of a functional class F is a function F such that |f (x)| ≤ F (x) for
all x and f ∈ F.
Lemma S1.1 (Theorem 2.14.1 in (Van Der Vaart and Wellner (1996))). Let P0 be the
distribution of the underlying observation and let F be a P0 -measurable class with anenvelope function F . We have(S1.1)

E sup |Gn f | . J(1, F) k F kP0 ,2
f ∈F

Lemma S1.2. [Lemma 3.6.7 in (Van Der Vaart and Wellner (1996))] Let Zn1 , ..., Znn
be arbitrary stochastic processes and (Mn1 , ..., Mnn )0 be any exchangeable random vector12

independent of Zn1 , ..., Znn . For any n0 > 0 and n > n0 , we have
!


n
X
1
EM (max1≤i≤n |Mni |)√√EZM≤ n0 EZ (kZn1 k)
Mni Zni
n i=1
n

Z ∞
p
PM (Mn1 ≥ u)du EZ
max+
0n0 <i≤n

i
1 X
√
Znj
n j=n +1!
.

0

We need to apply the following well-known entropy bounds concerning monotone functions or functions of bounded variation repeatedly. The bounds actually hold for thebracketing entropy uniformly over the underlying probability measure, which will be usedin Appendix C as well. We refer readers to Theorem 2.7.5 on [p.159] of Van Der Vaart andWellner (1996) or Lemma 3.8 on [p.36] of Van de Geer (2000) for the proofs.
Lemma S1.3 (Entropy Bounds). Let AC be the class of monotone functions with values
in [0, C], then for all δ > 0,√
(S1.2)J(δ, AC ) . δ.
Let BC be the class of functions of bounded variation with values in [0, C], then for all
δ > 0,√(S1.3)J(δ, BC ) . δ.
Now we obtain the entropy bounds for the key functional class in our context and prove
the asymptotic characterizations for terms appearing in our Theorem 4.1 and Theorem 4.2.
Lemma S1.4. The functional class G defined byno(S1.4)G ≡ (x, δ3 ) 7→ (1 − δ3 − F (x0 β + α)) : (α, β) ∈ Θ, F (·) ∈ A
has bounded entropy integral. Therefore, we have the following Glivenko-Cantelli resultsh i(Pn − P ) F̂n X 0 β̂n + α̂n ; β̂n = op (1), andih (P∗n − P ) F̂n X 0 β̂n∗ + α̂n∗ ; β̂n∗ = oPM (1),in PZ -probability. Moreover, we obtain the stochastic equicontinuity ash i(S1.5)Gn F̂n X 0 β̂n + α̂n ; β̂n − F0 (X 0 β0 + α0 ) = op (1), and(S1.6)h iG∗n F̂n X 0 β̂n∗ + α̂n∗ ; β̂n∗ − F0 (X 0 β0 + α0 ) = oPM (1),in PZ -probability.
Proof. We first verify that the uniform entropy integral J(1, G) is bounded. Because theisotonic estimator F̂n (t, β) is a monotonically increasing function for any given β, G is the3class of composite functions involving an monotonically increasing link/ridge function and alinear index x0 β +α with parameters (α, β) belonging to a compact Euclidean space. Hence,
by Lemma 2.3 in Baladbaoui, Groeneboom, and Hendrickx (2017) we get the following
bound on the uniform entropy log N (ε, G) . 1/ε, so the uniform entropy integral J(1, G)
is indeed bounded. Therefore, the functional class G is P-Donsker, which directly implies
the stated Glivenko-Cantelli properties.
Regarding the stochastic equicontinuity, consider the following class:noG ≡ x 7→ (F (x0 β + α) − F0 (x0 β0 + α0 )) : (α, β) ∈ Θ, F (·) ∈ A, |α−α0 |∨ k β−β0 k ∨ k F −F0 k∞ ≤  ,for some small positive
. Again G has bounded entropy integral similarly as G. Moreover,0F̂n X β̂n + α̂n ; β̂n − F0 (X 0 β0 + α0 ) belongs to G with probability tending to 1, because(S1.7)
k F̂n X 0 β̂n + α̂n ; β̂n − F0 (X 0 β0 + α0 ) k∞






≤k F̂n X 0 β̂n + α̂n ; β̂n − F0 X 0 β̂n + α̂n ; β̂n k∞ + k F0 X 0 β̂n + α̂n ; β̂n − F0 (X 0 β0 + α0 ) k∞
→ 0.
The first term on the right hand side of the inequality follows from the uniform consistency of NPMLE as in (S1.16), whereas the convergence of the second term is due to thesmoothness of F0 (u; β) (w.r.t. both u and β) and consistency of α̂n and β̂n . Thereafter,
the desired stochastic equicontinuity follows from applying (S1.1) to the class G .
When it comes to the bootstrap version, applying the multiplier inequality in Lemma S1.2
to G , we get(S1.8)

1
EZM |k G∗n k| . EZ |G | √ EM max Mni + EZ max k Gk k ,in0 ≤k≤nn

where G is the corresponding envelope function. The first term is of smaller order since1√ EM max Mni = op (1),inunder our assumptions on the bootstrap weights. Meanwhile, the Levy inequality (Proposition A.1.2 of Van Der Vaart and Wellner (1996)) implies:
P {max k Gk k> λ} ≤ 2P {k Gn k> λ}k≤nwhich makes the second term negligible.

An immediate consequence of Lemma S1.4 is Lemma S1.5 showing the negligibility ofcc∗related terms in our Theorem 4.2 and Theorem 4.5. The claims on I3nand I3ndirectly

4bb∗follow from results (S1.5) and (S1.6). Regarding I2n
and I2n, the proofs are even easier
because both α0 and β0 are fixed in those two terms.

Lemma S1.5. Suppose Conditions (1)-(8) hold. We have the characterization of following
smaller order terms:√ b√ c(S1.9)
nI2n = op (1)nI3n = op (1), and
√ c∗
√ b∗
(S1.10)nI2n = oPM (1)nI3n = oPM (1),in PZ -probability.
a
Now we prove several preparatory lemmas related to the linear representation of I2n.

Lemma S1.6. Suppose Conditions (1)-(8) hold. The following representations hold:Za(S1.11)I2n = − φα0 (u)(F̂n (u; β0 ) − δ1 )dP (u, δ1 ),where(S1.12)

φα0 (u) = g0 (u − α0 )/g0 (u).

Proof. The result follows similar argument as in Lemma 4.1 of Groeneboom, Jongbloed,
and Witte (2010):
Z a(S1.13)I2n = −
F̂n (u + α0 ; β0 ) − F0 (u + α0 ) g0 (u)duZ 
=−F̂n (u; β0 ) − F0 (u) g0 (u − α0 )duZ= − φα0 (u) F̂n (u; β0 ) − F0 (u) dG0 (u)Z= − φα0 (u)(F̂n (u; β0 ) − δ1 )dP (u, δ1 ),where in the last equality we have used the fact that δ1 dP = F0 (u)g0 (u), since the probability density function of the binary choice data (U, ∆1 ) is(S1.14)p(u, δ1 ) = F0 (u)δ1 (1 − F0 (u))1−δ1 g0 (u).


We consider the piece-wise constant version of φF0 which is constant on the same intervals
where the NPMLE F̂n (·; β) remains constant. Denote those intervals by [τi , τi+1 ). We define(S1.15)φ̄α0 (u) = φα0 (Ân (u; β)),

5

where
Ân (u; β) =




τi ,

if ∀t ∈ Ji : F0 (t) > F̂n (τi ; β),s,

if ∃s ∈ Ji : F0 (s) = F̂n (s; β),τi+1 ,

if ∀t ∈ Ji : F0 (t) < F̂n (τi ; β),for u ∈ Ji .
The following uniform convergence results are available in Lemma 5.9 in Groeneboom andWellner (1992), Lemma 3.1 of Groeneboom and Hendrickx (2018), and the corollary ofLemma 3.1 in Groeneboom and Hendrickx (2017).
Lemma S1.7. Suppose Conditions (1)-(8) hold, then we have
(S1.16)
P lim sup F̂n (u; β) − F0 (u; β) = 0 = 1,n→∞ β∈B,uand
sup F̂n (u; β) − F0 (u; β) = Op log n × n−1/3 ,β∈B,usup Ân (u; β) − u = Op log n × n−1/3 .
β∈B,uFurthermore, for the bootstrap version, we have∗
(S1.17)PM lim sup F̂n (u; β) − F0 (u; β) = 0|Z = 1,n→∞ β∈B,uand
sup F̂n∗ (u; β) − F0 (u; β) = Op log n × n−1/3 ,

β∈B,u


sup Â∗n (u; β) − u = Op log n × n−1/3 .

β∈B,uin PZ -probability.
Lemma S1.8. Suppose Conditions (1)-(8) hold, then we have(S1.18)k φ̄α0 − φα0 k∞ .k F̂n (u; β0 ) − F0 (u) k∞ .

Proof. Our proof adapts the argument in Lemma A.4 of Groeneboom, Jongbloed, andWitte (2010). Let ξ to denote some intermediate value in Taylor expansions whose value
may change from line to line. For given u in the support of F0 (·; β), then the interval Ji it
belongs to is of one of the following three types:(i) F0 (x; β) > F̂n (τi ; β) for all x ∈ Ji ;(ii) F0 (x; β) > F̂n (x; β) for some x ∈ Ji ;(iii) F0 (x; β) < F̂n (τi ; β) for all x ∈ Ji ;

6

see Figure 5 in Groeneboom, Jongbloed, and Witte (2010). First of all, when F0 (u; β) =F̂n (u; β), then by definition of φ̄α0 (u) = φα0 (u) so both the left-hand side and right-handside of (S1.18) are zero.
Next, for the case where F0 (u; β) 6= F̂n (u; β). For v, ξ ∈ Ji , we get by Taylor expansion:(S1.19) |F̂n (u; β)−F0 (u; β)| = |F̂n (v; β)−F0 (u; β)| = |F̂n (v; β)−F0 (v; β)−(u−v)f0 (ξ; β)|.
There are three possibilities. If Ân (u; β) = τi , then F0 (τi ; β) − F̂n (τi ; β) > 0 which gives
rise to(S1.20) |F̂n (u; β) − F0 (u; β)| = |(u − τi )f0 (ξ; β) + F0 (τi ; β) − F̂n (τi ; β)| ≥ |u − τi |f0 (ξ; β),where we apply (S1.19) while letting v = τi .
If Ân (u; β) = v for some v 6= u ∈ Ji , then F̂n (v; β) = F0 (v; β) so(S1.21) |F̂n (u; β) − F0 (u; β)| = |F̂n (v; β) − F0 (v; β) − (u − v)f0 (ξ; β)| = |u − v|f0 (ξ; β).
If Ân (u; β) = τi+1 , then F̂n (τi=1 −; β) > F0 (τi+1 −; β) ≥ 0 giving us(S1.22)|F̂n (u; β) − F0 (u; β)| = |(τi+1 − u)f0 (ξ; β) + F̂n (τi+1 −; β) − F0 (τi+1 ; β)| ≥ |τi+1 − u|f0 (ξ; β).
In sum, we obtain|F̂n (u; β) − F0 (u; β)| ≥ |v − u|f0 (ξ; β) ≥ c|u − v|,(S1.23)where we use the fact that the density function is bounded below by a positive constantover the compact interval.
Finally, the imposed smoothness condition on φα0 leads to|φ̄α0 (u) − φα0 (u)| = |φα0 (v) − φα0 (u)| ≤ c|v − u|,(S1.24)

where v ∈ [τi , τi+1 ]. Hence, the desired claim follows by combining (S1.23) and (S1.24). aGiven the above lemmas, we get the following characterization of I2n
and its bootstrapped∗aversion I2n .

Lemma S1.9. Suppose Conditions (1)-(8) hold, then we have the following linear representations:
h
i
√ a
√
nI2n = nP F0 − F̂n (·, β0 ) = Gn ψF0 + op (1).
(S1.25)and(S1.26)

√

∗a
nI2n
=

in PZ -probability.

h
i √
√
∗
nP F0 − F̂n (·, β0 ) = n [P∗n − P ] ψF0 + opM (1),7aProof. We only proof the claim regarding I2n
to avoid repetition. Given the characterizationof our isotonic estimate F̂n (u; β0 ) and the piece-wise constant nature of φ̄α0 , we get
Z(S1.27)φ̄α0 [F̂n (u; β0 ) − δ1 ]dPn = 0,

by equality (8.15) in Groeneboom and Jongbloed (2014). Therefore, starting with theaalternative representation of I2n
we getZ
a(S1.28)= φ̄α0 (F̂n (u; β0 ) − δ1 )d(Pn − P )I2nZ(S1.29)+ [φ̄α0 − φα0 ](F̂n (u; β0 ) − δ1 )dP (u, δ1 ).
In the next lemma, we show thatZ
Z
φ̄α0 (F̂n (u; β0 ) − δ1 )d(Pn − P ) = φα0 (F0 (u) − δ1 )d(Pn − P ) + op (n1/2 ),(S1.30)andZ(S1.31)[φ̄α0 − φα0 ](F̂n (u; β0 ) − δ1 )dP (u, δ1 ) = Op (n−2/3 ),which lead to the desired conclusion.

Lemma S1.10. Suppose Conditions (1)-(8) hold, then the following hold:Z
Rn ≡ φ̄α0 (F̂n (u; β0 ) − F0 (u))d(Pn − P ) = op (n−1/2 ),andZ(S1.32)Sn ≡[φ̄α0 − φα0 ](F̂n (u; β0 ) − δ1 )dP (u, δ1 ) = op (n−1/2 ).

Proof. We first handle the term Sn as follows.
Z(S1.33)
Sn = [φ̄α0 − φα0 ](F̂n (u; β0 ) − F0 (u))dG(u)(S1.34).k F̂n (u; β0 ) − F0 (u) k2∞ = Op (n−2/3 × log2 n),

where in the second step we used (S1.8).
Referring to the term Rn , we introduce some notations adapted from Groeneboom, Jongbloed, and Witte (2010). Define(S1.35)ξB (u) = φ̄α0 (u)B(u)where B is a function of bounded variation and with bounded superior norm C. And let(S1.36)

GC ≡ {ξB (u) : B ∈ BC }.

8

Given the result in (S1.7), for any small γ > 0 we can find finite constant term C such thatfor all n sufficiently large:
P r{Υn,C } ≡ P r{sup |F̂n (u; β) − F0 (u; β)| ≤ Cn−1/3 log n}u,β≥ 1 − γ/2.
Now for the vanishing sequence νn to be specified later, it is straightforward to arrive at
P r{|n1/2 Rn | > νn } = P r{|n1/2 Rn | > νn ∩ Υn,C } + P r{|n1/2 Rn | > νn ∩ Υcn,C }≤ νn−1 E |n1/2 Rn |1{Υn,C } + γ/2,
for any small γ. Again by (S1.7), we have
E |n1/2 Rn |1{Υn,C } ≤ E sup n1/2−1/3 log nZ
φ̄α0 (u)B(u)d(Pn − P )

B∈BC

≤n

−1/3

Z
log nE sup

ξ(u)dGn (u) .

ξ∈GC

The rest of our proof is to utilize Theorem 2.14.1 in Van Der Vaart and Wellner (1996)
to bound the expectation in the last display. Following the construction in Groeneboom,
Jongbloed, and Witte (2010), the entropy integral of GC is bounded above by a finite
constant, i.e., J(1, GC ) < ∞. And the L2 -norm of the envelope function is also bounded.
Thus applying (S1.1), we haveZ
−5/6
ξ(u)dGn (u) . n−5/6 log n,(S1.37)
E|Rn | ≤ n
log n × E supξ∈GC

which immediately leads to Rn = op (n−1/2 ).
Lemma S1.11. Suppose Conditions (1)-(8) hold. We have the characterization of following
smaller order terms:√ b√ b∗(S1.38)nI3n = oP (1) andnI3n = oPM (1),in PZ -probability.

b
Proof. Again we only prove the part involving I3n
as the other one follows analogously.
First of all, recall that
h
i
h
i
b
0
0
0
0
I3n = P F̂n (X β̂n + α̂n ; β̂n ) − F0 (X β̂n + α̂n ; β̂n ) −P F̂n (X β0 + α0 ; β0 ) − F0 (X β0 + α0 ) .

9Following the arguments in Lemma S1.9, we get(S1.39)
b
=
I3n

Z "

g0 (u − α̂n ; βˆn )[F0 (u; β̂n ) − δ1 ]g0 (u; βˆn )!


−g0 (u − α0 )[F0 (u) − δ1 ]g0 (u)

#d(Pn − P )+ op (n−1/2 ).
The smoothness assumption in Condition (5) implies that the function in the bracket of(S1.39) belongs to a P-Donsker class by Example 19.7 in Van Der Vaart (1998). Now the√ bconvergence of α̂n and β̂n leads to the desired conclusion that nI3n
= oP (1).
S2. Appendix C: Technical Proofs Related to Joint Estimation
We first prove the main theorem which shows the root-n consistency and asymptotic
normality of our estimators for the finite dimensional parameters. The overall structureof the proof is similar to the one of Theorem 4.2. Thus, we skip some straightforward
intermediate steps and highlight the changes that we make. Some technical details are
separately verified in lemmas of this section.
Related to the P -Glivenko-Cantelli or P -Donsker property, it turns out that it is moreconvenient to work with the bracketing entropy bounds. For that purpose, we collect thenecessary definitions from Van Der Vaart and Wellner (1996) as follows. The bracketing number N[] (, F, k·k2 ) for subclass F is defined to be the minimum of m such thatU
L
for ∀f ∈ F, fjL ≤ f ≤ fjU for some j, and fjU − fjL 2 ≤ . Denote, fm∃ f1L , f1U , ..., fm
H[] (, F, k·k2 ) ≡ log N[] (, F, k·k2 ).
We need a few more notations in order to characterize the consistency and rate of convergence. Let p0 be the true density function and p be a member of the class of densities,
denoted by P. We denoteq + q0q̄ =(S2.1),
2(S2.2)P̄ ≡ {q̄ : q ∈ P},(S2.3)

P̄ 1/2 ≡ {q̄ 1/2 : q̄ ∈ P̄}.

We consider a ball (measured according to the Hellinger distance) around the true density
q0 , intersected with P̄ 1/2 byP̄ 1/2 (δ) ≡ {q̄ 1/2 ∈ P̄ 1/2 : h(q̄, q0 ) ≤ δ},

10for some small positive δ. We refer to H[] (u, P̄ 1/2 (δ)) as the local entropy with bracketingand its corresponding entropy integral is given byZ δ
1/2
J[] (δ, P̄ (δ)) ≡H[] (u, P̄ 1/2 (δ))du ∨ δ.
δ 2 /213

Similarly as in Geskus and Groeneboom (1997), we define the following density functions:(S2.4)qF (u, δ1 , δ2 , δ3 ; α, β) = δ1 F (u; α, β)+δ2 {F (u+α; α, β)−F (u; α, β)}+δ3 {1−F (u+α; α, β)}.
We record Theorem 7.4 in Van de Geer (2000) which is needed to obtain the rate ofconvergence in our context.
Lemma S2.1 (Theorem 7.4 in Van de Geer (2000)). Given a upper bound Ψ(δ) for theentropy integral function in such a way that Ψ(δ)/δ 2 is a non-increasing function of δ. Thenfor a universal constant c and for δn → 0 such that√ 2nδn ≥ cΨ(δn ),then the likelihood estimator q̃n converges with rate Op (δn ). In fact, there exist constants
L0 , C0 , s.t. for all L ≥ L0
Pr {h(q̄, q0 ) ≥ Lδn } ≤ C0 exp −L2 δn2 n .
The next lemma delivers the rate of convergence for the NPMLE in terms of L2 norm
uniformly over the finite dimensional parameter.
Lemma S2.2. Regarding the convergence by L2 norm, we have(S2.5)

sup k F̃n (u; θ) − F0 (u; θ) k2 = Op log2 n × n−1/3 .
θProof. In order to obtain the rate of convergence, we need to first bound the entropy number
for the likelihood function:∆∆

F ∆1i (Xi0 β; θ) × (F (Xi0 β + α; θ) − F (Xi0 β; θ)) 2i × (1 − F (Xi0 β + α; θ)) 3i .
The only complication comes from the termnpo
F (x0 β + α; θ) − F (x0 β; θ) : (α, β, F ) .
(S2.6)FD ≡Its entropy number can be bounded in the following way; see Example 3.3(b) of Van deGeer (1993). If F (x0 β + α; θ) − F (x0 β; θ) > δ or F̄ (x0 β̄ + ᾱ; θ̄) − F̄ (x0 β̄; θ̄) > δ, thenqp00
F (x β + α; θ) − F (x β; θ) − F̄ (x0 β̄ + ᾱ; θ̄) − F̄ (x0 β̄; θ̄)
1 
F (x0 β + α; θ) − F̄ (x0 β̄ + ᾱ; θ̄) + F (x0 β; θ) − F̄ (x0 β̄; θ̄)
<√
δ

.

11

If both F (x0 β + α; θ) − F (x0 β; θ) ≤ δ or F̄ (x0 β̄ + ᾱ; θ̄) − F̄ (x0 β̄; θ̄) ≤ δ, obviously one has
q
p
√
0
0
F (x β + α; θ) − F (x β; θ) − F̄ (x0 β̄ + ᾱ; θ̄) − F̄ (x0 β̄; θ̄) ≤ 2 δ.
In sum, for any probability measure Q, we get√
N (4 δ, FD , L2 (Q)) ≤ N (δ, F0 , L2 (Q)),
where F0 ≡ {F (x0 β + α; θ) : (α, β, F )}. Compared with the calculation in Van de Geer
(1993), one needs to account for the presence of finite dimensional parameter which incurs
an additional log n factor. Therefore, one can apply Lemma S2.1 to get
sup h(qF̃n ,θ , qF0 ,θ ) = Op (log2 n × n−1/3 ).

(S2.7)

θ

Also, not thatqp 2
(F̃n − F0 ) ≤ 4
F̃n − F02q2pand (F̃n − F0 ) ≤ 4
1 − F̃n − 1 − F0 ,2

we get thatsup k F̃n (·; θ) − F0 (·; θ) k= Op (log2 n × n−1/3 ).
θNext, we fill the details related to our proofs of the consistency and asymptotic normality for our joint estimation method. Specifically, we prove the existence of zero-crossing
points and show the stochastic equicontinuity of negligible terms related to our estimating
equations.
Existence of Zero-crossing Points. Recall that θ = (α, β 0 )0 . The first coordinate of β isnormalized to be 1; i.e., the overall number of unknown parameters is equal to K. The
uniform convergence of the estimating equation leads toΦn (θ) = Φ̇θ0 (θ − θ0 ) + rn (θ),(S2.8)

where rn (θ) = op (1) + o(θ − θ0 ). We now define for h > 0,the function
Φn,h (θ) = Φ̇θ0 (θ − θ0 ) + r̃n,h (θ),(S2.9)
with(S2.10)−dr̃n,h (θ) = hZkh (u1 − α) · · · kh (uK − βK )rn (u1 , · · · , uK )du1 · · · duK ,where k(·) is a standard kernel density function supported on [−1, 1] and β 0 = (β1 , · · · , βd )0 .
Note that limh→0 r̃n,h (θ) = rn (θ).

12

We re-parameterize as follows by defining(S2.11)γ = Φ̇θ0 θ,and γ0 = Φ̇θ0 θ0 .

This gives(S2.12)Φn,h (θ) = γ − γ0 + r̃n,h (Φ̇−1
θ0 γ).

Given the result in S2.8, the mapping(S2.13)γ 7→ γ0 − r̃n (Φ̇−1
θ0 γ)

maps, for each δ > 0, the ball Bδ (γ0 ) = {γ : |γ − γ0 | ≤ δ} into Bδ/2 (γ0 ) = {γ : |γ −γ0 | ≤ δ/2} with probability approaching to 1. Therefore by Brouwer’s fixed point theorem(Groeneboom and Hendrickx (2018)), the mapping(S2.14)

γ 7→ γ0 − r̃n,h (Φ̇−1
θ0 γ),

has a fixed point which we denote by γn,h . Let θn,h ≡ Φ̇−1
θ0 γn,h , then we have(S2.15)Φn,h (θn,h ) = 0

By compactness of the parameter space, the sequence (θn,1/k )∞k=1 must have a subsequence(θn,1/kl ) with a limit point θ̄n as l → ∞. Finally, we prove as in Groeneboom and Hendrickx(2018) that Φn (θ) has a crossing of zero at θ̄n by contradiction.
Suppose that the j-th component Φjn of Φn does not have a crossing of zero at θ̄n . Then
there must be an open ball Bδ (θ̄n ) = {θ : |θ − θ̄n | < δ} of θ̄n such that Φjn has a constantsign in Bδ (θ̄n ), say Φjn (θ) ≥ c > 0 for all θ ∈ Bδ (θ̄n ) and some constant c > 0. Arguing as
Groeneboom and Hendrickx (2018), the j-th component of Φjn,h of Φn,h satisfies(S2.16)cΦjn,h (θ) ≥ ,2

for sufficiently small h and all θ ∈ Bδ (θ̄n ), which contradicting S2.15, since θn,h for h = 1/kl
belongs to Bδ (θ̄n ) for large kl .
Lemma S2.3. Under our conditions, we get(Pn − P )ζ(Z; α̃n , β̃n , F̃n (·; α̃n , β̃n )) = Pn ζ(Z; α0 , β0 , F0 ) + op (n−1/2 ),andP [ζ(Z; α̃n , β̃n , F̃n (·; α̃n , β̃n )) − ζ(Z; α0 , β0 , F̃n (·; α0 , β0 ))]= P [ζ(Z; α̃n , β̃n , F0 (·; α̃n , β̃n )) − ζ(Z; α0 , β0 , F0 (·))] + op (n−1/2 ).
Proof. The proof essentially follows from Lemma S1.4 by the stochastic equicontinuity of
the related P-Donsker classes. The only minor change applies to the functional class that

13

the NPMLE F̃n belongs to, given that the NPMLE is a sub-distribution or a defective
distribution. However, the entropy bound for the monotone functions in Lemma S1.3 does
not depend on the range of the function, as long as it is finite. Hence, the results follow. 
From now on, we drop the additional index θ in defining the NPMLE or its probability limit with out confusion. Because the remaining proofs reduce to characterizing theasymptotic property of the linear functional for the NPMLE given the data and the true
unknown parameter θ0 . Also, we denote the empirical probability measure of the ordered
response data by Qn and its population version by QF0 where the distribution is set to be
the true unknown F0 .
We mentioned in the main text that the NPMLE could be a defective distribution (or asub-distribution) in finite sample; i.e., F̃n (u) < 1 for any u in the support. However, thisplays a minor role regarding the large sample properties we study because the defectivenessdoes not occur with probability 1 as sample size goes to infinity.
Lemma S2.4 (Proposition 1 in Geskus and Groeneboom (1997)). We havelim Pr{F̃nn→∞

is defective} = 0.

We need the following lemma which characterizes the NPMLE F̃n borrowed from Corollary 1 in Geskus and Groeneboom (1997).
Lemma S2.5. Any function σ that is constant at the same intervals as F̃n satisfiesZδ11 − δ1 − δ2
δ2
(S2.17)σ(u)−+dQn (u, δ1 , δ2 ) = 0.
F̃n (u) 1 − F̃n (u + α0 ) F̃n (u + α0 ) − F̃n (u)We also need to consider the piece-wise constant version of ξF which is constant on thesame intervals where the NPMLE F̃n (·; α, β) remains constant. Denote those intervals by[τi , τi+1 ). We define¯ = ξ (Ãn (u; α, β)),ξ(u)F̃(S2.18)whereÃn (u; α, β) =




τi ,

if ∀t ∈ Ji : F0 (t) > F̃n (τi ; α, β),s,

if ∃s ∈ Ji : F0 (s) = F̃n (s; α, β),τi+1 ,

if ∀t ∈ Ji : F0 (t) < F̃n (τi ; α, β),for u ∈ Ji .
A key element in determining the asymptotic distribution of (α̃n , β̃n ) is the linear functional of NPMLE.

14

Lemma S2.6. Under our conditions, we haveZ√
(S2.19)n(κ(F̃n ) − κ(F0 )) = φF0 d(Qn − QF0 ) + op (n−1/2 ).

Proof. The proof of this result requires several intermediate lemmas that we present afterwards. Here we describe the crux of the arguments divided into the following main four
steps.
Step 1. The first step is to rewrite the effect from estimating the distribution using
NPMLE in terms of its linear functional:
Z√√(S2.20)n(κ(F̃n ) − κ(F0 )) = n κ̃F0 d(F̃n − F0 ).
Step 2. The second step is similar as in the proof of our Lemma S1.6 where we apply
the integration-by-parts. Now we haveZ
Z
(S2.21)
κ̃F0 d(F̃n − F0 ) = − φF̃n dQF0 ,in our Lemma S2.9, utilizing the relationship between the linear operator L and its adjointL∗ .
Step 3. Next, we need to consider the piecewise approximation of φF which has thesame jump locations as the NPMLE. By the characterization lemma S2.5, one getsZ
φ̄F̃n dQn = 0.
Thus, we haveZ

Z
(S2.22)−φF̃n dQF0 = −φ̄F̃n d(Qn − QF0 ) + op (n−1/2 ),in which we have used thatZφ̄F̃n − φF̃n dQF0 = op (n−1/2 ),as proved in our Lemma S2.10.
Step 4. In the last step, we proceed with the following decomposition
Z
Z
Z


(S2.23)− φ̄F̃n dQF0 = − φF0 d(Qn − QF0 ) +φ̄F̃n − φF0 d(Qn − QF0 ),where the second term is shown to be negligible asymptotically. Note that φF is theintegrated score function, therefore it is absolutely continuous with respect to F ; also see
our Lemma S2.7. Given that all F are (sub) distribution functions, then φF is of boundedvariation, so is its piece-wise constant approximation φ̄F . Therefore, one can show that the
random entropy integral as a function of δ is of order Op (δ 1/2 ) for the functional class that

15

includes (φ̄F̃n − φF0 ). And then by the uniform consistency of F̃n , we getZ2
φ̄F̃n − φF0 dQF0 → 0,R
with probability 1. Now,φ̄F̃n − φF0 d(Qn − QF0 ) = op (n−1/2 ) follows from the stochastic
equicontinuity of the related P-Donsker class. In the end, we arrive at
Z
Z(S2.24)− φ̄F̃n dQF0 = − φF0 d(Qn − QF0 ) + op (n−1/2 ).
Next, we list a few intermediate lemmas that are similar as the ones in Appendix B, so
instead of presenting the full proofs we only highlight the necessary changes. First, the¯following lemma states that the piece-wise constant function ξ(u)is absolutely continuous
w.r.t the NPMLE F̃n , combining the closed-form expression of φF in Van de Geer (1997)
and the proof of Lemma 4 in Geskus and Groeneboom (1996).
Lemma S2.7. The derivative of φF at the points of continuity is bounded, uniformly over
F and the points of continuity; i.e.,(S2.25)|φF (y) − φF (x)| ≤ C1 |y − x|,for y and x in the same interval between jumps and a finite positive constant C1 . Moreover,
the jumps satisfy(S2.26)|φF (x) − φF (x−)| ≤ C2 |F (x) − F (x−)|,with a finite positive constant C2 .
The next lemma controls the approximation error for a function ξ with respect to its
piece-wise version ξ¯ determined by the NPMLE. The proof is identical to the one of Lemma
S1.8.
Lemma S2.8. Suppose our Conditions hold, then we have(S2.27)k ξ¯F̃n (u) − ξF̃n (u) k2 .k F̃n (u; α0 , β0 ) − F0 (u) k2 .

Lemma S2.9. Under our conditions, We haveZ
Z
(S2.28)κ̃F0 d(F̃n − F0 ) = − φF̃n dQF0 .
Proof. The idea of the proof follows from Lemma 1 in Geskus and Groeneboom (1997)[p.251].
Let 1 ∈ L2 (F ) denote the constant function 1(u) ≡ 1 for any u ∈ [CL , CU ]. Note that for16

this constant function, L(1) = 1 after applying the L transform. Thus, we haveZ
φF dQF0 = hφF , L(1)iQF0
Z∗= hL φF , 1iF0 = L∗ (φF )dF0 .
The desired conclusion immediately follows if we can show
Z∗
(S2.29)L (φF ) = κ̃F0 − κF0 dF,
and then set F equal to F̃n . To see this, recall the solution φF0 is obtained by differentiatingL∗ φF (u) = κ̃F (u) with φF = φF0 . Then by integrating, we get
L∗ φF (u) = κ̃F0 (u) + C,
for some constant C. To pin down this constant term, we use the fact that F is non-defective
so thatZ
C = CdFZ
ZZ∗∗= L φF dF − κ̃F0 dF = hL φF , 1iF − κ̃F0 dF.
Now it is immediate that φF is contained in L02 (QF ). We also havehL∗ φF , 1iF = hφF , L(1)iF = hφF , 1iF = 0,
which completes the proof.

√The following lemma characterizes some smaller order term in our Step 3 while analyzingn(κ(F̃n ) − κ(F0 ).

Lemma S2.10. Under our conditions, we haveZ
(S2.30)φ̄F̃n − φF̃n dQF0 = op (n−1/2 ).

Proof. We start by defining the function ϕnϕn (u) = − φ̄F̃n − φF̃n (u, 0, 1)F0 (u)− φ̄F̃n − φF̃n (u, 1, 0)[F0 (u + α0 ) − F0 (u)] + φ̄F̃n − φF̃n (u, 0, 0)[1 − F0 (u + α0 )],17

and the function ξF through the relationship ξF (u) ≡ φF (u)F (u)(1 − F (u + α0 )).
Then we obtain
1 − F̃n (u)ϕn (u) =ξ¯F̃n (u) − ξF̃n (u)F̃n (u + α0 ) − F̃n (u)hi
× F0 (u + α0 )(F̃n (u) − F0 (u)) + F0 (u)(F0 (u + α0 ) − F̃n (u + α0 ))−F̃n (u + α0 )ξ¯F̃n (u + α0 ) − ξF̃n (u + α0 )F̃n (u + α0 ) − F̃n (u)hi
× (1 − F0 (u + α0 ))(F̃n (u) − F0 (u)) + (1 − F0 (u))(F0 (u + α0 ) − F̃n (u + α0 )) .

We apply the Cauchy-Schwarz inequality to getZφ̄F̃n − φF̃n dQF0 ≤ C k ξ¯F̃n − ξF̃n k2 × k F̃n − F0 k2 .
(S2.31)
Following the analogous argument as in the proof of our Lemma S1.8, we get|ξ¯F̃n (u) − ξF̃n (u)| ≤ C|F̃n (u) − F0 (u)|.
Now the result follows from the convergence result that k F̃n (u) − F0 (u) k2 = Op (log2 n ×
n−1/3 ).

We complete the section by computing the Hessian matrix related to our joint estimation
method.
Lemma S2.11. Recall that the Hessian matrix isH(α, β) ≡!
∂∂E[−X ∂αF (X 0 β; α, β)]E[−X ∂βF (X 0 β; α, β)],∂E[− ∂αF (X 0 β + α; α, β)] E[− ∂β∂ 0 F (X 0 β + α; α, β)]then we have!
E[(X − E[X|X 0 β0 ])f0 (X 0 β0 )]E[(X − E[X|X 0 β0 ])⊗2 f0 (X 0 β0 )]H0 ≡ H(α0 , β0 ) = −.
E[f0 (X 0 β0 + α0 )]E[(X − E[X|X 0 β0 ])0 f0 (X 0 β0 + α0 )]Proof. In order to avoid repetition, we only show thatE[∂
F (X 0 β + α; α, β)]|α=α0 ,β=β0 = E[(X − E[X|X 0 β0 ])0 f0 (X 0 β0 + α0 )].
0∂βFirst of all, we have0

F (u; α, β) ≡ E[∆3 |X β+α = u] =Z

F0 (u+x0 (β0 −β)+α0 −α)fX|(X 0 β+α) (x|X 0 β+α = u)dx.

Because the first slope coefficient is normalized to be 1. We denote the conditional density
function of (X2 , · · · , XK ) given X 0 β + α = u by hθ (·|u). We make the following change ofvariable by taking t1 = x0 β + α and tj = xj for j = 2, · · · , K.

18Then we can write
!
Z
K
K
X
X
F (x0 β+α; α, β) = F0 (x0 β + α −βj x̃j ) + α0 +β0j x̃j hθ (x̃2 , · · · , x̃K |x0 β+α)ΠK
˜j
j=2 dxj=2

j=2

Now we take partial derivative w.r.t. βj for j = 2, · · · , K:
(S2.32)∂F (x0 β + α; α, β)
∂βj
!
Z
K
K
X
X
= (xj − x̃j )f0 (x0 β + α −βj x̃j ) + α0 +β0j x̃j hθ (x̃2 , · · · , x̃K |x0 β + α)ΠK˜j
j=2 dx
j=2

Z
+

F0

0

(x β + α −

K
X
j=2

j=2

βj x̃j ) + α0 +

K
X
j=2

!
β0j x̃j∂hθ (x̃2 , · · · , x̃K |x0 β + α)ΠK˜j .
j=2 dx∂βj

It is evident that the first term on the right-hand side of (S2.32) is equal to E[(X −E[X|X 0 β0 ])0 f0 (X 0 β0 + α0 )]. Because the function hθ (·|u) is a conditional density function that integrates to 1, the second term on the right-hand side of (S2.32) is zero, whenevaluated at θ = θ0 . Therefore, the desired result follows.


References
Baladbaoui, F., P. Groeneboom, and K. Hendrickx (2017): “Score estimation in
the monotone single index model,” working paper.
Geskus, R., and P. Groeneboom (1996): “Asymptotically optimal estimation of smooth functionals for interval censoring, I,” Statistica Neerlandica, 50, 69–88.
(1997): “Asymptotically optimal estimation of smooth functionals for interval
censoring, II,” Statistica Neerlandica, 51, 201–219.
Groeneboom, P., and K. Hendrickx (2017): “The nonparametric bootstrap for thecurrent status model,” Electronic Journal of Statistics, 11, 3446–3483.
(2018): “Current status linear regression,” The Annals of Statistics, 46, 1415–
1444.
Groeneboom, P., and G. Jongbloed (2014): Nonparametric estimation under shape
constraints. Cambridge University Press.
Groeneboom, P., G. Jongbloed, and B. I. Witte (2010): “Maximum smoothed
likelihood estimation and smoothed maximum likelihood estimation in the current status
model,” The Annals of Statistics, 38, 352–387.

19

Groeneboom, P., and J. A. Wellner (1992): Information Bounds and NonparametricMaximum Likelihood Estimation. Birkhauser.
Van de Geer, S. (1993): “Hellinger-consistency of certain nonparametric maximum likelihood estimators,” The Annals of Statistics, 21, 14–44.
(1997): “Asymptotic normality in mixture models,” ESAIM: Probability andStatistics, 1, 17–33.
(2000): Empirical processes in M-estimation. Cambridge University Press.
Van Der Vaart, A. (1998): Asymptotic statistics. Cambridge University Press.
Van Der Vaart, A., and J. A. Wellner (1996): Weak convergence and empirical
processes. Springer.