Discussion Paper No. 552

MAXIMUM LIKELIHOOD ESTIMATION
AND SPECIFICATION TEST OF
THE PROPORTIONAL HAZARD MODELUSING GROUPED DURATIONSKeunkwan Ryu

August 2001

The Institute of Social and Economic Research
Osaka University
6-1 Mihogaoka, Ibaraki, Osaka 567-0047, Japan

Maximum Likeliho od Estimation and Sp ecification Test of
the Proportional Hazard Model Using Grouped Durations

Keunkwan Ryu*

August 1997
Revised Aug. 2001

Abstract
This article develops parametric and semi-parametric, maximum likelihood estimation
methods of the proportional hazard model for the case when durations are grouped andheterogeneity is not fully measured. This article also extends Ryu’s (1994b) specification
tests for the proportional hazard model. The proposed test statistic is easy to implement,takes into account the form of alternative hypothesis to increase power, and identifies
the source and direction of non-proportionality. This article identifies the nature of thebias resulting from neglected heterogeneity. This bias is essentially a sample selection
bias, well-known in the literature. Gamma and discrete distributions are used to model
unobserved heterogeneity. By updating the heterogeneity distribution at each stage of
duration process, this article addresses the sample selection bias. The results obtained inthis article are applied to tackle left-censoring problem.
KEY WORDS: Binary choice model; Grouped duration; Hausman’s specification test;
Invariance property; Maximum likelihood estimation; Proportional hazard model; Unobserved heterogeneity; Left-censoring.

* Keunkwan Ryu is Associate Professor, Department of Economics, Seoul National University, Seoul 151-742, KOREA, and a Visiting Foreign Scholar, ISER, Osaka University.
The author is grateful to Sunil Sharma for constructive suggestions and comments. This
research was supported by a Korea Research Foundation grant.

1. INTRODUCTION

In duration analysis, we often face a situation with continuous models and discretedata. It is natural to think of real time as continuous and that events can occur at any
moment in time. However, duration is often measured in intervals, not the exact time
elapsed, by the nature of the observation scheme. For example, many economic duration
variables constructed from longitudinal surveys are at most known only up to weeklyintervals.
The proportional hazard model (PHM) is one of the most widely used continuous-time
duration models (Cox 1972, 1975, Cox and Oakes 1984, Kalbfleisch and Prentice 1980).
Under the PHM with unobserved heterogeneity, the hazard rate is specified as a productof three separate terms: a baseline hazard function describing the overall shape of thehazard rate over time, a proportionality factor capturing the covariate (regression) eﬀects
across diﬀerent individuals, and a random variable representing unobserved heterogeneity(hereafter, heterogeneity).
Maximum likelihood estimation (MLE) of the grouped duration models has been suggested by Thompson (1977), Prentice and Gloeckler (1978), Kiefer (1988), and Sueyoshi(1991). The first paper views a grouped duration as a sequence of binary survival indicators that follow an independent Logit probability model. As discussed in Ryu (1994b),
Thompson’s parameterization is quite diﬀerent from the conventional proportional hazard
model. The other three papers consider MLE of the PHM using grouped duration data.
Additionally, Kiefer and Sueyoshi develop likelihood ratio and Lagrange multiplier tests
for the proportionality assumption.
In this article, we develop a maximum likelihood estimation method of the PHM
for the case in which durations are grouped and unmeasured heterogeneity exists. The
motivation is that for much of the available survey data, duration variables are often
interval-censored, while many of the covariates are unobserved. This article also extends
Ryu’s (1994b) specification tests for proportionality. The suggested tests are easy to use,1take alternative hypotheses into account to increase their power, and identify the sourceand direction of non-proportionality without imposing a priori restrictions.
In addition, this article suggests a flexible, and tractable parameterization of thebaseline hazard function. In the grouped duration context, the suggested parameterization
can nest the non-parametric baseline hazard function as a special case. This article also
investigates the nature of the bias resulting from neglected heterogeneity. It is essentiallya sample selection bias, well-known in the literature. Gamma and discrete distributions
are used to capture heterogeneity.
Let us briefly sketch the estimation and specification test ideas. Grouped duration canbe viewed as a sequence of binary variables indicating whether the duration survives eachinterval or not. By constructing a synthetic binary data set treating each combination(individual, interval) as a new unit of indexation, we can reduce a grouped duration
analysis to a sequential binary choice analysis.
If there exists heterogeneity, one needs to integrate the random variable representing
heterogeneity out of each interval survival probability. Here, one has to take into accountselection over time of the underlying heterogeneity. It is because the sixth year graduate
students cannot be the same as the entering graduate students in terms of their underlying
type distribution. This selection should be accounted for to avoid sample selection bias, abias due to neglected heterogeneity.
To test the proportionality assumption, we can further aggregate the already grouped
data. If proportionality holds, the two estimators, one from the further grouped data andthe other from the original grouped data, will converge to the same quantity; however,
if proportionality is violated, they will diverge from each other. Therefore, a test of theequality of these two sets of estimators will oﬀer a test for the proportionality assumption.
Finally, to illustrate the detailed aspects of the new suggested estimation and test
procedures, this article explicitly considers the following five combinations regarding (i)whether to use parametric or non-parametric baseline hazard specification and (ii) whether
and how to allow for unobserved heterogeneity: (parametric baseline, Gamma heterogene2ity), (parametric baseline, discrete heterogeneity), (non-parametric baseline, no heterogeneity), (non-parametric baseline, Gamma heterogeneity), and (non-parametric baseline,
discrete heterogeneity).
The rest of the article is organized as follows. In Section 2, after setting up a framework for discussion, the relationship between group duration analysis and sequential binary
choice analysis is shown. Also, the basics of MLE are provided in the grouped duration
context. In Section 3, heterogeneity is introduced, and operationalized within the sequential binary choice representation of the grouped duration. Proportionality tests aredeveloped in Section 4. With flexible specification of the baseline hazard function, estimation and test procedures are illustrated in detail in Section 5. In Section 6, left-censoring
will be addressed using the suggested framework. Concluding remarks follow in Section 7.
2. FRAMEWORKLet T ∈ R+ represent a duration variable of interest. Leth(t|x, v) = h0 (t) exp(x0 β)v, v ∼ g(v), v > 0,(1)

be the hazard rate of duration T , where h0 (t) is a baseline hazard function, exp(x0 β) is
a proportionality factor, x and β are k × 1 vectors of observed covariates and regression
coeﬃcients, and v captures unmeasured heterogeneity through density g(v). This model
can be termed as a PHM with unobserved heterogeneity.
To complete the model, one has to specify the baseline hazard function and the heterogeneity distribution. A simple way of treating the baseline hazard function is to makea parametric assumption. In the literature, exponential and Weibull functional forms havebeen most popular. In this article, we oﬀer an alternative flexible parameterization together with the so called non-parametric treatment. Within the current grouped duration
setting, our suggested parameterization can nest the non-parametric baseline as a special
case.
3The proportionality of the specification (1) refers to the constancy of β. That is,the covariate x increases or decreases the hazard rate by the same proportion throughout
duration. This assumption may be too strong in some situations. Ryu (1994b) examines
two empirical studies where the proportionality assumption is skeptical, and proposesa new proportionality test within a minimum χ2 estimation context. The test in thisarticle extends Ryu’s (1994b) by generalizing the observation and aggregation schemes, byintroducing heterogeneity, and by considering general covariate types.
If v degenerates to a constant, there is no unobserved heterogeneity. It is well known
in the literature that neglected or mis-specified heterogeneity leads to biased estimation.
To capture heterogeneity, Gamma and discrete distributions have been most widely used(see Heckman and Singer 1984 and Nickel 1979). These distributions are easy to use
within the proportional hazard specification in (1). In particular, the discrete distribution
can approximate any unknown distribution as the number of mass points increases withsample size.
For the identification of level in the hazard specification in (1), we need to fix level
for two of three terms in the specification, h0 (t), exp(x0 β), and v. We will leave h0 (t) free,
and fix a level for the other two terms by (i) excluding a constant term from x0 β, and by(ii) imposing one restriction on g(v) (For details, see Section 4.)
A discrete observation scheme can often be represented as an equi-spaced partition
Q of the support R+ : Q = {0, l, 2l, · · · , rl, ∞}. Under Q, the researcher keeps a record ofindividuals’ status at every l time units, until time rl elapses. Without loss of generality,
let us assume l = 1. This assumption corresponds to taking l as the measurement unit
for duration. Let Ij = [j − 1, j), j = 1, · · · , r, and Ir+1 = [r, ∞). For each observation
falling within one of r non-right-censored intervals I1 , · · · , Ir , we know its duration up to aunit interval; for each observation falling within the right-censored interval Ir+1 , we only
knows its lower bound.
Let αj (x, v) be the probability that T survives Ij conditional on that it has already4

survived all previous intervals. Then we have
Z jh(t|x, v) dt] = exp[− exp(x0 β + γj )v],αj (x, v) = exp[−

(2)

j−1

where
γj = log[

Z j

j−1

h0 (t) dt], j = 1, · · · , r.

(3)

These formulas were originally provided by Prentice and Gloeckler (1978) for the casev = 1 with probability one.
Under a parametric baseline hazard specification, the γj ’s will be functionally related.
On the other hand, under a non-parametric baseline, the γj ’s will be treated as r free parameters. Sometimes, grouped duration data involve intervals of unequal length, typicallywith wider intervals at longer durations. Our method are easily extended to this case byredefining the γj ’s in equation (3).
Let us assume that there are n independent observations. Let i index each diﬀerent
observation: i = 1, · · · , n. Define Ti as the ith duration variable, xi as the covariate ofindividual i, and dji = 1 if Ti survives Ij conditional on Ti > j − 1, dji = 0 otherwise,
j = 1, · · · , r. Then, a grouped duration can be considered as a sequence of binary indicatorvariables, that is, Ti ⇐⇒ (d1i , · · · , dri ). The eﬀective number of terms in the sequencevaries depending on at which interval the individual dies. Note that d2i ’s are meaningfully
defined only for those who have survived I1 . By constructing a synthetic binary data set
treating each combination (individual, interval) as a new unit of indexation, we can reducea grouped duration analysis to a sequential binary choice analysis (Kiefer 1988; Prenticeand Gloeckler 1978; Ryu 1994b; Sueyoshi 1991; Thompson 1977). For each combination(individual, interval), a survivor of the jth interval receives the probability
αj (x) = P (T > j|T > j − 1, x) = Ej αj (x, v)(4)

if he or she has covariate x, where Ej denotes taking expectation over v using the jth
stage heterogeneity density, say, gj (v). Since αj (x) is only defined for those who satisfy
T > j − 1, the jth stage density gj (v) should take this selection into account.
5The log-likelihood contribution of the ith individual is:
li =

r
X
j=1

sji [dji log αji + (1 − dji ) log(1 − αji )],(5)

where sji is defined by sji = d1i × · · · × dj−1i with the convention s1i = 1, and αji is a
short-hand notation for αj (xi ). This representation shows a similarity between grouped
duration analyses and sequential binary choice analyses. Passing I1 , one of two outcomesoccurs: to survive or not; Conditional on that the individual has survived I1 (conditionalon d1i = s2i = 1), again one of the same binary outcomes happens, and so forth. Here,the cumulative survival indicator sji controls whether or not the jth interval is eﬀective
for the ith individual: eﬀective if sji = 1 and not eﬀective otherwise.
By adding li over i, we derive the log-likelihood function for the whole sample:
L(θ) =

n
X
i=1

li =

n X
r
X
i=1 j=1

sji [dji log αji + (1 − dji ) log(1 − αji )],(6)

where θ denotes the collection of all model parameters such as β, and parameters in
h0 (t) and g(v). The maximum likelihood estimator of θ, say θ̂, is defined as the argument maximizing (6), or alternatively as the argument solving the first-order condition:∂L(θ̂)/∂θ = 0. By expanding the first-order condition through Taylor series and rearranging terms, we derive∂ 2 L −1 ∂L(θ)],(7)∂θ∂θ 0∂θ
where ‘=D ’ means that both hand sides of =D have the same asymptotic distribution as√
n → ∞ up to n order (for details, see Amemiya 1985).
Dθ̂ − θ = −[By taking the derivative of (5) and adding up, we obtain the score function as
n
n
r
X
∂L
∂li X Xdji − αji ∂αji
(θ) =
=
sji
.
∂θ
∂θ
αji (1 − αji ) ∂θ
i=1
i=1 j=1(8)

By exploiting (i) information matrix equality, (ii) law of large numbers, (iii) E(dji −αji )2 =αji (1 − αji ), and (iv) independence of dji across both j and i, we can derivennr∂ 2 L p X ∂li ∂li p X X
1∂αji ∂αji
=
=
−
sji
,
0
0
∂θ∂θ
∂θ ∂θ
αji (1 − αji ) ∂θ ∂θ 0
i=1
i=1 j=1

6(9)where ‘=p ’means that both sides of ‘=p ’ have the same probability limit as n → ∞ after
being divided by n.
By plugging in (8) and (9) into (7), we can represent the maximum likelihood estimator θ̂ in an asymptotically equivalent form. As is clear from the above calculation, all wehave to know is the interval survival probability αji and its derivative ∂αji /∂θ. For thispurpose, we need to complete the model by specifying the heterogeneity density g(v) andthe baseline hazard function h(t).
The next section introduces heterogeneity within the sequential binary choice framework. This oﬀers an operationally convenient way of dealing with heterogeneity. Also, thistreatment will make clear the nature of the bias resulting from neglected heterogeneity.
The bias is a sample selection bias, well-known in the literature.

3. UNOBSERVED HETEROGENEITY

We keep the sequential binary choice representation of the grouped duration. Accordingly, we are interested in updating the heterogeneity density gj (v) as time passes, and incomputing each interval survival probability αj (x) using the jth stage (interval Ij ) density
gj (v), j = 1, · · · , r.
Given a density specification for the heterogeneity v ∼ g(v), we can use this density
to integrate out the heterogeneity v from those unconditional quantities such as survival
function, distribution function, and density function. However, we cannot use the samedensity g(v) to integrate out the heterogeneity term v from those conditional quantities
such as hazard rate, and interval survival probability. It is because the information contained in the conditioning statement has an implication on the heterogeneity. For example,
those who are staying in the graduate program longer will be diﬀerent from the enteringclass in terms of their latent type distribution, that is heterogeneity. Over time, “diligent”
students will finish the program, whereas “lazy” students will still hang around, so called
weeding out eﬀect. This eﬀect reflects selection over time.
7

Defining cj = exp(x0 β + γj ), we have from (2)αj (x, v) = P (T > j|T > j − 1, x, v) = exp[− exp(xβ + γj )v] = exp(−cj v).

(10)

Due to the selection over time, the jth stage heterogeneity density gj (v) will be diﬀerent
from the initial density, g(v). The conditional density gj (v) is given by
P (T > j − 1|x, v)g(v)e−(c1 +···+cj−1 )v g(v)gj (v) = g(v|T > j − 1, x) = R ∞=,Mv (c1 + · · · + cj−1 )0 P (T > j − 1|x, v)g(v) dv(11)

where Mv (t) = Ev (e−tv ) is the moment generating function of g(v). Note that g1 (v) =g(v). Since, each cj is positive, we can easily see that gj (v) is first-order stochasticallydecreasing in j. Let vj be a random variable having density function gj (vj ). Then, vj
captures the heterogeneity in the jth interval, and vj first-order stochastically dominates
vj 0 for all j < j 0 .
Alternative way of deriving gj (v) is to rely on a pure mathematical identity. Let
S(j|x, v) and S(j|x) be the survival probability that T exceeds j, conditional on (x, v) andx, respectively. One has (i) S(j|x, v) = α1 (x, v)×· · ·×αj (x, v), (ii) S(j|x) = Ev S(j|x, v) =R∞
0 S(j|x, v)g(v) dv, and (iii) S(j|x) = E1 α1 (x, v) × · · · × Ej αj (x, v), where Ej denotes the
expectation taken with respect to gj (v). Using (i) and (ii) and arranging terms, we derive
Z ∞S(j|x) =[α1 (x, v) × · · · × αj (x, v)]g(v) dv
0
Z ∞
Z ∞
α1 (x, v) × · · · × αj−1 (x, v)g(v)=α1 (x, v)g(v) dv × · · · ×
αj (x, v) R ∞dv.
α1 (x, v) × · · · × αj−1 (x, v)g(v) dv
0
0
0
(12)
Now, by matching (iii) and (12) term by term, we can easily see thatZ ∞α1 (x, v) × · · · × αj−1 (x, v)g(v) dv
gj (v) = α1 (x, v) × · · · × αj−1 (x, v)g(v)/0Z ∞= exp[−(c1 + · · · + cj−1 )v]g(v)/exp[−(c1 + · · · + cj−1 )v]g(v) dv
0= exp[−(c1 + · · · + cj−1 )v]g(v)/Mv (c1 + · · · + cj−1 ).
8(13)

Using the updated heterogeneity density in (13), we can compute the marginal (withrespect to v) interval survival probability αj (x):Z ∞αj (x) = Ej αj (x, v) =αj (x, v)gj (v) dv
0Z ∞=exp[−(c1 + · · · + cj−1 + cj )v]g(v) dv/Mv (c1 + · · · + cj−1 )

(14)

0

= Mv (c1 + · · · + cj−1 + cj )/Mv (c1 + · · · + cj−1 )with α1 (x) = Mv (c1 ). The interval survival probability in (14) is represented as a ratio of
the moment generating functions of g(v), evaluated at two diﬀerent points. Therefore, we
readily notice that any heterogeneity density g(v) will yield an analytical solution for theinterval survival probability αj (x) insofar as g(v) admits an analytical moment generating
function.
Now, we would like to specialize the above formulas (13) and (14) for the Gamma
and discrete distributions.

Let us consider the Gamma heterogeneity first.

Assumethat v follows a Gamma distribution with parameters a and b, denoted Gamma(a, b):R∞v ∼ g(v) = ba v a−1 e−bv /Γ(a), where Γ(a) = 0 ta−1 e−t dt. It is easy to show that theb amoment generating function is Mv (t) = ( b+t) . In particular, the mean and variance are

Ev = a/b and V ar(v) = a/b2 . To fix the level of v, let us assume that Ev = 1. This
identification condition imposes a = b. That is, the heterogeneity is modeled through a
Gamma distribution with mean one and variance 1/b.
From (13),gj (v) =bbj v b−1 e−bj v
e−(c1 +···+cj−1 )v g(v)bb e−(b+c1 +···+cj−1 )v vb−1==,
Mv (c1 + · · · + cj−1 )Γ(b)Mv (c1 + · · · + cj−1 )
Γ(b)

(15)

where bj = b +c1 +· · ·+cj−1 with b1 = b. Note that gj (v) is the density corresponding to a
Gamma distribution with parameters b and bj , denoted Gamma(b, bj ). Since Ej (v) = b/bj
and V arj (v) = b/b2j , we observe that gj (v) exhibits decreasing mean and variance as jincreases. That is, adverse but homogenizing selection is going on over time in terms ofheterogeneity distribution. Now, from (14),Mv (c1 + · · · + cj−1 + cj )αj (x) ==Mv (c1 + · · · + cj−1 )9µb + c1 + · · · + cj−1
b + c1 + · · · + cj−1 + cj

¶b

.

(16)

Next, let us consider a discrete heterogeneity distribution following Heckman andSinger (1984).

Assume that v takes M finite values v1 , · · · , vM with probabilities

p1 , · · · , pM . This distribution describes that there are M unobserved types. For level
identification, let us fix vM = 1. For vm , the following parameterization is convenient:vM = 1 (level normalization), vm = ewm , m = 1, · · · , M − 1,(17)

where wm ’s are unrestricted, −∞ < wm < ∞. For pm , the following parameterization isagain useful:eπm, pm =, m = 1, · · · , M − 1,
pM =
PM−1
PM−1
1 + m=1 eπm
1 + m=1 eπm
1(18)where πm ’s are unrestricted too, −∞ < πm < ∞. The moment generating function of thisPMdiscrete distribution is Mv (t) = Ev e−tv = m=1 e−tvm pm .
Again, from (13), we derivee−(c1 +···+cj−1 )vm pm
gj (vm ) = PM
, m = 1, · · · , M.
−(c1 +···+cj−1 )vm p
m
m=1 e

(19)

Also, from (14), we obtain

Mv (c1 + · · · + cj−1 + cj )=αj (x) =Mv (c1 + · · · + cj−1 )PM−(c1 +···+cj−1 +cj )vm
pm
m=1 e
,
PM
−(c1 +···+cj−1 )vm p
m
m=1 e

j = 1, · · · , r.

(20)

We can carry out the maximum likelihood estimation of the grouped PHM with heterogeneity by maximizing the log-likelihood function (6) after plugging in αji = αj (xi ).
Under a set of quite general regularity conditions (see Amemiya 1985, ch. 4), the maximum
likelihood estimator will converge to a normal distribution with mean equal to true parameters and variance matrix equal to the inverse of the information matrix. The information
matrix can be consistently estimated using (9) by replacing the unknown parameters with
their estimates.
4. TEST FOR PROPORTIONALITY
10

The popularity of the PHM has made the issue of model checking extremely important.
This section extends Ryu’s (1994b) test for proportionality to the general grouped duration
framework with discrete/continuous covariates and observed/unobserved heterogeneity.
Ryu’s (1994b) framework is rather limited in the sense that it only considers observed
categorical covariates, neglecting unmeasured heterogeneity and continuous covariates.
By further aggregating the already grouped duration data, we can artificially generate
another coarser set of grouped duration data. The aggregation can be represented asa contiguous grouping of integers, {1, · · · , r}. Let Q∗ = {G1 , · · · , Gr∗ }, where Gk ={g1 + · · · + gk−1 + 1, · · · , g1 + · · · + gk−1 + gk } with gk being the number of elements inthe subset Gk , gk ≥ 1, k = 1, · · · , r∗ . Then, Q∗ corresponds to the following aggregation
scheme: aggregate the first g1 intervals into a big interval, say I1∗ ; aggregate the next g2

intervals into a big interval, say I2∗ ; continue in the same way until we obtain the last big
interval Ir∗∗ . At least one of the gk ’s must be greater than one. Otherwise, no aggregation
occurs. Under Q∗ , the duration data will take the form d∗ki = Πj∈Gk dji . Here, d∗kitakes value one if individual i survives the big interval Ik∗ = ∪j∈Gk Ij and zero otherwise.

Moreover, α∗ki = Πj∈Gk αji = exp[− exp(xβ + γk∗ )] is the survival probability of the big
RPinterval Ik∗ , where γk∗ = log t∈I ∗ h0 (t) dt = log( j∈Gk eγj ), k = 1, · · · , r∗ . Obviously, theknew coarser data set Q∗ contains less information than the original finer data set Q. Let β̂be the estimate of β obtained by using the original data Q, and β̂ ∗ the estimate obtainedby using the further aggregated data Q∗ . By comparing these two estimates, we can design
a new proportionality test statistic.
If the PHM holds, both β̂ and β̂ ∗ will converge to the same β. However, if the PHM
fails to hold, β̂ and β̂ ∗ will converge to diﬀerent quantities. To see this, let us consider asimple case where r = 2, r∗ = 1, Q = {0, 1, 2, ∞}, Q∗ = {G1 }, and G1 = {1, 2}. Assume
that all covariates have decreasing impacts on the hazard rates:h(t, x) =½

h0 (t) exp(xβ 1 ),h0 (t) exp(xβ 2 ),for t ∈ I1 ;for t ∈ I2 ,with β 1 > β 2 > 0. Here, heterogeneity is assumed away by taking v = 1 with probability
11one. The aggregate estimator β̂ ∗ , obtained from the coarser data I1∗ = I1 ∪ I2 , is symmetrically aﬀected by β 1 and β 2 since I1∗ is the symmetric aggregation of I1 and I2 . On theother hand, the disaggregate estimator β̂, obtained from the finer data set I1 and I2 , isasymmetrically aﬀected by β 1 and β 2 since there are more observations in I1 than in I2(note that some individuals have died in I1 ). Therefore, the influence of β 1 relative to β 2
is stronger on β̂ than on β̂ ∗ . Since β 1 is larger than β 2 , β̂ should be stochastically larger
than β̂ ∗ . As a result, the diﬀerence between β̂ and β̂ ∗ converges to a zero vector under
the PHM, but to a non-zero vector under non-proportionality. The proposed test statistic
uses this disparate convergence pattern: if the diﬀerence β̂ ∗ − β̂ is significantly diﬀerent
from zero, reject the PHM; otherwise, do not.
The test statistic
R = (β̂ − β̂ ∗ )0 [V ar(β̂ − β̂ ∗ )]−1 (β̂ − β̂ ∗ )(21)follows a chi-square distribution with degrees of freedom equal to the number of parametersin β, say k. The variance inside the bracket can be expanded as
V ar(β̂ − β̂ ∗ ) = V ar(β̂) + V ar(β̂ ∗ ) − Cov(β̂, β̂ ∗ ) − Cov(β̂, β̂ ∗ )0 .

(22)

In calculating the above test statistic, the diﬃculty usually lies in computing the covariancematrix between the two estimators. It is so because variances can be easily estimatedthrough the inverse of the observed information matrices.
Then, how to compute the covariance matrix? It depends on whether we adopt aparametric baseline specification or not. Obviously, the original finer data set contains
bigger amount of information than the aggregated coarser data set. However, this ranking
in data information contents does not necessarily yield an eﬃciency ranking between β̂
and β̂ ∗ . If we make parametric baseline hazard assumption and thus estimate the same
fixed number of baseline hazard parameters in both data set-ups, we obtain a result thatβ̂ is more eﬃcient than β̂ ∗ (Ryu 1993a) and that their covariance reduces to V ar(β̂)(Hausman 1978). Obviously, the number of free parameters in the parametric specification
12

should not exceed the number of non-right-censored intervals r∗ in the coarser data set(see sub-section 5.5.). On the other hand, if we do not make a parametric functional
form assumption on the baseline hazard, the number of baseline hazard parameters being
estimated under the two data set-ups are diﬀerent: when we are estimating the PHM using
the finer data set, we are estimating more steps regarding the baseline hazard function(r > r∗ ). In this non-parametric case, the ranking in data information contents does not
imply an eﬃciency ranking between β̂ and β̂ ∗ . This lack of eﬃciency ordering prohibits
one from using Hausman’s (1978) result to simplify the covariance between β̂ ∗ and β̂ inthe non-parametric baseline case. In the next section, we will explain how to compute thecovariance matrix in the non-parametric baseline case.
Besides the overall chi-square test, we can also conduct individual t-tests. Under the
proportionality assumption,
tj = (β̂j − β̂j∗ )/σ(β̂−β̂ ∗ )j , j = 1, · · · , k(23)will have an asymptotic standard normal distribution, where β̂j and β̂j∗ are the jth elementsof β̂ and β̂ ∗ , and σ(β̂−β̂ ∗ )j is the square root of the jth diagonal element of V ar(β̂ − β̂ ∗ ).
The advantage of individual t-tests is to separately identify those covariates which exhibit
non-proportional eﬀects, and to inform the direction of those non-proportional eﬀects.
For instance, if tj is significantly positive (negative), then we can conclude that the jthcomponent of x has a non-proportional eﬀect on the hazard rate and that its coeﬃcient islarger (smaller) in the early intervals than in the later intervals. If the individual t-tests
detect non-proportionality, then it would be a better idea to estimate a non-proportional
hazard model by allowing diﬀerent β’s for those non-proportional covariates across diﬀerent
intervals.
If we have more than two non-right-censored intervals in the original data set (r > 2),
we have a lot more flexibility in choosing a further aggregation. For testing purposes, theselection of an optimal aggregation may be guided by whatever alternative hypothesis one
has in mind. For example, if there are three non-right-censored intervals in the original
13

data set (r = 3) and if one suspects that the covariate impact is weaker in the third interval(if there is any diﬀerence), one may put the first two intervals together, leaving the thirdinterval alone. That is, take r∗ = 2, Q∗ = {G1 , G2 }, G1 = {1, 2}, G2 = {3}, I1∗ = I1 ∪ I2 ,and I2∗ = I3 . This aggregation will yield a higher power against the suspected alternative
than any other aggregation. Of course, without a clear alternative in mind, we cannot
design an optimal aggregation.
The suggested test is very easy to implement, and allows one to take into account theform of suspected alternative hypothesis to increase power of the test. Most of all, the test
is much more convenient to use compared with other existing tests. To apply the likelihood
ratio test, one has to estimate the model under a non-proportional alternative. If there
are ten covariates (k = 10) and ten non-right-censored intervals (r = 10) in the original
data, the most general form of non-proportional hazard model will include 100(= 10 × 10)
β’s, too many parameters! To reduce the number of parameters, one has to introduce
a priori restrictions by assuming either that some covariates have proportional eﬀects orthat the non-proportional eﬀects satisfy certain parametric restrictions. The situation
is practically no better for the Lagrange multiplier test. Even though we only need toestimate the model under proportionality, we have to consider a very long vector of scorefunctions corresponding to 100 β’s. Again, we are forced to adopt the same a priori
assumptions as in the likelihood ratio test to solve the dimensionality problem. However,
our test identifies the source and nature of non-proportionality without imposing any apriori restrictions, a useful property not shared by the existing tests.

5. BASELINE HAZARD FUNCTION

This section specifies the baseline hazard function, and details the estimation andtest procedures. Both parametric and non-parametric specifications are introduced andcompared with each other. We first introduces the general framework in sub-section 5.1.
In sub-sections 5.2 through 5.4, we specialize the general framework to the cases of (non14

parametric baseline, no unmeasured heterogeneity), (non-parametric baseline, Gamma
heterogeneity), and (non-parametric baseline, discrete heterogeneity). Finally, parametricspecifications are studied in sub-section 5.5.
5.1. A general framework
In the nonparametric baseline case, γj ’s (γk∗ ’s) can be treated as r (r∗ ) free parametersin the finer (coarser, respectively) data set. It is because h0 (t) is left unspecified. On
the other hand, in the parametric case, γj ’s (γk∗ ’s) are functionally related due to theparametric restriction imposed on h0 (t). To be able to use the same parametric baseline
hazard specification across the finer and the coarser data sets, we should restrict the
number of free parameters in the parametric specification not to exceed r∗ , the number of
non-right-censored intervals in the coarser data set.
Let θ denote the vector of all model parameters in the finer data set, which includeβ, parameters in unmeasured heterogeneity if any, and parameters in the baseline hazard
function. Let θ∗ denote the vector of all parameters in the model applicable to the coarser
data set. In the non-parametric baseline case, the number of parameters are reduced
due to aggregation in the coarser data set, resulting in dim(θ∗ ) < dim(θ). But, in theparametric baseline case, we have θ = θ ∗ .
Under general regularity conditions (see Amemiya 1985, ch. 4), we have from (7)-(9)in Section 2.
nr1∂αji ∂αji −1 X X
dji − αji ∂αji
sji
[
sji
]
θ̂ − θ = [0 ]αji (1 − αji ) ∂θ ∂θ
αji (1 − αji ) ∂θ
i=1 j=1i=1 j=1
D

n X
r
X

n X
r
X

1
∂αji ∂αji −1∼ N (0, [
sji
).
0 ]
α
(1
−
α
)
∂θ
∂θ
ji
ji
i=1 j=1(24)

where the variance-covariance matrix can be evaluated at the estimated parameter values
for practical use.
Next, let us consider the coarser data set Q∗ = {G1 , · · · , Gr∗ } introduced in Section

4. Note that d∗ki = Πj∈Gk dji denotes the outcome whether the ith individual survives
15

Ik∗ or not. Accordingly, α∗ki = Πj∈Gk αji denotes the probability of surviving the interval
Ik∗ . Now θ∗ is of a smaller dimension than θ due to aggregation, with the diﬀerence beingr − r∗ . By a similar procedure, we obtain
∗

n X
r
X
∗
θ̂ − θ = [s∗kiD

i=1 k=1

∗

n X
r
∗
∗
X
1
∂αki
∂αki
d∗ki − α∗ki ∂α∗ki
−1
∗
]
[
s
ki ∗
∗ ) ∂θ ∗ ∂θ ∗0∗ ) ∂θ ∗ ]α∗ki (1 − αki
α
(1
−
α
ki
ki
i=1
k=1

n X
r∗
X

(25)

∗
1
∂α∗ki ∂αki
−1s∗ki ∗).
∼ N (0, [0 ]∗∗∗αki (1 − αki ) ∂θ ∂θ
i=1 k=1

Finally, let us compute the covariance between θ̂ and θ̂ ∗ . For this purpose, let us
rewrite the original estimator θ̂ as follows:
∗

n X
r
X

n r
1
∂αji ∂αji −1 X X X
dji − αji ∂αji=].
θ̂ − θ [
sji
]
sji
0
αji (1 − αji ) ∂θ ∂θ
αji (1 − αji ) ∂θ
i=1 j=1
i=1
D

(26)

k=1 j∈Gk

Now note that dji and d∗ki0 are independent either i 6= i0 or j ∈/ Gk . This is due tothe independence of survival indicators either for diﬀerent individuals or across nonoverlapping intervals. So, let us compute the covariance between dji and d∗ki0 for thecase i = i0 and j ∈ Gk . Observing that dji d∗ki = d∗ki for all j ∈ Gk , we can easily obtaincov(dji , d∗ki ) = α∗ki (1 − αji ).
Therefore,∗cov(θ̂, θ̂ ) = [n XrX
i=1 j=1

n X
r∗
X

×[

i=1 k=1sji

1∂αji ∂αji −1]
αji (1 − αji ) ∂θ ∂θ0∗n r∗
1 ∂αji1∂α∗ki X X ∗1∂α∗ki ∂αki
(
sji
)
][
s
]−1 .
ki ∗
∗
∗
∗0
∗ ∂θ ∗0
α
∂θ
1
−
α
∂θ
α
(1
−
α
)
∂θ
ji
ki
ki
ki
i=1
j∈G

X

(27)

k=1

k

Using (24), (25), and (27), we can obtain the variance matrices of β̂, β̂ ∗ , and theircovariance matrix, by figuring out the relevant blocks corresponding to β’s out of θ’s.
Then, these results can be used to compute the specification test statistics R and tj ’s in
Section 4.
5.2. Non-parametric baseline, no unmeasured heterogeneity16

This sub-section analyzes the case of non-parametric baseline hazard function without
unmeasured heterogeneity. Without imposing a parametric assumption on the baseline
hazard function h0 (t), we can still estimate the integrated baseline hazard over each of r
non-right-censored intervals I1 , · · · , Ir . In the current setting, we are able to consistently
R1
Rr
estimate γ1 , · · · , γr or, equivalently, exp(γ1 ) = 0 h0 (t) dt, · · · , exp(γr ) = r−1 h0 (t) dt.

This means that we can approximate the unknown function h0 (t) up to a step function
with r diﬀerent steps, insofar as r is either finite or increasing at a rate slower than
the sample size n. The model parameters can be represented as a (k + r) × 1 vector,θ = (β 0 , γ1 , · · · , γr )0 .
Note that0∂αji /∂θ = −αji exi β+γj zji ,(28)

where zji = (x0i , e0j )0 is a (k + r) × 1 column vector and ej = (0, · · · , 0, 1, 0, · · · , 0)0 is an
r × 1 unit column vector with one occupying the jth location. Thus,
n X
r
r
X
X
∂L
dji − αji x0i β+γj(θ) = −sjiezji = −
Zj0 Sj Ω−1
j uj ,∂θ1 − αji
i=1 j=1j=1(29)

where Zj = (zj1 , · · · , zjn )0 is an n × (k + r) matrix, Ωj is an n × n diagonal matrix of0(1 − αji )/(αji (exi β+γj )2 ), Sj is an n × n diagonal matrix of sji , and uj is an n × 1 vector0of (dji − αji )/(αji exi β+γj ). Using V ar(uj ) = Ωj , j = 1, · · · , r, the information matrix canbe approximated as
r
r
r
X
p X 0
∂ 2 L p ∂L ∂L p X 0
−1
−1
0
0
=
=[
−
Zj Sj Ωj uj ][Zj Sj Ωj uj ] =Zj Sj Ω−1
j Zj .
∂θ∂θ0∂θ ∂θ 0
j=1
j=1
j=1

(30)

Using these results, we haveDθ̂ − θ = −[r
X

−1
Zj0 Sj Ω−1
j Zj ]

j=1∼ N (0, [r
X

Zj0 Sj Ω−1
j uj

j=1

r
X

−1
Zj0 Sj Ω−1).
j Zj ]

j=1

17

(31)

Similarly for the coarser data set Q∗ = {G1 , · · · , Gr∗ }, we obtainD∗θ̂ ∗ − θ = −[

r
X

∗

0

Zk∗ Sk∗ Ω∗−1
Zk∗ ]−1
k

k=1

r
X

k=1

∗

0

Zk∗ Sk∗ Ω∗−1u∗k ∼ N (0, [
k

r
X

0

Zk∗ Sk∗ Ω∗−1
Zk∗ ]−1 ), (32)
k

k=1
0

∗
∗ 0
∗
where Zk∗ = (zk1
, · · · , zkn) is an n×(k +r∗ ) matrix with zki= (x0i , e∗k )0 being a (k +r∗ )×1

vector (e∗k is an r∗ × 1 unit vector with one in the kth location), Ω∗k is an n × n diagonal
0

∗

∗
∗
matrix of (1 − αki)/(αki(exi β+γk )2 ), Sk∗ is an n × n diagonal matrix of s∗ki , and u∗k is an0∗∗
n × 1 vector of (d∗ki − αki)/(α∗ki exi β+γk ). Again, note that S1∗ is an identity matrix of ordern, and that V ar(u∗k ) = Ω∗k , k = 1, · · · , r∗ .

Finally, the covariance between θ̂ and θ̂ ∗ iscov(θ̂, θ̂∗ )
r∗
r∗
X
¡ X
¢ ∗−1 ∗ X
0
−1
0
−1
0
∗ −1
=[
Zj Sj Ωj Zj ] [fj Zj Sj Ωk Zk ][Zk∗ Sk∗ Ω∗−1,
k Zk ]
j=1
k=1 j∈Gk
k=1
r
X

(33)

∗−1
∗
where Sj Ω∗−1(for all j ∈ Gk ) has been used, a property resulting fromk Sk = Sj Ωksji s∗ki = sji for all j ∈ Gk .

5.3. Non-parametric baseline, Gamma heterogeneity
This sub-section analyzes the case of non-parametric baseline with unmeasured heterogeneity modeled according to a Gamma distribution.

Assume that vi follows a

Gamma distribution with parameters a and b, denoted Gamma(a, b): vi ∼ i.i.d. g(v) =R∞ba v a−1 e−bv /Γ(a), where Γ(a) = 0 ta−1 e−t dt. It is easy to show that the moment generat-b a
ing function is Mv (t) = ( b+t) . The mean and variance are Ev = a/b and V ar(v) = a/b2 .

To fix the level of v, let us assume that Ev = 1, rendering a = b. That is, the unmeasured
heterogeneity is modeled through a Gamma distribution with mean one and variance 1/b.
The model parameters can be represented as a (k + 1 + r) × 1 vector, θ = (β 0 : b :b bγ1 , · · · , γr )0 . Using the moment generating function, Mv (t) = ( b+t) , we obtain

αji = (

b + c1i + · · · + cj−1i)b , j = 1, · · · , r,b + c1i + · · · + cj−1i + cji
18

(34)

where cji = exp(x0i β + γj ). Now, it is straightforward to compute ∂αji /∂θ = (∂αji /∂β 0 :∂αji /∂b : ∂αji /∂γ1 , · · · , ∂αji /∂γr )0 . We have
b−1
∂αji
cji
= −b2 [αji ] b
xi ,∂β(b + c1i + · · · + cj−1i + cji )2
b−1
∂αji
αji
cji
,
=
log αji + b[αji ] b∂bb(b + c1i + · · · + cj−1i + cji )2
b−1
∂αji
1
= b[αji ] b∂γk(b + c1i + · · · + cj−1i + cji )2(35)

× [cji (c1i δ1k + · · · + cj−1i δj−1k ) − cji δjk (b + c1i + · · · + cj−1i )],where δjk takes value one if j = k, and zero otherwise.
Under the coarser grouping, the model parameters becomes a (k + 1 + r∗ ) × 1 vector,θ ∗ = (β 0 : b : γ1∗ , · · · , γr∗∗ )0 . Similarly, we have
∗
αki
=(

b + c∗1i + · · · + c∗k−1i)a , k = 1, · · · , r∗ ,b + c∗1i + · · · + c∗k−1i + c∗ki

(36)where c∗ki = exp(x0i β + γk∗ ). From this, we obtain
b−1
∂α∗ki
c∗ki
= −b2 [α∗ki ] b
xi ,∂β(b + c∗1i + · · · + c∗k−1i + c∗ki )2
α∗
c∗ki
∂α∗ki
∗ b−1
= ki log α∗ki + b[αki
] b
,
∂b
b
(b + c∗1i + · · · + c∗k−1i + c∗ki )2
b−1
∂α∗ki
1
= b[α∗ki ] b∗∂γj(b + c1i + · · · + c∗k−1i + cki )2(37)

× [c∗ki (c∗1i δ1j + · · · + c∗k−1i δk−1j ) − c∗ki δkj (b + c∗1i + · · · + c∗k−1i )],where δkj takes value one if k = j and zero otherwise.
By combining the above results with the general framework introduced in sub-section
5.1, we can make inferences on θ and θ∗ , and carry out the proportionality test.
5.4. Non-parametric baseline, discrete heterogeneity
This sub-section analyzes the case of non-parametric baseline hazard functionwith
unmeasured heterogeneity modeled according to a discrete distribution. Assume that vi
takes M finite values v1 , · · · , vM with probabilities p1 , · · · , pM . This distribution describes
19that there are M unobserved types. For level identification, let us fix vM = 1. For vm ,the following parameterization is convenient:vM = 1 (level normalization), vm = ewm , m = 1, · · · , M − 1,(38)

where wm ’s are unrestricted, −∞ < wm < ∞. For pm , the following parameterization isuseful:
pM =

1

1+

PM−1 πmm=1 e, pm =eπm, m = 1, · · · , M − 1,
PM−1
1 + m=1 eπm(39)where πm ’s are unrestricted too, −∞ < πm < ∞.

The model parameters can be represented as a (k + 2M − 2 + r) × 1 vector, θ =(β 0 : w1 , · · · , wM−1 : π1 , · · · , πM−1 : γ1 , · · · , γr )0 . Using the moment generating function,PM
Mv (t) = Ev e−tv = m=1 e−tvm pm , we obtainPM−(c1i +···+cj−1i +cji )vmpm
Mv (c1i + · · · + cj−1i + cji )
m=1 e
= PM
αji =
, j = 1, · · · , r. (40)+···+c
)v
−(c
1i
j−1i
mp
Mv (c1i + · · · + cj−1i )m
m=1 e
From this, we can compute ∂αji /∂θ = (∂αji /∂β 0 : ∂αji /∂w1 , · · · , ∂αji /∂wM−1 :∂αji /∂π1 , · · · , ∂αji /∂πM−1 : ∂αji /∂γ1 , · · · , ∂αji /∂γr )0 . For notational convenience, letPM
us define aji = c1i + · · · + cji and δji = m=1 e−aji vm pm , we have
M
∂αji
1 X
=
pm vm (αji aj−1i e−aj−1i vm − aji e−aji vm ),∂βδji m=1∂αji
1=
pk vk (αji aj−1i e−aj−1i vk − aji e−aji vk ),
∂wk
δji

M
X
1
1
∂αji
=
pk (e−aji vk − αji e−aj−1i vk ) +pkpm (αji e−aj−1i vm − e−aji vm )∂πkδjiδji m=1

(41)

M
∂αji
cki X
=
pm vm [αji e−aj−1i vm 1(k≤j−1) − e−aji vm 1(k≤j) ],∂γkδji m=1

where a0i = 0.

Under the coarser grouping, the model parameters becomes a (k + 2M − 2 + r∗ ) × 1

vector, θ∗ = (β 0 : w1 , · · · , wM −1 : π1 , · · · , πM −1 : γ1∗ , · · · , γr∗∗ )0 . Similarly, we have
PM
∗
∗
−(c∗
1i +···+ck−1i +cki )vm p
m
∗
m=1 e
αki = PM, k = 1, · · · , r∗ , k = 1, · · · , r∗ ,+···+c∗
)vm
−(c∗
1i
k−1i
pm
m=1 e
20

(42)

where c∗ki = exp(x0i β + γk∗ ).

∗From this, we can compute ∂αji/∂θ∗ = (∂α∗ji /∂β 0 :∗/∂w1 , · · · , ∂α∗ji /∂wM−1 : ∂α∗ji /∂π1 , · · · , ∂α∗ji /∂πM −1 : ∂α∗ji /∂γ1∗ , · · · , ∂α∗ji /∂γr∗ )0 . For
∂αji
PM
∗
∗
notational convenience, let us define a∗ji = c∗1i + · · · + c∗ji and δji= m=1 e−aji vm pm . We

have

M
∂α∗ji
∗
∗
1 X
= ∗
pm vm (α∗ji a∗j−1i e−aj−1i vm − a∗ji e−aji vm ),∂βδji m=1

∂α∗ji
∗
∗
1
= ∗ pk vk (α∗ji a∗j−1i e−aj−1i vk − a∗ji e−aji vk ),
∂wk
δji

M
X
∂α∗ji
∗
∗
1
1
vk
vk
−a∗
∗ −a∗
ji
j−1i
= ∗ pk (e− αji e) + ∗ pkpm (α∗ji e−aj−1i vm − e−aji vm )∂πkδjiδji m=1

(43)

M
∂α∗ji
∗
∗
c∗ki X=pm vm [α∗ji e−aj−1i vm 1(k≤j−1) − e−aji vm 1(k≤j) ],∗∗∂γkδji m=1

where a∗0i = 0.
Again, by combining the above results with the general framework introduced in
sub-section 5.1, we can make inferences on θ and θ ∗ , and carry out the proportionality
test.
5.5. Parametric baseline, Gamma and discrete heterogeneity
Once we parameterize h0 (t) using m free parameters, we have a fully parametric
duration model. This article proposes one to use the following flexible specification.
γj = log

Z jj−1

h0 (t) dt = δ0 + δ1 j + · · · + δm−1 j m−1 .

(44)This parameterization is polynomial in time, and it guarantees positivity of the integrated
baseline hazard over each interval, Ij . The exponential distribution (constant baseline
hazard) is a special case when m = 1. The monotonic hazard feature of the Weibull
distribution is well captured by choosing m = 2. As m increases, the resulting parametricspecification becomes more flexible.
21

As a practical guide, if the baseline hazard function is expected to be monotone, use
m = 2. If U or inverted U-shaped, use m = 3. Without any clear idea on the shape, applymodel selection criteria to the choice of m. If data are rich enough, use the non-parametricbaseline specification.
If the number m of free parameters is equal to the number r of non-right-censoredintervals, then both the parametric and non-parametric baseline specifications would yield
the same results. This is because of the invariance property of the maximum likelihood
estimation. If m is greater than r, then we are not able to identify all these m parametersfrom grouped duration data with just r non-right-censored intervals. We can only identify
r restrictions on m parameters, because the baseline hazard parameters enter only the
likelihood function through γ1 , · · · , γr . On the other hand, if m is smaller than r, then weare virtually imposing r − m parametric restrictions on those r integrated baseline hazard
rates. For identification of the parametric baseline hazard function (44), we readily notethat the number of free parameters in h0 (t) should not exceed the number of non-rightcensored intervals in the data set, m ≤ r.
In the parametric baseline cases, the maximum likelihood estimation and specification
tests are straightforward: (1) obtain two sets of parametric maximum likelihood estimators, one from the original data set and the other from the new coarser data set, (2)estimate each variance as the inverse of the observed information matrix, and (3) estimate
their covariance as the variance of the estimator from the original data set, a result due
to Hausman (1978).
From the results in the non-parametric baseline case, we can derive the partial deriva∗
with respect to each model parameters. The partial derivatives withtives of αji and αki

respect to β and parameters in unmeasured heterogeneity are the same as before. Onlychanges occur in the partial derivatives of αji and α∗ki with respect to the parameters inthe baseline hazard function. But, this can be handled easily using the chain rule of differentiation. For example, ∂αji /∂δk = (∂γ1 /∂δk , · · · , ∂γr /∂δk )(∂αji /∂γ1 , · · · , ∂αji /∂γr )0 .
∗Similarly for ∂αji/∂δk .

22

To allow maximum flexibility within the parametric baseline specification (44), one
can use m = r∗ , resulting in
γj = log

Z j

j−1

∗

h0 (t) dt = δ0 + δ1 j + · · · + δr∗ −1 j r −1 , j = 1, · · · , r.

(45)

Under this parameterization, the coarser estimator β̂ ∗ is in fact equivalent to thesemi-parametric estimator, whereas the original estimator β̂ is an estimator obtained after
imposing r − r∗ restrictions on the baseline hazard function.
6. APPLICATION TO LEFT-CENSORING
The suggested framework can be used to deal with left-censoring issue in duration
analysis. If every individual is observed from the start of his or her episode (called, flow
sampling), there is no problem of left-censoring. Often, however, individuals are observed
to be already in the middle of an episode (called, stock sampling). The resulting duration
variable is said to be left-censored, and is known to complicate the estimation (see, forexample, Amemiya 1985, ch. 11).
Considering that a longer duration is more likely to be in progress at a random start
time of observation, the fact that an episode is left-censored implies that the corresponding
duration is more likely to be longer than a typical duration (described as length biased
sampling). We would like to address this problem by assuming that there is a latent
unobserved variable aﬀecting duration, say v. The essential feature of left-censoring canbe described by the changing distribution of v over the duration process. Therefore, by
updating the distribution of v at each stage of the process as suggested in this article, we
can account for the selectivity, thus for the left-censoring.
To be concrete, for left-censored observations, we replace αj (x) in eq. (14) with thefollowing modified interval survival probability:
αj (x|s) = Es+j αj (x, v) =

Z ∞
0

23

αj (x, v)gs+j (v) dv,(46)

where gs+j (v) is the s + jth stage heterogeneity density, gs+j (v) = g(v|T > s + j − 1, x)(see eq. (11)). Note that to be able to compute αj (x|s), we need information on how longthe process has already been in progress, that is, elapsed duration s, by the time of firstobservation.
7. CONCLUDING REMARKS
Often, duration data are available in a grouped form due to a certain discrete observation mechanism, while many covariates are unobservable. This article develops a generalsequential binary choice framework of grouped duration model with unobserved heterogeneity. Considering that many economic duration data are grouped and many variables
are missing, the proposed methods will prove useful in many situations.
This article introduces a new flexible parameterization of the baseline hazard function. Non-parametric baseline hazard function is covered as a special case of this flexible
parameterization. This article proposes a new, and operationally convenient, way of tackling unobserved heterogeneity from a sample selection perspective. Gamma and discrete
distributions have been used to capture heterogeneity. By updating the heterogeneity
distribution, we can account for selection of type over duration, and thus can avoid theneglected heterogeneity bias. In fact, the bias resulting from neglected heterogeneity isessentially a sample selection bias. This interpretation allows us to keep the sequential
binary choice representation of the grouped duration model, and still enables us to estimate the model easily. Also, this article extends Ryu’s (1994b) proportionality test
to a general grouped duration setting. The test can detect the source and direction of
non-proportionality without adopting a priori restrictions, not shared by existing test procedures.
The general framework suggested in this article proves useful in addressing leftcensoring problem.

24

REFERENCES
Amemiya, T. (1985), Advanced Econometrics, Boston, MA: Harvard University Press.
Amemiya, T. (1973), “Regression Analysis when the Dependent Variable isTruncated Normal,” Econometrica, 41, 997—1017.
Berndt, E. K., Hall, B. H., Hall, R. E., and Hausman, J. A. (1974), “Estimation andInference in Nonlinear Structural Models,” Annals of Economic and Social
Measurement, 3, 653—665.
Cox, D. R. (1972), “Regression Models and Life-Tables,” Journal of the Royal StatisticalSociety B, 34, 187—220.
Cox, D. R. (1975), “Partial Likelihood,” Biometrika, 62, 269-279.
Cox, D. R., and Oakes, D. (1984), Analysis of Survival Data, London: Chapman & Hall.
Han, A. K., and Hausman, J. A. (1990), “Flexible Parametric Estimation of Duration
and Competing Risk Models,” Journal of Applied Econometrics, 5, 1—28.
Hausman, J. A. (1978), “Specification Tests in Econometrics,” Econometrica, 46,
1251-1272.
Hausman, J. A., and Wise (1979), Econometrica.
Heckman, J., 1979, “ Sample Selection Bias as a Specification Error,” Econometrica,
47, 153—161.
Heckman, J., 1976, “The Common Structure of Statistical Models of Truncation,
Sample Selection and Limited Dependent Variables and a Simple Estimator for
Such Models,” The Annals of Economic and Social Measurement, 5, 475—492.
Heckman, J., and Singer, B. (1984), “A Method for Minimizing the Impact ofDistributional Assumptions in Econometric Models for Duration Data,”
Econometrica, 52, 271-320.
Jennrich, R., 1969, “Asymptotic Properties of Non-Linear Least Squares Estimators,”
The Annals of Mathematical Statistics, 40, 633—643.
Kalbfleisch, J. D., and Prentice, R. L. (1980), The Statistical Analysis of Failure Time25

Data, New York, NY: John Wiley & Sons.
Kiefer, N. M. (1988), “Analysis of Grouped Duration Data,” Contemporary Mathematics,
80, 107—137.
Lee, L. F., 1982, “Some Approaches to the Correction of Selectivity Bias,” Review of
Economic Studies, 49, 355—372.
Lin, D. Y. (1991), “Goodness-of-Fit Analysis for the Cox Regression Model Based on aClass of Parameter Estimators,” Journal of the American Statistical Association,
86, 725—728.
Newey, W., and D. McFadden (1994), “Large Sample Estimation and Hypothesis Testing,”
in R. F. Engle and D. McFadden, eds., Handbook of Econometrics, IV, Elsevier.
Moﬃtt, R. (1985), “Unemployment Insurance and the Distribution of Unemployment
Spells,” Journal of Econometrics, 28, 85—101.
Newey, W. (1984), “A Method of Moments Interpretation of Sequential Estimators,”
Economics Letters, 14, 201—206.
Newey, W. and K. West, 1987, “A Simple Positive Semi-Definite, Heteroscedasticity
and Auto-correlation Consistent Covariance Matrix,” Econometrica, 55, 703—708.
Nickell, S. (1979), “Estimating the Probability of Leaving Unemployment,” Econometrica,
47, 1249-1266.
A. Pagan, 1986, “Two Stage and Related Estimators and Their Applications,” Review
of Economic Studies, 53, 517—538.
A. Pagan, 1984, “Econometric Issues in the Analysis of Regressions with Generated
Regressors,” International Economic Review, 25, 221—247.
Prentice, R., and Gloeckler, L. (1978), “Regression Analysis of Grouped Survival Datawith Application to Breast Cancer Data,” Biometrics, 34, 57—67.
Ryu, Keunkwan (1994a), “Modified Heckman’s Two-step Estimation and Specification
Test in Sample Selection Models,” mimeo., UCLA.
Ryu, K. (1994b), “Group Duration Analysis of the Proportional Hazard Model:
Minimum Chi-square Estimators and Specification Tests,” Journal of the26

American Statistical Association, 89, 1386—1397.
Ryu, Keunkwan (1993a), “Monotonicity of the Fisher Information and the KullbackLeibler Divergence Measure,” Economics Letters, 42, 121-128.
Ryu, K. (1993b), “An Extension of Marshall and Olkin’s Bivariate Exponential
Distribution,” Journal of the American Statistical Association, 88, 1458—1465.
Sueyoshi, G. T. (1991), “Evaluating Simple Alternatives to the Proportional Hazard
Model: Unemployment Insurance Receipt and the Duration of Unemployment,”
mimeo., San Diego, CA: University of California.
Thompson, W. A. (1977), “On the Treatment of Grouped Observations in Life Studies,”
Biometrics, 33, 463—470.
White, H., 1980, “A Heteroscedasticity-Consistent Covariance Matrix Estimator and A
Direct Test for Heteroscedasticity,” Econometrica, 48, 817—838.

27