Geert Dhaene
email: geert.dhaene@econ.kuleuven.be.

K U Leuven

Koen Jochmans
email: koen.jochmans@uclouvain.be.

AN ADJUSTED PROFILE LIKELIHOOD FOR NON-STATIONARY PANEL DATA MODELS WITH FIXED EFFECTS

Keywords: JEL classification: C13, C22, C23 adjusted likelihood, bias correction, dynamic panel data, fixed effects, non-stationarity

), with fixed effects integrated out using a new data-independent prior. It appears that ρ and β are identified as the unique point where the large N adjusted profile likelihood reaches a local maximum (or a flat inflection point, as a limiting case) inside or on an ellipsoid centered at the maximum likelihood estimator. We prove this when p = 1 and report numerical calculations that support it when p > 1. The global maximum of the adjusted profile likelihood lies at infinity for any N .

Introduction

Maximum likelihood estimates of dynamic linear fixed effect models from short panels are known to be severely biased and inconsistent because of the incidental parameter problem first described by [START_REF] Neyman | Consistent estimates based on partially consistent observations[END_REF]; see [START_REF] Nerlove | Experimental evidence on the estimation of dynamic economic relations from a time series of cross-sections[END_REF][START_REF] Nerlove | Further evidence on the estimation of dynamic economic relations from a time series of cross-sections[END_REF] and [START_REF] Nickell | Biases in dynamic models with fixed effects[END_REF], for example. The failure of maximum likelihood has shifted interest towards GMM-based estimation in these models. An overview and discussion is provided by Arellano (2003, Chapter 6).

In an interesting paper, [START_REF] Lancaster | Orthogonal parameters and panel data[END_REF] argued for a return to likelihood-based inference and suggested a Bayesian resolution of the incidental parameter problem in the non-stationary AR(1) model. 1 The fixed effects are first orthogonalized (in the information sense) to the AR(1) model without covariates, [START_REF] Moreira | A maximum likelihood method for the incidental parameter problem[END_REF] maximal invariant likelihood differs from the adjusted profile likelihood.

Perhaps our main (and unexpected) finding is that one should look for a local maximizer of the adjusted likelihood/marginal posterior because the global maximum is reached at infinity. The reason for this is that the profile log-likelihood is log-quadratic in the parameters while the adjustment term is a polynomial, with positive coefficients, in the autoregressive parameters. Thus, for large enough autoregressive parameter values, the adjustment term dominates, and the adjusted profile log-likelihood must be re-increasing. This finding complements observations by [START_REF] Lancaster | Orthogonal parameters and panel data[END_REF], who noticed in Monte Carlo experiments that his marginal posterior was sometimes re-increasing. This phenomenon is not just a small sample problem or an artifact of an unbounded parameter space. The adjusted likelihood/marginal posterior has its global maximum at infinity for any sample size, and it may already be re-increasing in the stationary parameter region and reach its maximum at the boundary. Consequently, identification is not achieved by global maximization, and estimation based solely on the first order condition is incomplete.

To examine identification, we profile out the covariate slopes, which are unique given the autoregressive parameters, so the analysis reduces to that of the model without covariates.

In the AR(1) model, we show that the autoregressive parameter is identified as the unique local maximizer-or flat inflection point, as a limiting case-of the (large sample) adjusted profile likelihood in the interval where the unadjusted profile score is downward sloping. The (large sample) maximum likelihood estimate is the midpoint of this interval. Identification as a flat inflection point is an instance of first order lack of identification, in [START_REF] Sargan | Identification and lack of identification[END_REF] terminology. The leading case of interest is a unit root without deterministic trends, which is identified as a flat inflection point regardless of the number of time periods, in line with the rank deficiency of the expected Jacobian associated with the moment conditions of [START_REF] Ahn | Efficient estimation of models for dynamic panel data[END_REF]; see [START_REF] Alvarez | Robust likelihood estimation of dynamic panel data models[END_REF]. For the unit root case the flat inflection point property has been noted independently by [START_REF] Ahn | Likelihood based inference for dynamic panel data models[END_REF].

Identification in the AR(p) model via the adjusted likelihood is a harder problem. It appears that the autoregressive parameter vector is identified as the unique local maximizer or flat inflection point of the adjusted profile likelihood on a region bounded by a well-identified ellipsoid that is centered at the maximum likelihood estimate. We have no proof of this conjecture, but we report numerical calculations that support it.

The model parameters can be estimated by the local maximizer of the adjusted likelihood over a plug-in estimate of the ellipsoid (or, when no local maximizer exists, by the minimizer of the norm of the adjusted profile score). When the parameters are identifed as local maximizers, the estimator converges at the parametric rate, and its limit distribution is normal and correctly centered.

[3]

In Section 1, we present the model and derive the bias of the profile score. The adjusted likelihood is given in Section 2, where we also connect with [START_REF] Lancaster | Orthogonal parameters and panel data[END_REF], [START_REF] Arellano | Robust priors in nonlinear panel data models[END_REF], and [START_REF] Bun | Bias-corrected estimation in dynamic panel data models[END_REF]. Global properties of the adjusted likelihood are discussed in Section 3. In Section 4, we define the adjusted likelihood estimator and derive its main asymptotic properties. Section 5 discusses some simulation results. Technical proofs are collected in an appendix.

1 Bias of the profile score 1.1 Model and profile likelihood Suppose we observe a scalar variable y, the first p ≥ 1 lags of y, and a q-vector of covariates x (which may include lags), for N units i and T periods t. Assume that y it is generated by

y it = y ′ it-ρ + x ′ it β + α i + ε it , i = 1, ..., N ; t = 1, ..., T, (1.1)
where y it-= (y it-1 , ..., y it-p ) ′ and the ε it are identically distributed with mean 0 and variance σ 2 and are independent across i and t and also of x i ′ t ′ for all i ′ and t ′ . Let y 0 i = (y i(1-p) , ..., y i0 ) ′ , X i = (x i1 , ..., x iT ) ′ , and ε i = (ε i1 , ..., ε iT ) ′ . We make no assumptions about how (y 0 i , α i , X i ), i = 1, ..., N , are generated, except for the assumption that the covariates are not linearly dependent. The unknown parameters are θ = (ρ ′ , β ′ ) ′ , σ 2 , and α 1 , ..., α N . Let θ 0 and σ 2 0 be the true values of θ and σ 2 . Our interest lies in consistently estimating θ 0 under large N and fixed T asymptotics. We do not impose the stationarity condition on ρ 0 , i.e., we allow any ρ 0 ∈ R p .

Let

z it = (y ′ it-, x ′ it ) ′ , Y i-= (y i1-, ..., y iT -) ′ , Z i = (Y i-, X i )
, and y i = (y i1 , ..., y iT ) ′ , so that M y i = M Z i θ + M ε i where M = I T -T -1 ιι ′ and ι is a conformable vector of ones. We assume that N -1 N i=1 Z ′ i M Z i and its probability limit as N → ∞ are nonsingular. The Gaussian quasilog-likelihood, conditional on y 0 1 , ..., y 0 N and normalized by the number of observations, is given by

- 1 2N T N i=1 T t=1 log σ 2 + 1 σ 2 (y it -z ′ it θ -α i ) 2 + c,
where, here and later, c is a non-essential constant. Replacing the nuisance parameters α 1 , ..., α N and σ 2 with their MLE for given θ gives the (normalized) profile log-likelihood,

l(θ) = - 1 2 log 1 N N i=1 (y i -Z i θ) ′ M (y i -Z i θ) + c. [4]
The profile score, s(θ) = ∇ θ l(θ), has elements

s ρ j (θ) = N i=1 (y i -Z i θ) ′ M y i,-j N i=1 (y i -Z i θ) ′ M (y i -Z i θ) , j = 1, ..., p, s β j (θ) = N i=1 (y i -Z i θ) ′ M x i,j N i=1 (y i -Z i θ) ′ M (y i -Z i θ) , j = 1, ..., q,
where y i,-j is the jth column of Y i-and x i,j is the jth column of X i .

Bias of the profile score

It is well known that the profile score is asymptotically biased, that is,

plim N →∞ s(θ 0 ) = 0,
where, here and later, probability limits and expectations are taken conditionally, given (y 0 i , α i , X i ), i = 1, ..., N . Hence, θ 0 = arg max θ plim N →∞ l(θ) and the MLE, solving s(θ) = 0, is inconsistent.

The following calculation yields the bias of the profile score, either asymptotically or, when ε it ∼ N (0, σ 2 ), for fixed N . Rewrite (1.1) as

Dy i = Cy 0 i + X i β + ια i + ε i , i = 1, ..., N,
where D = D(ρ) and C = C(ρ) are the T × T and T × p matrices

D =              1 0 • •
• • • 0 -ρ p • • • -ρ 1 1              , C =                ρ p • •
0 • • • • • • 0 ρ p 0 • • • • • • • • • 0 . . . . . . . . . . . . . . . 0 • • • • • • • • • 0                . Then y 0 i y i = ξ i + F ε i , ξ i = y 0 i D -1 (Cy 0 i + X i β + ια i ) , F = 0 D -1 ,
and y i,-j = S j (ξ i + F ε i ), where S j = (0 T ×(p-j) : I T : 0 T ×j ), a selection matrix. Therefore,

plim N →∞ s ρ j (θ 0 ) = plim N →∞ N -1 N i=1 ε ′ i M S j (ξ 0i + F 0 ε i ) plim N →∞ N -1 N i=1 ε ′ i M ε i = E(ε ′ i M S j F 0 ε i ) E(ε ′ i M ε i ) = trM S j F 0 T -1 , plim N →∞ s β j (θ 0 ) = 0, [5]
where ξ 0i and F 0 are ξ i and F , evaluated at θ 0 . If, in addition, ε it ∼ N (0, σ 2 ), then

E[s ρ j (θ 0 )] = E N i=1 ε ′ i M S j (ξ 0i + F 0 ε i ) N i=1 ε ′ i M ε i = E N i=1 ε ′ i M S j F 0 M ε i N i=1 ε ′ i M ε i = E(ε ′ i M S j F 0 M ε i ) E(ε ′ i M ε i ) = trM S j F 0 T -1 , (1.2) E[s β j (θ 0 )] = 0,
by well-known properties of the normal distribution and the following geometric argument, which goes back to [START_REF] Fisher | The moments of the distribution for normal samples of measures of departure from normality[END_REF] and [START_REF] Geary | A general expression for the moments of certain symmetrical functions of normal samples[END_REF]. Let v ∼ N (0, σ 2 I g ) and let Q be a g × h matrix such that Q ′ Q = I h , so QQ ′ is idempotent. Write I g -QQ ′ as P P ′ , where

P ′ P = I g-h . Transform v into m = P ′ v, the radius r = (v ′ QQ ′ v) 1/2
, and the h -1 polar angles a of Q ′ v. Then the elements of (m ′ , r, a ′ ) ′ are independent. Therefore, for any g

× g matrix W , if A = v ′ QQ ′ W QQ ′ v and B = v ′ QQ ′ v,
then the ratio A/B depends on v only through a and hence is independent of B, which implies that

E(A) = E(A/B)E(B) and E(A/B) = E(A)/E(B). 3
The transition to (1.2) now follows from applying this property to the ratio Then plim N →∞ s(θ 0 ) = b(ρ 0 ); if, in addition, ε it ∼ N (0, σ 2 0 ), then E[s(θ 0 )] = b(ρ 0 ). Two remarks are worth making. First, the bias of the profile score, b(ρ 0 ), depends only on ρ 0 and T . It is independent of the initial observations, the fixed effects, and the covariates (recall that the expectations are conditioning on these). This is in sharp contrast with the bias of the ML (the within-group) estimator, which was first derived by [START_REF] Nickell | Biases in dynamic models with fixed effects[END_REF] for the first order autoregressive model under the assumption of stationarity of the initial observations. The bias of the MLE depends on the initial observations, the fixed effects, and the covariate values; it is smaller when the initial observations are more outlying relative to the stationary distributions. Second, the Nickell bias concerns a probability limit as N → ∞ whereas here, when the errors are normal, E[s(θ 0 )] = b(ρ 0 ) is a finite sample result holding for fixed N and T and may therefore be of independent interest in a time series setting.

N i=1 ε ′ i M S j F 0 M ε i N i=1 ε ′ i M ε i = ε ′ (I N ⊗ M )(I N ⊗ S j F 0 )(I N ⊗ M )ε ε ′ (I N ⊗ M )ε with v = ε = (ε ′ 1 , ..., ε ′ N ) ′ , QQ ′ = I N ⊗ M ,
For k = (k 1 , ..., k p ) ′ ∈ N p , let ρ k = p j=1 ρ k j j , and let τ = (1, ..., p) ′ . It is shown in the Appendix that

D -1 =      1 0 • • • 0 ϕ 1 . . . . . . . . . . . . . . . . . . 0 ϕ T -1 • • • ϕ 1 1      , ϕ t = τ ′ k=t (ι ′ k)! k 1 ! • • • k p ! ρ k , (1.3)
where the summation is over all k ∈ N p satisfying τ ′ k = t. Therefore, since

S j F = 0 0 D -1 j 0 where D -1 j is the leading (T -j) × (T -j) block of D -1 , we obtain b j (ρ) = - ι ′ D -1 j ι T (T -1) = - T -j-1 t=0 T -j -t T (T -1) ϕ t , j = 1, ..., p,
where ϕ 0 = 1. Hence, for j = 1, ..., p, b j (ρ) is a multivariate polynomial of degree Tj -1 in ρ 1 , ..., ρ p with non-positive coefficients.

2 Adjusted profile likelihood 2.1 Centered profile score and adjusted profile likelihood By construction, the centered (or adjusted) profile score,

s A (θ) = s(θ) -b(ρ),
is asymptotically unbiased, i.e., plim N →∞ s A (θ 0 ) = 0. In the Appendix, we show that the differential equation ∇ θ a(ρ) = b(ρ) has a solution which, up to an arbitrary constant of integration, is

a(ρ) = S∈S a S (ρ), a S (ρ) = - T -1 t=|S| T -t T (T -1) k∈K S :τ ′ k=t (ι ′ k -1)! k 1 ! • • • k p ! ρ k S S , (2.1)
where S is the collection of the non-empty subsets of {1, ..., p}; |S| is the sum of the elements of S; K S = {k ∈ N p |k j > 0 iff j ∈ S}; and ρ S = (ρ j ) j∈S and k S = (k j ) j∈S are subvectors of ρ and k determined by S.

We call the function l A (θ) = l(θ)a(ρ), [7] which has derivative ∇ θ l A (θ) = s A (θ), an adjusted profile log-likelihood. Every subvector ρ S of ρ contributes to l A (θ) an adjustment term, -a S (ρ), which takes the form of a multivariate polynomial in ρ j , j ∈ S, with positive coefficients that are independent of p.

(ρ) = - T -1 t=1 T -t T (T -1) ρ t-1 , a(ρ) = - T -1 t=1 T -t T (T -1)t ρ t .
With independent uniform priors on the reparameterized effects η i = α i e -(T -1)a(ρ) and on θ and log σ 2 , Lancaster's posterior density for ϑ 2l(θ) . Integrating over σ 2 gives f (θ|data) ∝ e -N (T -1)a(ρ) (Q 2 (θ)) -N (T -1)/2 and, hence,

= (θ ′ , σ 2 ) ′ is f (ϑ|data) ∝ σ -N (T -1)-2 exp -N (T -1)a(ρ) -Q 2 (θ)σ -2 /2 , where Q 2 (θ) = N i=1 (y i -Z i θ) ′ M (y i -Z i θ) ∝ e -
f (θ|data) ∝ e N (T -1)l A (θ) .
(2.2) Thus, the posterior density and the adjusted likelihood for θ are equivalent as a basis for inference. More generally, for any p and q, independent uniform priors on η 1 , ..., η N , θ, log σ 2 , with η i = α i e -(T -1)a(ρ) and a(ρ) as in (2.1), yield a posterior f (θ|data) that is related to l A (θ) as in (2.2).

Lancaster's choice of independent uniform priors on η i has its foundation in the AR(1) model without covariates, where η i is information orthogonal to ϑ and the posterior density (hence also e l A (θ) ) has an interpretation as a [START_REF] Cox | Parameter orthogonality and approximate conditional inference[END_REF] approximate conditional likelihood (see also Arellano, 2003, pp. 103-106). Information orthogonalization to a multidimensional parameter is generally not possible (see, e.g., Severini, 2000, pp. 340-342). Here, orthogonalization is not possible when the AR(1) model is augmented with covariates, as shown by Lancaster, or when the autoregressive order, p, is greater than one, as shown in the Appendix. From a bias correction perspective, however, information orthogonality is sufficient but not necessary. For any p and q, s A (θ) = 0 is an unbiased estimating equation, and the bias calculation underlying it is immune to the non-existence of orthogonalized fixed effects. [START_REF] Arellano | Robust priors in nonlinear panel data models[END_REF] approach shares the integration step with [START_REF] Lancaster | Orthogonal parameters and panel data[END_REF] but allows non-uniform priors on fixed effects or, equivalently, non-orthogonalized fixed effects. Of interest are bias-reducing priors, i.e., weighting schemes that deliver an integrated likelihood whose score equation has bias o(T -1 ) as opposed to the standard O(T -1 ). The linear dynamic model (with general p, q) illustrates an interesting result of Arellano and Bonhomme that generalizes the scope of uniform integration to situations where orthogonalization is impossible. For a given prior π i (α i |ϑ), the (normalized) log integrated likelihood is ρ) is a data-independent bias-reducing (in fact, bias-eliminating) prior in the sense of Arellano and Bonhomme. 6 Now, π i (α i |ϑ) ∝ e -(T -1)a(ρ) is equivalent to π i (η i |ϑ) ∝ 1, i.e., to a uniform prior on η i = α i e -(T -1)a(ρ) , leading to the same l I (ϑ). Arellano and Bonhomme (2009, Eq. (11)) give a necessary and sufficient condition for a uniform prior to be bias-reducing.

l I (ϑ) = 1 N T N i=1 log σ -T /2 exp - 1 2σ 2 T t=1 (y it -z ′ it θ -α i ) 2 π i (α i |ϑ) dα i + c. Choosing π i (α i |ϑ) ∝ e -(T -1)a(ρ) yields l I (ϑ) = - T -1 2T log σ 2 - T -1 T a(ρ) - Q 2 (θ) 2N T σ 2 + c. Profiling out σ 2 gives σ 2 (θ) = arg max σ 2 l I (ϑ) = Q 2 (θ)/(N (T -1)) and l I (θ) = max σ 2 l I (ϑ) = T -1 T l A (θ) + c, so l I (θ) and l A (θ) are equivalent. Because a(ρ) does not depend on true parameter values, π i (α i |ϑ) ∝ e -(T -1)a(
With ℓ i (ϑ, η i ) denoting i's log-likelihood contribution in a parametrization η i , the condition is that

plim N →∞ 1 N N i=1 ∇ η i (A -1 i B i ) = o(1) as T → ∞, (2.3)
where

A i = A i (ϑ, η i ) = -E ϑ,η i ∇ η i η i ℓ i (ϑ, η i ), B i = B i (ϑ, η i ) = E ϑ,η i ∇ ϑη i ℓ i (ϑ, η i ), and ∇ η i (A -1 i B i
) is evaluated at the true parameter values. When η i and ϑ are information orthogonal, B i = 0 and (2.3) holds. However, condition (2.3) is considerably weaker than parameter orthogonality. In the present model, when p > 1 or q > 0, and thus no orthogonalization is possible, it follows from our analysis and Arellano and Bonhomme ( 2009) that (2.3) must hold for η i = α i e -(T -1)a(ρ) . Indeed, as we show in the Appendix,

∇ η i (A -1 i B i ) = 0 (2.4)
6 Our focus is on estimating θ 0 but, noting the degrees-of-freedom correction in σ 2 (θ), the estimating equation derived from l I (ϑ) is also unbiased for σ 2 0 .

[9]

because A -1 i B i is free of η i . Finally, the centered profile score relates directly to [START_REF] Bun | Bias-corrected estimation in dynamic panel data models[END_REF]

. Observe that s(θ) = N i=1 Z ′ i M (y i -Z i θ)/Q 2 (θ). Further, M y i = M Z i θ + M ε i where θ is the within-group estimator, with residuals ε i satisfying N i=1 Z ′ i M ε i = 0. Therefore, solving s A (θ) = 0 is equivalent to solving θ -θ = N i=1 Z ′ i M Z i -1 b(ρ)Q 2 (θ).
(2.5)

In the AR(1) model, (2.5) corresponds to [START_REF] Bun | Bias-corrected estimation in dynamic panel data models[END_REF] proposal for bias-correcting the within-group estimator.

3 Global properties of the adjusted profile likelihood when N is large

At this point it is tempting to anticipate that θ 0 maximizes plim N →∞ l A (θ). However, as the following analysis shows, -a(ρ) dominates plim N →∞ l(θ) as ρ → ∞ in almost all directions and, as a result, plim

N →∞ l A (θ) is unbounded from above. Let h(θ) = ∇ θ ′ s(θ), c(ρ) = ∇ θ ′ b(ρ), and
L A (θ) = L(θ) -a(ρ), L(θ) = plim N →∞ l(θ), S A (θ) = S(θ) -b(ρ), S(θ) = plim N →∞ s(θ), H A (θ) = H(θ) -c(ρ), H(θ) = plim N →∞ h(θ). Using M (y i -Z i θ) = -M Z i (θ -θ 0 ) + M ε i , we have L(θ) = - 1 2 log plim N →∞ 1 N N i=1 (ε ′ i M ε i -2(θ -θ 0 ) ′ Z ′ i M ε i + (θ -θ 0 ) ′ Z ′ i M Z i (θ -θ 0 )) + c. Let b 0 = b(ρ 0 ) = S(θ 0 ) and note that plim N →∞ 1 N N i=1 Z ′ i M ε i = plim N →∞ 1 N N i=1 ε ′ i M ε i b 0 = σ 2 0 (T -1) b 0 . Hence, defining V 0 = V (θ 0 ) by plim N →∞ 1 N N i=1 Z ′ i M Z i = plim N →∞ 1 N N i=1 ε ′ i M ε i V 0 = σ 2 0 (T -1) V 0 , we can write L(θ) = - 1 2 log (1 -2(θ -θ 0 ) ′ b 0 + (θ -θ 0 ) ′ V 0 (θ -θ 0 )) + c by absorbing the term -1 2 log (σ 2 0 (T -1)) into c. As N → ∞, the ML estimator of θ converges in probability to θ M L = arg max θ L(θ) = θ 0 +V -1 0 b 0 and has asymptotic bias V -1 0 b 0 . This expression [10]
generalizes the bias calculations in [START_REF] Nickell | Biases in dynamic models with fixed effects[END_REF]. Note that (θ

0 -θ M L ) ′ V 0 (θ 0 -θ M L ) = b ′ 0 V -1 0 b 0 . Furthermore, L(θ) = - 1 2 log 1 -b ′ 0 V -1 0 b 0 + (θ -θ M L ) ′ V 0 (θ -θ M L ) + c, S(θ) = - V 0 (θ -θ M L ) 1 -b ′ 0 V -1 0 b 0 + (θ -θ M L ) ′ V 0 (θ -θ M L ) , H(θ) = - V 0 1 -b ′ 0 V -1 0 b 0 + (θ -θ M L ) ′ V 0 (θ -θ M L ) + 2S(θ)S(θ) ′ .
Note that L(•) and H(•) are even and S(•) is odd about θ M L and that

H(θ 0 ) = 2b 0 b ′ 0 -V 0 and H A (θ 0 ) = 2b 0 b ′ 0 -V 0 -c 0 , where c 0 = c(ρ 0 ). Since L(θ) is log-quadratic in θ and a(ρ) is a multivariate polynomial with negative coefficients, L A (θ) = L(θ) -a(ρ) is unbounded from above.
For example, if we put ρ = kr with r in the positive orthant of R p and let k → ∞, then the term -a(ρ) dominates and L A (θ) → ∞. 7 Hence, θ 0 = arg max θ L A (θ) and θ 0 has to be identified as a functional of L A (θ) other than its maximizer. Because S A (θ 0 ) = 0, we need to select θ 0 from the set of stationary points of L A (θ), that is, from the set of zeros of S A (θ). In general, this set is not a singleton. Indeed, whenever θ 0 is a local maximizer of L A (θ) (which will often be the case), L A (θ), being smooth and unbounded, must also have at least one local minimum. A moment-based estimation strategy based solely on solving s A (θ) = 0 is, therefore, incomplete.

AR(1) without covariates

Consider the AR(1) model without covariates, i.e., p = 1 and q = 0. Let

ζ 2 0 = (V 0 -b 2 0 ) /V 2 0 . Then, L(ρ) = - 1 2 log ζ 2 0 + (ρ -ρ M L ) 2 + c, S(ρ) = - ρ -ρ M L ζ 2 0 + (ρ -ρ M L ) 2 , H(ρ) = - ζ 2 0 -(ρ -ρ M L ) 2 (ζ 2 0 + (ρ -ρ M L ) 2 ) 2 , by absorbing -1 2 log V 0 into c. Note that ζ 2 0 = -1/H(ρ M L ). Recall that S(ρ) is odd about ρ M L = ρ 0 + b 0 /V 0 . The zeros of H(ρ) are ρ = ρ M L -ζ 0 and ρ = ρ M L + ζ 0 , so S(ρ) decreases on [ρ, ρ]
and increases elsewhere. Note that all of ρ, ρ, ρ M L , and ζ 0 are identified by S(•). Further, ρ M L and ζ 0 act as location and scale parameters of S(•). For any given ρ 0 , ρ M L and ζ 0 are determined by V 0 . As V 0 increases, |b 0 /V 0 | and ζ 0 decrease, that is, the bias of ρ M L decreases in absolute value, the length of [ρ, ρ] shrinks, and S(ρ) becomes steeper on (ρ, ρ).

7 It has already been observed that, in finite samples, a bias-adjusted likelihood may be re-increasing. Moreover, [START_REF] Lancaster | Orthogonal parameters and panel data[END_REF] presented graphs of a re-increasing marginal posterior log-density for ρ in an AR(1) model with time dummies. [START_REF] Arellano | Robust priors in nonlinear panel data models[END_REF], using a slightly different prior, presented a graph of a re-increasing integrated likelihood for an AR(1) model. Because l(θ) is log-quadratic for any N ≥ 1 and a(ρ) does not depend on the data, l A (θ) is re-increasing regardless of the sample size. [11] There is a sharp lower bound on V 0 . From y i,-1 = S 1 (ξ 0i + F 0 ε i ) and the independence between ξ 0i and ε i , we obtain

V 0 = plim N →∞ 1 N N i=1 y ′ i,-1 M y i,-1 σ 2 0 (T -1) = V LB 0 + V ξξ ,
where

V LB 0 = trF ′ 0 S ′ 1 M S 1 F 0 T -1 , V ξξ = plim N →∞ 1 N N i=1 ξ ′ 0i S ′ 1 M S 1 ξ 0i σ 2 0 (T -1)
.

So V 0 ≥ V LB 0
and this lower bound implies an upper bound on |b 0 /V 0 | and on the length of [ρ, ρ], and a lower bound on the steepness of S(ρ) on (ρ, ρ).

Lemma 1. V LB 0 is given by

V LB 0 = 1 T -1   T -2 j=0 (T -j -1) ρ 2j 0 - 1 T T -2 j=0 j k=0 ρ k 0 2   and satisfies (i) V LB 0 ≥ 2b 2 0 ; (ii) V LB 0 ≥ 2b 2 0 -c 0 with equality if and only if T = 2 or ρ 0 = 1. By Lemma 1, H(ρ 0 ) = 2b 2 0 -V 0 ≤ 0 and, hence, (ρ -ρ M L ) 2 = V 0 -b 2 0 V 2 0 ≥ b 2 0 V 2 0 = (ρ 0 -ρ M L ) 2 .
Therefore, ρ 0 ∈ [ρ, ρ]. Since S(ρ) is a rational function that vanishes at ±∞ and b(ρ) is a polynomial, S A (ρ) has finitely many zeros. Thus, because S A (ρ 0 ) = 0 and, by Lemma 1, H A (ρ 0 ) = 2b 2 0 -V 0c 0 ≤ 0, it follows that L A (ρ) has a local maximum or a flat inflection point at ρ 0 . Our main result for the AR(1) model without covariates is the uniqueness of such a point in [ρ, ρ], thereby identifying ρ 0 as a functional of L A (ρ).

Theorem 1. ρ 0 is the unique point in [ρ, ρ] where L A (ρ) has a local maximum or a flat inflection point.

Equivalently, ρ 0 is the unique point in [ρ, ρ] where b(ρ) approaches S(ρ) from below.

L A (ρ) has a flat inflection point at ρ 0 if and only if

V 0 = V LB 0 = 2b 2 0 -c 0 .
The latter equality holds if and only if T = 2 or ρ 0 = 1. The former holds if and only if V ξξ = 0, which requires M S 1 ξ 0i to be negligibly small for almost all i. The elements of S 1 ξ 0i are ρ j-1

0 y 0 i + α i j-1 k=1 ρ k-1 0 , j = 1, ..., T , so M S 1 ξ 0i = 0 if and only if y 0 i (1 -ρ 0 ) = α i .
Thus, when ρ 0 = 1 and there are no deterministic trends, i.e., α i = 0, L A (ρ) has a flat inflection point at ρ 0 for any T . When

ρ 0 = 1, V 0 = V LB 0 = 2b 2
0c 0 only when T = 2 and a very strong condition holds on the initial observations and the fixed effects, which is unlikely to hold in situations where a fixed effect modeling approach is called for. Thus, when ρ 0 = 1, except in quite special circumstances, ρ 0 is [12] the unique point in [ρ, ρ] where L A (ρ) attains a strict local maximum. Note that, when ρ 0 is a local maximizer of L A (ρ), it need not be the global maximizer on [ρ, ρ], which may instead be ρ. To see why this may happen, interpret the situation where L A (ρ) has a flat inflection point at ρ 0 as a limiting case of the property that L A (ρ) is re-increasing.

Figure 1 illustrates how ρ 0 is identified by L A (ρ) for two cases, each with T = 4. The figures on the left correspond to the case ρ 0 = .5 with V 0 = V LB 0 + V ξξ and V ξξ corresponding to stationary initial observations. Those on the right correspond to the unit root case without deterministic trends, i.e., ρ 0 = 1 and V 0 = V LB 0 . In each case, the bottom figures show S(ρ) (solid line) and b(ρ) (dashed line); the top figures show L(ρ) (solid line), -a(ρ) (dashed line), and L A (ρ) = L(ρ)a(ρ) (thick line). In all the figures, vertical lines indicate ρ, ρ 0 , and ρ, from left to right. In the case of ρ 0 = .5, ρ 0 is the unique local maximizer of L A (ρ) on [ρ, ρ]. Note that there is a second solution of S A (ρ) = 0 on [ρ, ρ], which corresponds to a local minimum of L A (ρ). In the unit root case, ρ 0 is the unique flat inflection point of L A (ρ) on [ρ, ρ]. ρ 0 + b 0 /V 0 . The proof of Theorem 1, as a by-product, shows that if T is even, then b 0 < 0; and, if T is odd, then b(ρ) decreases and has a unique zero at some point ρ u ∈ [-2, -1), so b 0 has the same sign as ρ uρ 0 .

AR(1) with covariates

In the AR(1) model with covariates, profiling out β yields a profile likelihood of ρ with essentially the same properties as in the AR(1) model without covariates. Let

β(ρ) = arg max β L A (ρ, β). Clearly, β(ρ) = arg max β L(ρ, β) = arg min β (θ -θ M L ) ′ V 0 (θ -θ M L
). Partition V 0 , V -1 0 , and b 0 as

V 0 = V 0ρρ V 0ρβ V 0βρ V 0ββ , V -1 0 = V ρρ 0 V ρβ 0 V βρ 0 V ββ 0 , b 0 = b 0ρ 0 ,
where V 0ρρ , V ρρ 0 , and b 0ρ are scalars. Using

V ρρ 0 = (V 0ρρ -V 0ρβ V -1 0ββ V 0βρ ) -1 , we have V 0ββ (β(ρ) -β M L ) = -V 0βρ (ρ -ρ M L ), min β (θ -θ M L ) ′ V 0 (θ -θ M L ) = (ρ -ρ M L ) 2 /V ρρ 0 , 1 -b ′ 0 V -1 0 b 0 = 1 -b 2 0ρ V ρρ 0 .
The first of these equations, together with V 0 (θ 0θ M L ) = -b 0 , yields β(ρ 0 ) = β 0 , so β 0 is identified whenever ρ 0 is. Profiling out β from L(ρ, β) gives the limiting profile log-likelihood of ρ as

L(ρ) = L(ρ, β(ρ)) = - 1 2 log ζ 2 0 + (ρ -ρ M L ) 2 + c
(slightly abusing notation), where ζ 2 0 is redefined as

ζ 2 0 = 1 -b 2 0ρ V ρρ 0 V ρρ 0 and 1 2 log V ρρ 0 is absorbed into c. Lemma 2. (V ρρ 0 ) -1 ≥ V LB 0 , with V LB 0
as defined earlier and given in Lemma 1.

In view of Lemma 2, we can invoke the result for the AR(1) model without covariates. If

ρ = ρ M L -ζ 0 and ρ = ρ M L + ζ 0 with ζ 0 redefined as indicated, then ρ 0 is the unique point in [ρ, ρ] where L A (ρ) = L(ρ) -a(ρ)
has a local maximum or a flat inflection point. By the proof of Lemma 2, the conditions under which ρ 0 is a flat inflection point of L A (ρ) are the same as in the AR(1) model without covariates. Note, also, that the presence of covariates does not affect the sign of the asymptotic bias of the MLE of ρ. Furthermore, it follows from the proof of Lemma 2 that the inclusion of covariates in the model cannot increase V ρρ 0 , so the magnitude of ρ M Lρ 0 = V ρρ 0 b 0ρ can only decrease relative to the model without covariates. This generalizes an observation of [START_REF] Phillips | Bias in dynamic panel estimation with fixed effects, incidental trends and cross section dependence[END_REF] in the stationary AR(1) model, and contrasts with [START_REF] Nickell | Biases in dynamic models with fixed effects[END_REF]Nickell ( , p. 1424).

[14]

AR(p)

Consider first the AR(p) model with p > 1 and without covariates, i.e., q = 0. Then

L(ρ) = - 1 2 log (1 + (ρ -ρ M L ) ′ W 0 (ρ -ρ M L )) + c, W 0 = V 0 1 -b ′ 0 V -1 0 b 0 , S(ρ) = - W 0 (ρ -ρ M L ) 1 + (ρ -ρ M L ) ′ W 0 (ρ -ρ M L ) , H(ρ) = - W 0 1 + (ρ -ρ M L ) ′ W 0 (ρ -ρ M L ) + 2S(ρ)S(ρ) ′ ,
where

-1 2 log(1 -b ′ 0 V -1 0 b 0 ) is absorbed into c. Note that W 0 = -H(ρ M L ), so W 0 is identified by L(•).
As in the AR(1) case, there is a lower bound on V 0 . Recalling that Y i-= (y i,-1 , ..., y i,-p ) and y i,-j = S j (ξ 0i + F 0 ε i ), where ξ 0i and ε i are independent, we have

V 0 = plim N →∞ 1 N N i=1 Y ′ i-M Y i- σ 2 0 (T -1) = V LB 0 + V ξξ
where V LB 0 and V ξξ have elements

V LB 0 jk = trF ′ 0 S ′ j M S k F 0 T -1 , (V ξξ ) jk = plim N →∞ 1 N N i=1 ξ ′ 0i S ′ j M S k ξ 0i σ 2 0 (T -1)
, for 1 ≤ j, k ≤ p. Hence, V 0 -V LB 0 is positive semi-definite, which we write as V 0 ≥ V LB 0 . When p ≥ T , while V 0 is nonsingular by assumption, rank(V LB 0 ) ≤ T -1 because S j F 0 = 0 for j ≥ T , which implies that V LB 0 jk = 0 whenever j ≥ T or k ≥ T . Thus, when p ≥ T , although V 0 can be arbitrarily close to V LB 0 , V 0 = V LB 0 . Further, when p ≥ T , b j (ρ) = 0 for j ≥ T because the sum defining b j (ρ) is empty, and c ij (ρ) = 0 for i + j ≥ T . Hence, when p ≥ T , V LB 0 -2b 0 b ′ 0 and V LB 0 -2b 0 b ′ 0 + c 0 have only zeros beyond their leading (T -1) × (T -1) blocks. A proof of generalizations of (i)-(ii) of Lemma 1 and Theorem 1 to the AR(p) model would be desirable but is more difficult than in the AR(1) case and our attempts failed. 8 Numerical computations, however, suggest that

V LB 0 ≥ 2b 0 b ′ 0 , V LB 0 ≥ 2b 0 b ′ 0 -c 0 , (3.1) rank(V LB 0 -2b 0 b ′ 0 + c 0 ) = min(p, T -2) if p j=1 ρ 0j = 1 or T < p + 2, p -1 else. (3.2)
8 A major difficulty is the rapidly increasing complexity of ϕ t as p increases. For example, ϕ t = ⌊t/2⌋ k=0

(t-k)! (t-2k)!k! ρ t-2k 1 ρ k 2 when p = 2.
Compare this with ϕ t = ρ t 1 when p = 1 and note that V LB 0 , b 0 , and c 0 depend in non-trivial ways on the ϕ t 's.

[15] Specifically, we computed the eigenvalues of V LB 0 -2b 0 b ′ 0 and V LB 0 -2b 0 b ′ 0 + c 0 for p = 2, 3, 4; T = 2, ..., 7; and all ρ 0 in a subset of R p chosen as follows. For p = 4, we put a square grid on the Cartesian product of the two triangles defined by

-1 ≤ γ 2 ≤ 1, γ 2 -1 ≤ γ 1 ≤ 1 -γ 2 , -1 ≤ γ 4 ≤ 1, γ 4 -1 ≤ γ 3 ≤ 1 -γ 4 , (3.3)
which is the stationary region of the lag polynomial γ(L)

= (1 -γ 1 L -γ 2 L 2 )(1 -γ 3 L -γ 4 L 2 ).
For each point on this grid and for each of the values m = 1, 2, 4, ρ 0 was calculated by equating the coefficients on both sides of mρ 01 Lρ 02 L 2ρ 03 L 3ρ 04 L 4 = mγ(L). For m = 1, the stationary region is covered, while for larger m a larger region is covered, though less densely.

In addition to (3.3) we set γ 4 = 0 for p = 3, and γ 3 = γ 4 = 0 for p = 2. The grid points on the region defined by (3.3) were spaced at intervals of .002 when p = 2, .02 when p = 3, and .1 when p = 4. We found that, uniformly over this numerical design, the eigenvalues of

V LB 0 -2b 0 b ′ 0 and V LB 0 -2b 0 b ′ 0 + c 0 are non-negative and the rank of V LB 0 -2b 0 b ′ 0 + c 0 is
as given by (3.2). These findings, while obviously not a proof, support (3.1) and (3.2), and we shall proceed under the assumption that (3.1) and (3.2) hold. 9

Because

V 0 ≥ V LB 0 , (3.1) implies that V 0 ≥ 2b 0 b ′ 0 and that H A (ρ 0 ) = 2b 0 b ′ 0 -V 0 -c 0 ≤ 0. Pre-and postmultiplication of V 0 ≥ 2b 0 b ′ 0 by b ′ 0 V -1 0 and V -1 0 b 0 gives b ′ 0 V -1 0 b 0 ≤ 1 2 ≤ 1 -b ′ 0 V -1 0 b 0 . Recalling that (ρ 0 -ρ M L ) ′ V 0 (ρ 0 -ρ M L ) = b ′ 0 V -1 0 b 0 , we see that (ρ 0 -ρ M L ) ′ W 0 (ρ 0 -ρ M L ) ≤ 1.
Therefore, if (3.1) and (3.2) hold, ρ 0 is a point in the ellipsoidal disk E = {ρ : (ρρ M L ) ′ W 0 (ρρ M L ) ≤ 1} where L A (ρ) has a local maximum or a flat inflection point. We approached the question of uniqueness of such a point numerically. For the same numerical design as above and with V 0 = V LB 0 , we applied the Newton-Raphson algorithm to find a stationary point of L A (ρ), starting at ρ M L and using the Moore-Penrose inverse of H A (ρ) whenever H A (ρ) is singular. Uniformly over this design, the algorithm was found to converge to ρ 0 , thus supporting the conjecture that ρ 0 is the unique point in E where L A (ρ) has a local maximum or a flat inflection point. 10 In the AR(p) model with covariates, β can be profiled out of L A (θ) just as in the AR(1) model. Here, again, β 0 = β(ρ 0 ). Lemma 2 continues to hold for p > 1. Hence, if ρ 0 is identified in the model without covariates in the way we suggested, then it is identified in the model with covariates in exactly the same way, now with E defined through

W 0 = (1 -b ′ 0ρ V ρρ 0 b 0ρ ) -1 V 0ρρ
, in obvious notation. 9 The same computations but with T = 8, 9, 10 further supported the conclusions. Here, however, when m = 4 and p = 3, 4 the computations are numerically less stable because the polynomial terms may be extremely large and their sum numerically imprecise. It is natural to expect, however, that (3.1) and (3.2) hold for all T if they hold for the small values of T considered.

10 Computations with T = 8, 9, 10 gave the same results except in certain cases with m = 4 and p = 3, 4 where the algorithm failed to converge because the rank of H A (ρ) was underestimated.

[16]

Estimation

For a given ρ, we define

β AL (ρ) = arg max β l A (ρ, β) = N i=1 X ′ i M X i -1 N i=1 X ′ i M (y i -Y i-ρ)
as the adjusted likelihood (AL) estimator of β 0 . (Note that β AL (ρ) coincides with the MLE of β 0 for a given ρ.) The unadjusted and adjusted profile log-likelihoods for ρ are l(ρ) = l(ρ, β AL (ρ)) and l A (ρ) = l(ρ)a(ρ). Let s(ρ), s A (ρ), h(ρ), and h A (ρ) be the corresponding profile scores and Hessians. Let W = -h( ρ M L ), where ρ M L is the MLE of ρ 0 , and let

E = {ρ : (ρ -ρ M L ) ′ W (ρ - ρ M L ) ≤ 1}.
We define the AL estimator of ρ 0 as

ρ AL = arg min ρ∈ b E s ′ A (ρ)s A (ρ) s.t. h A (ρ) ≤ 0,
that is, as the strict local maximizer of l A (ρ) on the interior of E if such a maximizer exists and otherwise as the minimizer of the norm of s A (ρ) on E. The AL estimator of β 0 , then, is

β AL = β AL ( ρ AL ). Let N → ∞. Then, l A (ρ) converges to L A (ρ) uniformly in ρ since -a(ρ) is nonstochastic and sup ρ |l(ρ) -L(ρ)| = o p (1). Further, ρ M L p → ρ M L , W p → -H(ρ M L ) = W 0 , and E p → E in the sense that Pr[ρ ∈ E] → 1 {ρ∈E} for any ρ not on the boundary of E. It follows that θ AL = ( ρ ′ AL , β ′ AL ) ′ p → θ 0 . When H A (θ 0 ) is nonsingular,

by the usual argument involving a Taylor series expansion of s

A (θ) around θ 0 , √ N ( θ AL -θ 0 ) d → N (0, Ω) , Ω = H A (θ 0 ) -1 (V 0 -b 0 b ′ 0 )H A (θ 0 ) -1 . (4.1)
The asymptotic variance can be estimated in the usual way. We have not investigated the limit distribution of θ AL in the situation where H A (θ 0 ) is singular. Presumably this could be done along the lines of [START_REF] Rotnitzky | Likelihood-based inference with singular information matrix[END_REF].

Simulations

We used simulations to examine the finite sample properties of the AL estimator in AR(1) and AR(2) models without covariates and in an AR(1) model with a single stationary and exogenous regressor. We compared the AL with the one-step GMM estimator of [START_REF] Arellano | Some tests of specification for panel data: Monte Carlo evidence and an application to employment equations[END_REF], which leaves the initial observations unrestricted. 11 In all the designs, we set N = 100, generated

11 For the GMM estimator we used all the moments that follow from the exogeneity of x it , i.e. E [x ij (ε itε it-1 )] = 0, j = 1, ..., T ; t = 2, ..., T ; and E [y it-j (ε itε it-1 )] = 0, j = 2, ..., t; t = 2, ..., T . It should be noted that the GMM estimator remains consistent under heteroskedasticity while the bias calculation underlying the AL estimator requires time series homoskedasticity.

[17] ε it and α i as N (0, 1) variates, 12 and chose ρ 0 in the interior of the stationary region, which implies that y it is eventually stationary as t → ∞. We varied the information content of the data through the initial observations. Let µ i = lim t→∞ E(y it-|α i ) and Σ i = lim t→∞ Var(y it-|α i ), so, if y 0 i was drawn from the stationary distribution, we would just have µ i = E(y 0 i |α i ) and Σ i = Var(y 0 i |α i ). Let G i G ′ i = Σ i be the Cholesky factorization of Σ i . We set y 0 i = µ i + ψG i ι for some chosen scalar ψ ≥ 0, which is a p-variate version of setting the initial observations ψ standard deviations away from the stationary mean. So ψ controls the outlyingness of the initial observations relative to the stationary distributions. All else being equal, V 0 increases in ψ and V 0 → V LB 0 as ψ → 0, so the data carry less information as ψ gets smaller. The effect of strong inlying observations (small ψ) on the informativeness of the data is stronger when T is small because it takes time to revert to the stationary distribution. The effect of ψ is vanishingly small as ρ 0 moves to the boundary of the stationary region. 13 We set ψ = 0, 1, 2 in the AR(1) models and ψ = .3, 1, 2 in the AR(2) model. 14 In the models without a covariate, µ i and Σ i follow immediately from α i and ρ 0 . In the AR(1) model with a covariate, x it was generated by a stationary AR(1), x it = δα i + γx it-1 + u it with u it ∼ N (0, σ 2 u ) and x i0 drawn from the stationary distribution. Here,

µ i = α i 1 -ρ 0 1 + δβ 0 1 -γ , Σ i = 1 1 -ρ 2 0 1 + β 2 0 1 -γ 2 1 + γρ 0 1 -γρ 0 σ 2 u .
We set δ = γ = σ u = .5 and β 0 = 1ρ 0 , inducing dependence between the covariate and the fixed effect, and keeping the long-run multiplier of x on y constant at unity across designs, as in [START_REF] Kiviet | On bias, inconsistency, and efficiency of various estimators in dynamic panel data models[END_REF]. We varied ψ, T , and ρ 0 . Tables 1-3 present Monte Carlo estimates, based on 10,000 replications at each design point, of the bias and standard deviation (STD) of the AL and GMM estimators and of the coverage rates of the asymptotic and bootstrap 95% confidence intervals (CI A .95 and CI B .95 ). 15 In the AR(1) model (Table 1), for all designs ρ AL has less bias and less standard deviation than ρ GMM . When ρ 0 = .5, the AL estimator is virtually unbiased except when ψ = 0 and T = 2, in which case ρ 0 is weakly identified. 16 Both estimators deliver 95% confidence intervals with broadly correct coverage, although the coverage errors are somewhat larger for ρ GMM , 12 The AL estimator is invariant with respect to the fixed effects, but the GMM estimator is not. 13 In the AR(1) model without covariates, for example,

S 1 ξ 0i = g 0 (y i0 -µ i ) + ιµ i with g 0 = (1, ρ 0 , ..., ρ T -1 0 ) ′ , hence V ξξ = g ′ 0 M g0 (1-ρ 2 0 )(T -1) ψ 2 . As ρ 0 ↑ 1, V ξξ →
0 for any fixed ψ. 14 In the AR(2) model, when ψ = 0 the weighting matrix used by GMM is singular for all T and the likelihood function does not depend on ρ 2 for T = 2. Therefore, we set ψ > 0 although it is possible to work with ψ = 0 for T ≥ 3 by dropping the earliest observation from the set of instruments.

15 To compute CI B .95 , we used the bootstrap percentile method with 39 bootstrap samples formed by randomly drawing N units i with replacement from {1, ..., N }. A small number of bootstrap draws suffices for studying coverage rates via simulation.

16 When ψ = 0 and T = 2, l(ρ) depends on ρ 0 and ρ only through ρ 0ρ, so the distribution of ρ AL is

[18] [20] [21]

ψ T ρ 0 ρ AL ρ GMM ρ AL ρ GMM ρ AL ρ GMM ρ AL ρ GMM 0 2 .5 -.
where they increase in T . When ρ 0 is increased to .95, ρ GMM has a large negative bias and its standard deviation increases considerably. The AL estimator also deteriorates, but much less.

The confidence intervals based on ρ AL continue to provide good coverage, with the bootstrap being slightly superior, while those based on ρ GMM have large coverage errors for all ψ, especially when T ≥ 8. Turning to the AR(2) model (Table 2), the AL and GMM estimators perform well in terms of bias (for the chosen values of ρ 0 ), although there is a non-negligible bias when T = 2 and also when T = 4 and the initial observations are strong inlyers. As N = 100, the probability that the adjusted likelihood has no local maximum in the relevant region is fairly large when both T and ψ are small. In most designs, ρ AL has smaller standard deviation than ρ GMM , with the difference decreasing in T and ψ. For both estimators, the confidence intervals have very reasonable coverage, with those based on ρ AL being slightly more accurate in most designs.

The results for the AR(1) model with a covariate (Table 3) follow a similar pattern as in the AR(1) model without covariates. The coefficient on the exogenous variable, β 0 , is generally estimated with small bias by both estimators. When (ρ 0 , β 0 ) = (.5, .5), the biases and standard deviations are fairly small and, here, the AL estimator does not uniformly dominate GMM. When T = 2, θ GMM is less biased than θ AL while the opposite holds when T ≥ 4. When (ρ 0 , β 0 ) is changed to (.95, .05), the bias of β GMM increases only moderately, but ρ GMM is now heavily biased. As a result, GMM-based confidence intervals for ρ 0 tend to have poor coverage while those for β 0 have broadly correct coverage. Compared to GMM, the AL estimator is much less biased and yields confidence intervals with approximately correct coverage, especially when computed by bootstrapping.

Appendix

Proof of equations (1.3). For arbitrary ρ 1 , ..., ρ T -1 , D and its inverse are

D =      1 0 • • • 0 -ρ 1 . . . . . . . . . . . . . . . . . . 0 -ρ T -1 • • • -ρ 1 1      , D -1 =      1 0 • • • 0 φ 1 . . . . . . . . . . . . . . . . . . 0 φ T -1 • • • φ 1 1      ,
where φ 1 , ..., φ T -1 are recursively obtained as φ 1 = ρ 1 and φ j = ρ j + j-1 k=1 φ k ρ j-k , j = 2, ..., T -1. Recursive substitution gives

φ j = k 1 +2k 2 +•••+jk j =j (k 1 + ... + k p )! k 1 ! • • • k p ! ρ k 1 1 ρ k 2 2 • • • ρ k j j .
equivariant under translations of ρ 0 and, hence, the bias, standard deviation, and coverage rates are independent of ρ 0 .

[22]

Putting

ρ p+1 = • • • = ρ T -1 = 0 gives φ j = ϕ j .
Proof of equations (2.1). For j = 1, ..., p, let S j = {S ∈ S|j ∈ S}. Group terms by S ∈ S j to write b j (ρ)dρ j = S∈S j B j,S (ρ) + c, where

B j,S (ρ) = - T -j-1 t=0 T -j -t T (T -1) k∈K j,S :τ ′ k=t (ι ′ k)! k 1 ! • • • (k j + 1)! • • • k p ! ρ j ρ k S S and K j,S = {k ∈ N p | for all j ′ = j, k j ′ > 0 iff j ′ ∈ S} ⊃ K S . A change of variable from k j + 1 to k j gives B j,S (ρ) = - T -j-1 t=|S|-j T -j -t T (T -1) k∈K S :τ ′ k=t+j (ι ′ k -1)! k 1 ! • • • k p ! ρ k S S ,
where the lower limit in the first sum changed from 0 to |S|j because, when t < |S|j, no k ∈ K S satisfies τ ′ k = t + j. A further change of variable from t + j to t gives B j,S (ρ) = a S (ρ), with a S (ρ) as defined in (2.1). Therefore,

b j (ρ) = ∇ ρ j S∈S j a S (ρ) = ∇ ρ j S∈S a S (ρ) = ∇ ρ j a(ρ)
which completes the proof.

Proof of (2.4). In the parameterization η i = α i e -(T -1)a(ρ) , we have

ℓ i (ϑ, η i ) = - 1 2 log σ 2 - 1 2T σ 2 T t=1 (y it -z ′ it θ -η i e (T -1)a(ρ) ) 2 + c, ∇ η i ℓ i (ϑ, η i ) = e (T -1)a(ρ) T σ 2 (y i -Z i θ -η i e (T -1)a(ρ) ι) ′ ι, and
E ϑ,η i ∇ η i η i ℓ i (ϑ, η i ) = -σ -2 e 2(T -1)a(ρ) , E ϑ,η i ∇ σ 2 η i ℓ i (ϑ, η i ) = 0, E ϑ,η i ∇ θη i ℓ i (ϑ, η i ) = -σ -2 e (T -1)a(ρ) η i (T -1)b(ρ)e (T -1)a(ρ) + E ϑ,η i Z ′ i ι/T . The jth column of Y i-is y i,-j = S j (ξ i + F ε i ) , so the jth element of E ϑ,η i Y ′ i-ι is E ϑ,η i y ′ i,-j ι = ι ′ S j ξ i = ι ′ D -1 j ιη i e (T -1)a(ρ) + m j , m j = ι ′ S j y 0 i D -1 (Cy 0 i + X i β) . Hence, E ϑ,η i Z ′ i-ι/T = -η i (T -1)b(ρ)e (T -1)a(ρ) + m, [23]
where m = (m 1 , ..., m p , ι ′ X i /T ) ′ is free of η i . Consequently,

A -1 i B i = -e -(T -1)a(ρ) m 0 and ∇ η i (A -1 i B i ) = 0.
Proof that no orthogonalization exists when p > 1. In the original parameterization, if l i (ϑ, α i ) is i's log-likelihood contribution, we have

E ϑ,α i ∇ α i α i l i (ϑ, α i ) = -σ -2 , E ϑ,α i ∇ σ 2 α i l i (ϑ, α i ) = 0, E ϑ,α i ∇ θα i l i (ϑ, α i ) = -σ -2 E ϑ,α i Z ′ i ι/T,
and so, by the preceeding proof,

A -1 i B i = - E ϑ,α i Z ′ i ι/T 0 = - -(T -1)b(ρ)α i + m 0 .
Suppose some reparameterized fixed effect, say ζ i , is information orthogonal to ϑ. Then

α i = α i (ϑ, ζ i ) must satisfy the differential equation ∇ ϑ α i = A -1 i B i , that is, ∇ ρ j α i = (T -1)b j (ρ)α i -m j , j = 1, . . . , p, (A.1) ∇ β j α i = -m p+j , j = 1, . . . , q, (A.2)
and ∇ σ 2 α i = 0. We show that these equations are inconsistent. Suppose q > 0. Then (A.1) implies ∇ ρ j β j ′ α i = -∇ β j ′ m j , which is generally non-zero, while (A.2) implies ∇ ρ j β j ′ α i = 0, so the equations are inconsistent. Suppose q = 0. Then

m j = ι ′ S j I p D -1 C y 0 i , j = 1, . . . , p,
and, because ∇ ρ j ′ b j (ρ) = ∇ ρ j ′ ρ j a(ρ) = ∇ ρ j b j ′ (ρ), (A.1) will be inconsistent if ∇ ρ j ′ m j = ∇ ρ j m j ′ for some j, j ′ . Take j = p and j ′ = p -1. The first element of y 0 i appears in m p and m p-1 with coefficients γ p = 1 + ρ p T -p-1 t=0 ϕ t and γ p-1 = ρ p T -p t=0 ϕ t , respectively. Differentiating gives

∇ ρ p-1 γ p = ρ p T -p-1 t=0 ∇ ρ p-1 ϕ t = ρ p T -p t=1 ∇ ρp ϕ t , ∇ ρp γ p-1 = ρ p T -p t=1 ∇ ρp ϕ t + T -p k=0 ϕ t ,
using ϕ 0 = 1 and ∇ ρ p-1 ϕ t = ∇ ρp ϕ t+1 . The latter follows from differentiating ϕ t and a change of variable from k p-1 -1 to k p-1 , giving

∇ ρ p-1 ϕ t = τ ′ k=t-p+1 (ι ′ k + 1)! k 1 ! • • • k p ! ρ k ,
which is invariant under a unit shift of p and t. Therefore, ∇ ρ p-1 γ p = ∇ ρp γ p-1 , and (A.1) is inconsistent.

[24]

Proof of Lemma 1.

Let A = S 1 F 0 and B = ∇ ρ 0 A. Then b 0 = - ι ′ Aι T (T -1) , c 0 = - ι ′ Bι T (T -1) , V LB 0 = trA ′ M A T -1 = T trAA ′ -ι ′ AA ′ ι T (T -1) . Hence, V LB 0 ≥ 2b 2 0 and V LB 0 ≥ 2b 2 0 -c 0 if and only if T trAA ′ -ι ′ AA ′ ι - 2 (ι ′ Aι) 2 T (T -1) ≥ 0, (A.3) T trAA ′ -ι ′ AA ′ ι - 2 (ι ′ Aι) 2 T (T -1) -ι ′ Bι ≥ 0.
(A.4)

The matrix A = A T is

A = 0 0 1 0 , T = 2, A = A T -1 0 a ′ T 0 , a T = (ρ T -2 , ρ T -3 , ..., 1) ′ , T > 2,
where the subscript on ρ is omitted. By recursion, it can be deduced that as stated in the lemma. Now let r > 0 and use the equalities just obtained to see that if (A.4) holds for ρ = r, then (A.3) holds for ρ = r and (A.3) and (A.4) hold for ρ = -r, with strict inequalities for T ≥ 3. Hence, we only need to show that (A.4) holds for ρ ≥ 0, with equality if and only if T = 2 or ρ = 1. Write (A.4) as Q T ≥ 0. Because Q 2 = 0, to show that (A.4) holds, it suffices to show that ∆Q T ≥ 0 for T ≥ 2, where ∆ (•) T = (•) T +1 -(•) T . Write ∆Q T as

ι ′ Aι =
∆Q T = ∆ T trAA ′ -ι ′ AA ′ ι -2 (ι ′ Aι) 2 T (T -1) -ι ′ Bι T = (trAA ′ ) T +1 -2 (ι ′ Aι) 2 T +1 T (T + 1) + 2 (ι ′ Aι) 2 T T (T -1) -∆ (ι ′ Bι) T + T ∆ (trAA ′ ) T -∆ (ι ′ AA ′ ι) T [25]
and denote the quantities in braces as τ 1 , τ 2 , and τ 3 . Using T (T + 1) /2 = T -1 i=0 (Ti), we have

R =      1 1 • • • 1 ρ 2 ρ 2 • • • ρ 2 . . . . . . . . . . . . ρ 2T -2 ρ 2T -2 • • • ρ 2T -2      -      1 ρ • • • ρ T -1 ρ ρ 2 • • • ρ T . . . . . . . . . . . . ρ T -1 ρ T • • • ρ 2T -2      .
Consider the principal minors of R. Those of order 1 are 0; those of order 2 are det 0 ρ 2iρ i+j ρ 2jρ i+j 0 = ρ i+j ρ jρ i 2 ≥ 0, 0 < i < j < T, given ρ ≥ 0; and those of order greater than 2 are 0 because R is the sum of two matrices of rank 1 and, hence, rank (R) ≤ 2. Therefore, R is positive semi-definite and τ 1 ≥ 0. Furthermore, T (T -1) . Note that τ 2 is a polynomial of degree 2T -4 in ρ. When T = 2 or ρ = 1, τ 2 = 0. When ρ = 1,

τ
(ι ′ Aι) T = T -1 -T ρ + ρ T (1 -ρ) 2 , ∆ (ι ′ Bι) T = 1 -ρ T -T ρ T -1 + T ρ T (1 -ρ) 2 , [26]
and so

τ 2 = d T -1 -T ρ + ρ T 2 (1 -ρ) 4 - 1 -ρ T -T ρ T -1 + T ρ T (1 -ρ) 2 .
For T > 2, lim ρ→1 τ 2 (1ρ) -2 = 1 72 T (T -1) (T -2) (T + 1) > 0 and, therefore, τ 2 = (1ρ) 2 P (T, ρ), where P (T, ρ) is a polynomial of degree 2T -6. If all coefficients of P (T, ρ) =

2T -6 j=0 p j ρ j are positive, we conclude that τ 2 ≥ 0. Write τ 2 =

2T -4 j=0 q j ρ j , where q j is found as Equating the coefficients of τ 2 and (1ρ) 2 P (T, ρ) gives p k = k j=0 (k + 1j) q j . To show that p k > 0 for 0 ≤ k ≤ 2T -6, we only need to show that p k > 0 for k up to T -2 because for larger k, q k > 0 and so p k increases in k. For k up to min(T -2, 2T -6), we obtain

q j =
p k = d 12 (k + 1) (k + 2) (k + 3) (T -1) (2T -k -2) + k 10 (k -1) -T (T -1)
and, hence, p k > 0 because either k < T -2, implying 2Tk -2 > T , or k = T -2 ≤ 2T -6, implying T ≥ 4 and k ≥ 2. Therefore, τ 2 ≥ 0. This establishes Q T ≥ 0, that is, (A.4). Recall that Q 2 = 0 and note that ρ = 1 implies τ 1 = τ 2 = τ 3 = 0 and, hence, Q T = 0. Therefore, Q T = 0 if T = 2 or ρ = 1. If T ≥ 2 and ρ = 1, then ∆Q T > 0 because τ 3 > 0 when T = 2 and τ 2 > 0 when T > 2. Therefore, Q T = 0 only if T = 2 or ρ = 1.

Proof of Theorem 1. L A (ρ) having a local maximum or a flat inflection point at ρ 0 is equivalent to b(ρ) approaching S(ρ) from below as ρ approaches ρ 0 from the left. We will write this as b(ρ) ↑ S(ρ) at ρ 0 , and show that b(ρ) ↑ S(ρ) on [ρ, ρ] at most once. From [27]

∇ ρ H(ρ) = 2(ρ -ρ M L ) (3ζ 2 0 -(ρ -ρ M L ) 2 ) (ζ 2 0 + (ρ -ρ M L ) 2 )

  and W = I N ⊗ S j F 0 . To summarize the results so far regarding the bias of the profile score, let b(ρ) = (b 1 (ρ), ..., b p+q (ρ)) + 1, ..., p + q.

Figure!

  Figure 1: Identification

  where u = (T, T -1, ..., 1) ′ and

  i -1) (Tj -1) ρ i+ji -1) (j + 1) ρ j , where d = 2

d 6 {

 6 (j + 1) (j (j -1) + 6 (T -1) (Tj -1)) -3jT (T -1)} , j ≤ T -2, d 6 (2Tj -1) (2Tj -2) (2Tj -3) ,T -1 ≤ j.

3t

  it follows that S(ρ) is strictly concave on [ρ, ρ M L ] and strictly convex on [ρ M L , ρ]. Because ϕ t = ρ t , b(ρ) and its first two derivatives are b(ρ) = -(T -1t) T (T -1) ρ t-1 , d(ρ) = -T -2 t=2 t (t -1) (T -1t) T (T -1) ρ t-2 .

Table 1 :

 1 Simulations for the AR(1) model, N = 100

	bias	STD	CI A .95	CI B .95

Table 2 :

 2 Simulations for the AR(2) model, N = 100 AL ρ GMM ρ AL ρ GMM ρ AL ρ GMM ρ AL ρ GMM

	bias	STD	CI A .95	CI B .95
	ψ ρ .3 T ρ 0 2 .6 -.141 .2 -.129 1 2 .6 -.141 .2 -.124 2 2 .6 -.147 .2 -.130 .3 4 .6 -.071 -.279 .123 -.267 -.748 -.267 -.316 -.266 -.258 .2 -.033 -.125 .098 1 4 .6 -.001 -.043 .123 .2 -.002 -.023 .094 2 4 .6 .008 -.011 .082 .2 .003 -.007 .072 .3 8 .6 -.016 -.102 .065 .2 -.009 -.053 .056 1 8 .6 .006 -.033 .063 .2 .002 -.015 .052 2 8 .6 .000 -.011 .037 .2 .000 -.003 .037 .3 24 .6 .000 -.028 .024 .2 .000 -.018 .024 1 24 .6 .000 -.020 .022 .2 .000 -.011 .022 2 24 .6 .000 -.011 .020 .2 .000 -.003 .019 .3 2 1 -.144 -.265 -.2 -.083 -.640 1 2 1 -.145 -.265 -.2 -.073 -.234 2 2 1 -.147 -.266 -.2 -.078 -.165 .3 4 1 -.070 -.235 .118 -.2 .004 -.021 .087 1 4 1 .009 -.030 .115 -.2 .006 -.006 .081 2 4 1 .003 -.008 .059 -.2 .001 -.003 .060 .3 8 1 -.009 -.068 .059 -.2 .002 -.015 .051 1 8 1 .002 -.024 .051 -.2 .001 -.001 .042 2 8 1 -.001 -.009 .033 .001 .002 .033 -.2 .3 24 1 -.001 -.018 .022 -.2 .000 -.008 .022 1 24 1 -.001 -.015 .021 -.2 .000 -.004 .021 2 24 1 -.001 -.010 .019 -.2 .000 .000 .019		-.882 -.959 -.881 -.909 -.880 -.890 .280 .900 .133 .936 .120 .945 .089 .959 .064 .967 .064 .965 .088 .941 .056 .949 .055 .965 .045 .963 .037 .953 .036 .944 .025 .956 .023 .950 .023 .951 .022 .945 .020 .950 .019 .946 -.882 -.969 -.880 -.928 -.881 -.901 .256 .902 .084 .966 .098 .951 .078 .965 .051 .966 .059 .949 .072 .953 .045 .965 .048 .966 .041 .950 .034 .948 .034 .944 .023 .950 .022 .947 .022 .948 .021 .947 .020 .950 .019 .946	.920 .930 .945 .940 .916 .930 .920 .920 .916 .927 .918 .922 .782 .914 .808 .932 .921 .945 .931 .950 .939 .946 .943 .949 .773 .941 .833 .943 .897 .948 .929 .947 .932 .947 .943 .945 .795 .949 .869 .945 .859 .950 .914 .946 .913 .951 .944 .948 .918 .934 .968 .945 .913 .927 .921 .931 .918 .930 .918 .919 .808 .906 .934 .945 .930 .947 .939 .944 .940 .945 .941 .945 .830 .945 .926 .945 .916 .947 .945 .944 .936 .946 .945 .947 .872 .950 .927 .946 .898 .948 .940 .947 .920 .947 .947 .948	.977 .983 .977 .977 .977 .978 .698 .763 .895 .925 .933 .942 .536 .695 .816 .910 .913 .942 .607 .755 .732 .865 .856 .939 .977 .992 .979 .978 .977 .979 .721 .937 .908 .943 .936 .943 .653 .911 .856 .946 .922 .947 .770 .906 .818 .935 .879 .949
	'-' indicates non-existence of the corresponding population moment.

Table 3 :

 3 Simulations for the AR(1) model with a covariate, N = 100 GMM θ AL θ GMM θ AL θ GMM θ AL θ GMM

				bias	STD	CI A .95	CI B .95
	ψ T θ 0 θ 0 θ AL 2 .5 -.092 -.051 .266 .5 -.027 -.015 .249 1 2 .5 .019 .001 .269 .5 -.001 .000 .258 2 2 .5 .030 .002 .182 .5 -.010 -.002 .261 0 4 .5 .014 -.093 .141 .5 .004 -.008 .127 1 4 .5 .012 -.069 .119 .010 .126 .5 -.001 2 4 .5 .001 -.024 .063 .5 .000 .009 .129	-.900 -.965 -.939 -.968 -.966 -.949 .119 .957 .124 .959 .102 .968 .123 .951 .062 .968 .128 .944	.882 .943 .943 .945 .910 .947 .950 .939 .940 .949 .941 .942 .845 .945 .941 .947 .883 .948 .943 .949 .921 .944 .940 .941	.945 .962 .945 .949 .945 .941 .737 .946 .802 .947 .890 .941
	0	8	.5 -.001 -.053 .051 .5 .001 .004 .074	.050 .965 .075 .943	.813 .948 .943 .946	.582 .945
	1	8	.5 -.001 -.047 .045 .5 .001 .011 .075	.047 .959 .076 .940	.837 .945 .938 .943	.627 .937
	2	8	.5 -.001 -.026 .035 .5 .001 .011 .075	.036 .950 .075 .948	.884 .939 .945 .949	.754 .943
	0 24	.5 -.001 -.032 .020 .5 .000 .007 .038	.020 .948 .038 .950	.646 .943 .945 .950	.323 .941
	1 24	.5 .5	.000 -.031 .019 .000 .009 .038	.020 .946 .038 .950	.652 .945 .942 .952	.343 .937
	2 24	.5 .5	.000 -.026 .017 .000 .009 .038	.018 .947 .038 .949	.700 .948 .943 .952	.416 .935
	0 1 2 0 1	2 .95 -.145 -.904 .266 .05 -.007 -.030 .234 2 .95 -.124 -.975 .266 .05 .004 .046 .240 2 .95 -.068 -.534 .264 .05 .007 .063 .246 4 .95 -.086 -.641 .124 .05 -.002 -.013 .121 4 .95 -.067 -.595 .123 .05 .004 .041 .122	-.879 -.974 -.891 -.974 -.910 -.972 .253 .892 .115 .973 .252 .906 .114 .973	.809 .930 .977 .947 .811 .937 .977 .941 .842 .948 .968 .943 .267 .906 .944 .948 .312 .930 .938 .948	.926 .991 .923 .989 .936 .980 .123 .955 .148 .943
	2	4 .95 -.020 -.338 .122 .05 .003 .056 .128	.207 .939 .120 .971	.561 .950 .916 .945	.346 .917
	0 1	8 .95 -.043 -.347 .064 .05 .000 -.006 .073 8 .95 -.026 -.306 .064 .05 .003 .028 .075	.100 .888 .075 .970 .098 .911 .074 .972	.036 .903 .943 .948 .067 .942 .926 .944	.001 .946 .004 .933
	2	8 .95 .05	.003 -.155 .064 .000 .032 .075	.068 .947 .072 .967	.317 .939 .928 .951	.056 .918
	0 24 .95 -.007 -.116 .024 .05 .000 -.001 .038 1 24 .95 .000 -.101 .024 .05 .000 .009 .038	.022 .920 .039 .970 .021 .942 .039 .964	.000 .943 .946 .950 .000 .944 .942 .950	.000 .949 .000 .942
	2 24 .95 .05	.001 -.055 .017 .000 .010 .038	.015 .963 .038 .950	.024 .946 .941 .951	.000 .937
		'-' indicates non-existence of the corresponding population moment.

Our bias calculation is related to[START_REF] Alvarez | Robust likelihood estimation of dynamic panel data models[END_REF] in the case of time series homoskedasticity. We calculate the bias of the profile score for the regression coefficients, while they calculated the bias of the conditional score for all common parameters given maximum likelihood estimates of the fixed effects.

For a discussion and historical perspective on this device, see[START_REF] Conniffe | When moments of ratios are ratios of moments[END_REF].[6]

The adjustment, -a(ρ), is robust to cross-sectional heteroskedasticity (σ 2 i instead of σ 2 , say). We investigated this in an older version of this paper[START_REF] Dhaene | An adjusted profile likelihood for non-stationary panel data models with incidental parameters[END_REF], where we also derived the adjustment for the AR(1) model when there are incidental time trends. In the latter case, the adjustment is still a polynomial.

 5 Unbalanced panel data can be accommodated by writing l(θ) as a weighted average of K balanced subpanel profile log-likelihoods, say l k (θ), and adjusting each l k (θ). Missing data, with observations missing completely at random, fit into this scheme by forming complete subseries from each incomplete time series.[8]

When T ≤ 3, b(ρ) is linear and so, given that S(ρ) is concave-convex on [ρ, ρ], b(ρ) ↑ S(ρ) on [ρ, ρ] at most once. Suppose T ≥ 4. Then, b(ρ) is a polynomial of degree 2 or higher with negative coefficients, so b(ρ) is negative, decreasing, and strictly concave, on R + . Further, by Descartes' rule of signs, c(ρ) has one zero on R -when T is even and none when T is odd, and d(ρ) has no zeros on R -when T is even and one when T is odd. Suppose T is even. Then c(-1) = 0 and b(-1) = -1 2(T -1) < 0, so b(ρ) is negative and strictly concave on R, and, hence, its intersection with S(ρ) on [ρ, ρ] can only be on (ρ M L , ρ], where S(ρ) is strictly convex and is approached from below by b(ρ) at most once. Now suppose T is odd and T ≥ 5. Then,

If R is non-empty, to show that b(ρ) ↑ S(ρ) on [ρ, ρ] at most once, it suffices to show that S(ρ) decreases faster than b(ρ) on R, i.e., H(ρ) < c(ρ) for ρ ∈ R. We will show below that (i)

and so, H(ρ) < -1 2 for ρ ∈ [ρ M L , ρ w ]. On the other hand, T (1

So, c(ρ) > -1 2 for ρ ∈ [ρ u , ρ v ] and H(ρ) < c(ρ) for ρ ∈ R. We conclude that b(ρ) ↑ S(ρ) on [ρ, ρ] at most once, provided (i) and (ii) hold, which we now show. Write

2T -4 j=0 v j ρ j 0 , where

2

. We have

Only when 2j + 1 < T is it possible that v 2j + v 2j+1 < 0, so it suffices to show that d k ≥ 0 for 2k + 1 < T . We obtain, for 2k + 1 < T ,

Therefore, (ii) holds.

[29]

Proof of Lemma 2. Use y

where Ξ i = (S 1 ξ 0i , X i ) is independent of ε i . Proceeding as above, we have

where

) is positive semi-definite and V XX is positive definite by assumption. Therefore, V ξξ -V ξX V -1 XX V Xξ ≥ 0 and (V ρρ 0


Geert Dhaene
email: geert.dhaene@econ.kuleuven.be.

K U Leuven

Koen Jochmans
email: koen.jochmans@uclouvain.be.

PROFILE-SCORE ADJUSTMENTS FOR NONLINEAR FIXED-EFFECT MODELS

Keywords: JEL classification: C13, C15, C23, C25 adjusted profile score, bias correction, fixed effects, nonlinear models

Maximum-likelihood estimates of nonlinear panel data models with fixed effects are generally not consistent as the number of units, N , grows large while the number of time periods, T , stays fixed. The inconsistency can be viewed as a consequence of the bias of the score function, where the unit-specific parameters have been profiled out. We investigate ways of adjusting the profile score so as to make it unbiased or approximately unbiased. This leads to estimators, solving an adjusted profile score equation, that are fixed-T consistent or have less asymptotic bias, as T → ∞, than maximum likelihood. One approach to adjusting the profile score is to subtract its bias, evaluated at maximumlikelihood estimates of the fixed effects. When this bias does not depend on the incidental parameters, the adjustment is exact. Otherwise, it does not eliminate the bias entirely but reduces its order (in T ), and it can be iterated, reducing the bias order further. We examine a range of nonlinear models with additive fixed effects. In many of these, an exact bias adjustment of the profile score is possible. In others, suitably adjusted profile scores exhibit much less bias than without the adjustment, even for very small T .

Introduction

Consider the problem of inferring the value of a finite-dimensional parameter θ in a parametric model from a panel data set consisting of T observations on N units. In microeconometric models, unit-specific parameters, called fixed effects or incidental parameters, are often included to account for unobserved heterogeneity. For example, in the agricultural production-function application of [START_REF] Mundlak | Empirical production function free of management bias[END_REF], firm-specific intercepts serve to control for the impact of managerial ability and soil quality on firm output. Alternatively, [START_REF] Hausman | Econometric models for count data with an application to the patents-R& D relationship[END_REF] and [START_REF] Hospido | Modelling heterogeneity and dynamics in the volatility of individual wages[END_REF] introduced fixed effects to allow for heterogeneity in dispersion parameters in applications to the patents-R&D relationship and the volatility of wages, respectively. Unfortunately, including fixed effects generally renders the maximum-likelihood estimator (MLE) of θ inconsistent if T remains fixed while N → ∞ [START_REF] Neyman | Consistent estimates based on partially consistent observations[END_REF]. The problem is known as the incidental-parameter problem.

The profile score function replaces the fixed effects with their maximum-likelihood estimates for a given θ and, therefore, is a feasible version of the score function that would be used if the fixed effects were known. This replacement generally induces a bias of order O(T -1 ). The MLE, θ, sets the profile score to zero and therefore inherits this bias. There are important situations where alternative estimating equations are available that are free of fixed effects (see, e.g., [START_REF] Arellano | Panel data models: Some recent developments[END_REF], for an overview). However, there is no general method for deriving such estimating equations and they may not exist simply because θ may not be fixed-T pointidentified (see [START_REF] Chamberlain | Binary response models for panel data: Identification and information[END_REF]. In this paper we seek to adjust the profile score, following [START_REF] Mccullagh | A simple method for the adjustment of profile likelihoods[END_REF]. The key element is a calculation of the bias of the profile score, either analytically or via simulation, which is then evaluated at maximum-likelihood estimates of the fixed effects. If the bias is free of fixed effects, this leads to an unbiased estimating equation. Otherwise, it results in an estimating equation whose bias is O(T -2 ). We show that it is possible to iterate the adjustment, yielding adjusted profile scores with bias of successively smaller order, O(T -2 ), O(T -3 ), ... Depending on the situation at hand, the adjustments give rise to estimators that are either fixed-T consistent or have a smaller order of bias than the MLE. Our approach fits into the literature on bias-corrected fixed-effect estimation recently surveyed by [START_REF] Arellano | Understanding bias in nonlinear panel models: Some recent developments[END_REF] and inference from integrated likelihoods [START_REF] Lancaster | Orthogonal parameters and panel data[END_REF]Arellano and[START_REF] Arellano | Robust priors in nonlinear panel data models[END_REF], and the parallel developments in the statistics literature (e.g., [START_REF] Li | Efficiency of projected score methods in rectangular array asymptotics[END_REF]Waterman, 2003, and[START_REF] Sartori | Modified profile likelihood in models with stratum nuisance parameters[END_REF].

Focusing on the profile score rather than on θ directly has some advantages. First, it offers a direct way of verifying whether the presence of incidental parameters effectively leads to the inconsistency of the MLE. For example, in the fixed-effect Poisson and exponential-regression models a short calculation suffices to show that the MLE is consistent. On the other hand, verifying whether θ and the incidental parameters are likelihood orthogonal, which is sufficient for the consistency of the MLE, may be a cumbersome task, especially because that may be true in one parametrization but not in another (see [START_REF] Lancaster | The incidental parameter problem since 1948[END_REF]. For the Poisson model, for instance, the equivalence between maximum likelihood and the conditional-likelihood estimator introduced in [START_REF] Hausman | Econometric models for count data with an application to the patents-R& D relationship[END_REF] was not known until [START_REF] Lancaster | Orthogonal parameters and panel data[END_REF] and [START_REF] Blundell | Individual effects and dynamics in count data models[END_REF]. Second, there are several models of practical interest where the bias of the profile score, although non-zero, is free of incidental parameters while the magnitude of plim N →∞ θ depends on the distribution of incidental parameters and covariates. One important model where this is the case is the linear dynamic fixed-effect model [START_REF] Nickell | Biases in dynamic models with fixed effects[END_REF][START_REF] Dhaene | An adjusted profile likelihood for non-stationary panel data models with fixed effects[END_REF]. Weibull and gamma duration models are other examples; details are provided below. In such cases, unbiased estimating equations and fixed-T consistent estimators can be formed by centering of the profile score, a point already made by

[2] [START_REF] Neyman | Consistent estimates based on partially consistent observations[END_REF]. Consequently, adjusting the profile score can be much simpler than approaches targeted at adjusting the MLE directly; see the approaches described in MacKinnon and Smith (1998), for example. In addition, there is an asymptotic justification, complementing the discussion and motivations offered in [START_REF] Mccullagh | A simple method for the adjustment of profile likelihoods[END_REF], for using an adjusted profile score also in situations where the expected profile score does depend on the fixed effects. Fourth, the adjustments, including the iterated adjustment, are easy to carry out. They do not require explicit knowledge of the dependence of the bias on the fixed effects. This is in contrast to the bias correction methods in [START_REF] Arellano | A likelihood-based approximate solution to the incidental parameter problem in dynamic nonlinear models with multiple effects[END_REF], [START_REF] Carro | Estimating dynamic panel data discrete choice models with fixed effects[END_REF], and [START_REF] Bester | A penalty function approach to bias reduction in non-linear panel models with fixed effects[END_REF]. Fifth, the iterative procedure leads to higher-order bias adjustments, as does the jackknife (Dhaene and Jochmans, 2010a).

Section 2 presents the profile score adjustment and how it can be iterated. We discuss examples in Section 3, mostly to nonlinear models. Section 3 illustrates the gains of the adjustments by simulations in the context of static and dynamic binary-choice models.

Adjusting the profile score

We are given a panel data set (y it , x it ) where i = 1, . . . , N and t = 1, . . . , T . Assume independence across i. The conditional density of y it given x it , f (y it |x it ; θ, η i ), is known up to the common parameter θ and the unit-specific parameters η i . Both θ and η i may be vectors. We are interested in estimating θ. Since [START_REF] Neyman | Consistent estimates based on partially consistent observations[END_REF], it is known that the maximum-likelihood estimate (MLE) of θ need not be consistent as N → ∞ with T fixed. One may view the inconsistency as resulting from a biased profile score function. The (normalized) profile loglikelihood and score functions, and their ith contributions, are

l = 1 N i l i , l i = l i (θ) = 1 T t log f (y it |x it ; θ, η i ), s = 1 N i s i , s i = s i (θ) = 1 T t ∇ θ log f (y it |x it ; θ, η i ),
where η i is the MLE of η i for a given θ,

η i = η i (θ) = arg max η i 1 T t log f (y it |x it ; θ, η i ).
Assuming that f is sufficiently regular, the MLE solves s = 0 for θ. Let E = E θ,η i denote the expectation operator at true parameter values, with exogenous variates (and, possibly, initial observations) held fixed at observed values. As is well known, the expected score vanishes at the true value, i.e. E 1 T t ∇ θ log f (y it |x it ; θ, η i ) = 0. However, the profile score replaces η i with η i . Except in special cases, this makes Es i and its aggregate Es nonzero, causing s = 0 to be a biased estimating equation and the MLE to be inconsistent. Under regularity conditions, [3] Es i = O(T -1 ) = Es as T → ∞ and the O(T -1 ) bias of the profile score carries over to the MLE, that is, the asymptotic bias of the MLE, as N → ∞ and T → ∞ sequentially, is O(T -1 ) and can be large.

Our approach centers on a calculation of Es i , either analytically or numerically, and on the estimation of Es i for given θ. Three mutually exclusive cases arise:

(i) Es i = 0;

(ii) Es i = 0 but Es i is free of η i ;

(iii) Es i = 0 and Es i depends on η i .

In case (i), s is unbiased and the MLE is consistent. The interesting point here is that a simple calculation, that of Es i , will reveal so.

In case (ii), the MLE is inconsistent but the adjusted profile score s -Es is unbiased and free of fixed effects. This paves the way for fixed-T consistent estimation. As it turns out, this is the case, surprisingly, in a number of static nonlinear models and in the linear dynamic model.

In case (iii), [START_REF] Mccullagh | A simple method for the adjustment of profile likelihoods[END_REF] proposed using the adjusted profile score s -Es instead of s, where E = E θ, η i is E but with η i replacing η i . 1 The proposal was made in a more general context than the one considered here. McCullagh and Tibshirani discussed many examples, including several with incidental parameters, where the adjusted profile score improves on the profile score. In search of a general justification, they wrote (p. 342) "the centring of the profile log-likelihood function should improve the consistency of the maximizer of the likelihood" and yet, a few lines down, "We have no strong argument for this claim". We provide a large T asymptotic justification.

Let T → ∞. Consider Es i as a function of η i . Under regularity conditions, replacing η i with η i introduces a relative bias of O(T -1 ), i.e. E Es i = ( 1

E(s -Es) = O(T -2 ).
That is, the adjustment removes the first-order bias from s, leaving only bias of order O(T -2 ).

The adjustment can be iterated. The bias E(s -Es) of the adjusted profile score can be appoximated by E(s -Es), again with relative bias O(T -1 ), and subtracted from s -Es to give the second-order adjusted profile score s -2 Es + E Es 1 In addition to the centering step, McCullagh and Tibshirani also considered a rescaling of the adjusted profile score at to restore the information identity. We omit this step as our focus is on getting the estimating equation correctly centered. Similarly, an alternative adjustment would be to scale Es i by the fisher information matrix, as in [START_REF] Firth | Bias reduction of maximum likelihood estimates[END_REF]. [4] with bias

E(s -2 Es + E Es) = O(T -3 ).
Note that, unlike EEs = Es (since Es is a constant), E Es = Es because E Es is E Es (a constant) but with E evaluated at η i instead of η i for all i. The structure of the iterated adjustments is now apparent. Letting E (k) denote the kfold iteration of E, the jth order adjusted profile score is

j k=0 j k (-1) k E (k) s,
with bias O(T -j-1 ), given regularity conditions.

In general, the profile score adjustments, first and higher-order alike, have only an asymptotic justification. Note, however, that in case (ii) they all coincide with s -Es and the adjustment is "exact". Whether there are interesting cases where some jth order adjustment is exact only as of some j ≥ 2 is not known to us. In nonlinear models, the adjustments are generally not exact (though with important exceptions, as we shall see) but only approximations in the sense that they yield approximately unbiased estimating equations, to varying degrees of approximation. Nevertheless, it is hoped that, even when T is small, they yield improvements over the profile score. Whether that is true has to be examined on a case by case base. At the time of writing, our experience with the high-order adjustments is still limited, although we report on some simulations in Section 3.

Implementing the adjustments and solving the adjusted score equations requires evaluating Es i for given θ and η i . Often Es i is not available in closed form, but it can be approximated by the average of R simulations of s i . For large enough R, this average approximates Es i to any desired accuracy, but for the sake of adjusting bias any R suffices, even R = 1. We do not recommend setting R = 1, however, except perhaps in models where evaluating s i is computationally costly. A small R will only inflate the variance of the estimator somewhat. For the higher-order adjustments, which require evaluating terms like E Es i , we suggest using small values of R in all inner expectations, and possibly a larger R in the outermost expectation. Finally, when approximating expectations by simulations, we suggest to keep the basic stream of random numbers used to generate R data sets, which is essentially of dimension R × N × T × dim y it , constant for all values of θ and fixed effects, and for all levels of depth in E (k) .

The adjustments discussed above seek to alter the estimating equation s = 0 to make it unbiased or approximately unbiased. Extensions are possible, perhaps even outside the parametric setting. One variation, in case (iii), is to slightly modify s = 0 so as make it, in essence, a case (ii) problem, where an exact adjustment is feasible. This possibility arises in the two-period logit model, as we discuss in Section 2. More generally, when q = 0 is some other estimating equation that is free of fixed effects and has bias Eq = O(T -1 ), the type of adjustments discussed is [5] possible provided that the expectation E can be evaluated. Another extension is to quantities other than θ, for example

µ = 1 N i µ i , µ i = µ i (θ, η i ),
where µ i (•, •) is a known function, such as a marginal effect. Replacing η i with η i gives

m = 1 N i m i , m i = m i (θ) = µ i (θ, η i ),
with bias O(T -1 ), given regularity conditions. The bias of m i is Em i -µ i and can be appoximated as Em i -m i to give 2m -Em as a first-order adjustment of m, and j k=0 j+1 k+1 (-1) k E (k) m as a jth order adjustment.

Examples

Our examples are models in which the distribution of a scalar y it depends on a vector x it through η i + β x it or η i exp(β x it ). The common parameter consists of β and possibly an additional scale or shape parameter. Details about calculations of Es i are given in the Appendix.

Models where Es

i = 0
There are several models with fixed effects but where there is no incidental-parameter problem.

Example 1 (Poisson counts) Consider Poisson counts y it with mean λ it = η i exp(β x it ) and independence across t given x i1 , ..., x iT . Here f (y it |x it ; β, η i ) = exp(-λ it )λ y it it /y it ! and θ = β. [START_REF] Lancaster | Orthogonal parameters and panel data[END_REF] and [START_REF] Blundell | Individual effects and dynamics in count data models[END_REF] have shown that β and the η i are likelihood orthogonal after a parameter transformation. Alternatively, a calculation shows that Es i = 0. For given β, the MLE of η i is η i = t y it / t exp(β x it ). Letting λ it = η i exp(β x it ), the profile log-likelihood and score for unit i are

l i = T -1 t (-λ it + y it log λ it ) + c = T -1 t y it -log t exp (β x it ) + β x it + c, s i = T -1 t y it -t λ it x it t λ it + x it ,
where (here and later) c is an inessential constant. From Ey it = λ it , it follows that Es i = 0. One may view Es i = 0 as one implication of the conditional moment conditions given in [START_REF] Chamberlain | Feedback in panel data models[END_REF]. The solution to s = 0, that is, the MLE, achieves the semiparametric efficiency bound [START_REF] Hahn | A note on the efficient semiparametric estimation of some exponential panel models[END_REF].

[6]

Example 2 (Exponential durations)

Let y it be exponentially distributed with mean λ -1

it , where λ it = η i exp (β x it ), and independent across t given x i1 , ..., x iT . The density is f (y it |x it ; β, η i ) = λ it exp(-λ it y it ) and θ = β. Here,

η i = [T -1 t y it exp(β x it )] -1 .
Again letting λ it = η i exp(β x it ), we have

l i = T -1 t (log λ it -λ it y it ) = T -1 t -log t (y it exp (β x it )) + β x it + c, s i = T -1 t -t λ it y it x it t λ it y it + x it .
Conditionally on the x it , the λ it y it are i.i.d. unit-exponential variates. Therefore

E t λ it y it x it t λ it y it = t E λ it y it t λ it y it x it = t T -1 E t λ it y it t λ it y it x it = T -1 t x it
and it follows that Es i = 0. The MLE in the exponential regression model with fixed effects and exogenous regressors is fixed-T consistent. [START_REF] Greene | Fixed and random effects in nonlinear models[END_REF] simulations support this, but a proof seems to be new.

Models where Es

i = 0 but Es i is free of η i
Example 3 (Many normal means) This is [START_REF] Neyman | Consistent estimates based on partially consistent observations[END_REF] classic example of the incidental-parameter problem. The problem is to infer θ = σ 2 from independent observations y it ∼ N (η i , σ 2 ). For any given σ 2 , the MLE of η i is y i = T -1 t y it , so

l i = -(2T ) -1 t (log σ 2 + σ -2 (y it -y i ) 2 ), s i = -(2T ) -1 t (σ -2 -σ -4 (y it -y i ) 2 ),
and

Es i = -(2T σ 2 ) -1 . The MLE of σ 2 is (N T ) -1 i,t (y it -y i ) 2 and converges to σ 2 (1 -T -1 ). But since Es i is free of η i , a feasible unbiased estimating equation is s -Es = 0. Its root, (N T -T ) -1 i,t (y it -y i ) 2
, coincides with the outcome of many other approaches; see, e.g., [START_REF] Mccullagh | A simple method for the adjustment of profile likelihoods[END_REF]. A regression version of the model is y it |x i1 , ..., x iT ∼ N (η i + βx it , σ 2 ), where θ now is β and σ 2 . Here, the bias of the profile score for β is zero and for σ 2 it is -(2T σ 2 ) -1 , as before. Again, solving s -Es = 0 yields the standard solution: (i) the MLE of β (which is least-squares with unit-specific de-meaned data) is left unchanged; (ii) a [7] one degree-of-freedom correction is applied to the MLE of σ 2 . Even though the MLE of σ 2 is inconsistent and σ 2 and β are not profile likelihood orthogonal, the inconsistency does not carry over to the MLE of β. This is because the maximizer of l(β, σ 2 ) with respect to β does not depend on σ 2 (though not vice versa).

Example 4 (Dynamic linear regression) Since [START_REF] Nickell | Biases in dynamic models with fixed effects[END_REF], dynamic linear models have become another classic instance of the incidental-parameter problem. Consider the model y it = η i +βy it-1 +ε it with ε it ∼ N (0, σ 2 ) and unrestricted initial observations y i0 . Let y i = (y i1 , ..., y iT ) , y i-= (y i0 , ..., y iT -1 ) , and define the T × T matrix M = I -T -1 ιι where ι is a vector of ones. For given β and σ 2 , the MLE of η i is η i = T -1 ι (y i -βy i-). The profile log-likelihood and the elements of the profile score for unit i are

l i = - 1 2 log σ 2 + T -1 σ -2 (y i -βy i-) M (y i -βy i-) + c, s iβ = -T -1 σ -2 (y i -βy i-) M y i-, s iσ 2 = - 1 2 σ -2 -T -1 σ -4 (y i -βy i-) M (y i -βy i-) .
Using backward substitution it is easy to show that

Es iβ = -T -1 T t=1 (T -t)β t-1 , Es iσ 2 = -(2T σ 2 ) -1 ;
see, e.g., [START_REF] Alvarez | Robust likelihood estimation of dynamic panel data models[END_REF]. [START_REF] Cox | Parameter orthogonality and approximate conditional inference[END_REF] orthogonalization approach leads to essentially the same result; see [START_REF] Lancaster | Orthogonal parameters and panel data[END_REF]. While the adjusted score equation s-Es = 0 is unbiased it typically has more than one root, so the appropriate root has to be selected; see [START_REF] Dhaene | An adjusted profile likelihood for non-stationary panel data models with fixed effects[END_REF]. When the model is extended with exogenous covariates and p lags of y it , Es i is still available in closed form and remains free of η i [START_REF] Dhaene | An adjusted profile likelihood for non-stationary panel data models with fixed effects[END_REF]. Finally, note that s i depends only on the first two moments of the data, so the calculation of Es i is robust to non-normality.

Example 5 (Weibull durations) In this model, y κ it is exponentially distributed with mean λ -1

it , where λ it = η i exp (β x it ), and independent across t given x i1 , ..., x iT . The density of y it is f (y it |x it ; β, κ, η i ) = κy κ-1 it λ it exp (-λ it y κ it ) and θ = (β , κ) . For given β and κ, the MLE of η i is

η i = [T -1 t y κ it exp (β x it )] -1 . With λ it = η i exp(β x it )
, the profile log-likelihood and score for [8] unit i are

l i = T -1 t log κ + (κ -1) log y it + log λ it -λ it y it = T -1 t log κ + (κ -1) log y it -log t (y κ it exp (β x it )) + β x it + c, s iβ = T -1 t -t y κ it λ it x it t y κ it λ it + x it , s iκ = T -1 t κ -1 + log y it -t (log y it ) y κ it λ it t y κ it λ it .
Given the x it , the λ it y κ it are i.i.d., and so Es iβ = 0 by the same argument as in the exponential regression model. A calculation gives Es iκ = (κT ) -1 , free of η i , so the adjusted score s -Es is feasible and unbiased. Although the profile score for β is unbiased, the MLE of β is inconsistent because the profile score for κ is biased and β and κ are not information orthogonal. [START_REF] Lancaster | The incidental parameter problem since 1948[END_REF] showed that an information-orthogonal transformation of η i exists. Integrating the transformed effects from the likelihood using a uniform prior leads to [START_REF] Chamberlain | Heterogeneity, omitted variable bias, and duration dependence[END_REF] marginal-likelihood estimator.

Example 6 (Gamma durations) Here, y it is gamma distributed with shape paramater κ and scale λ -1

it , where λ it = η i exp (β x it ), and independent across t given x i1 , ..., x iT . The density function is f (y it |x it ; β, κ, η i ) = y κ-1 it λ κ it exp (-λ it y it ) /Γ(κ) and θ consists of β and κ. The MLE of η i for given β and κ is

η i = κ[T -1
t y it exp (β x it )] -1 . Letting λ it = η i exp(β x it ) as before, the profile log-likelihood for unit i is

l i = T -1 t -log Γ(κ) + (κ -1) log y it + κ log λ it -y it λ it or, equivalently, T -1 t -log Γ(κ) + (κ -1) log y it + κ log (κT ) -κ -κ log t (y it exp (β x it )) + κβ x it ,
with partial derivatives

s iβ = T -1 t κ -t y it λ it x it t y it λ it + x it , s iκ = T -1 t -ψ(κ) + log (κT ) + log y it -log t (y it exp (β x it )) + β x it ,
where ψ(κ) is the derivative of log Γ(κ). Again, Es iβ = 0, because the y it λ it are independently gamma distributed with scale one and shape κ, given the x it . A calculation shows that Es iκ =

[9] log (κT ) -ψ(κT ), again free of η i . In this model, as in the linear model, there is an incidentalparameter problem only for κ. The solutions of s = 0 and s -Es = 0 differ for κ but coincide for β because, similar to the linear model, the maximizer of l (β, κ) with respect to β does not depend on κ (though not vice versa). Using a similar argument as in the Weibull model, [START_REF] Chamberlain | Heterogeneity, omitted variable bias, and duration dependence[END_REF] derived a fixed-T consistent estimator.

2.3 Models where Es i = 0 and Es i depends on η i

Example 7 (Two-period negbin2 counts) In this model, y it is a negbin2 count with mean λ it = η i exp(β x it ), variance λ it + γ -1 λ 2 it , and there is independence across t given x i1 , ..., x iT . The probability mass function for y it is

f (y it |x it ; β, γ, η i ) = Γ (γ + y it ) Γ (γ) Γ (y it + 1) λ it λ it + γ y it γ λ it + γ γ ,
where γ > 0 is an overdispersion parameter (γ → ∞ yields Poisson counts). Unlike the fixed effects in the negbin1 model of [START_REF] Hausman | Econometric models for count data with an application to the patents-R& D relationship[END_REF], the fixed effects enter the negbin2 model in the standard way, as a means to control for omitted time-invariant covariates; see, e.g., [START_REF] References Allison | Fixed-effects negative binomial regression models[END_REF] and Winkelmann (2008, p. 227-228). The common parameter is θ = (β , γ) . For T = 2, the simulation results in [START_REF] References Allison | Fixed-effects negative binomial regression models[END_REF] suggest that the MLE of β is free of incidental-parameter bias. The analysis below shows that there is incidental-parameter bias for β when T = 2, but that it is very small and can be ignored for practical purposes. For γ the incidental-parameter bias is much larger. For general T , the MLE of η i for given β and γ satisfies

t y it -η i exp (β x it ) γ + η i exp (β x it ) = 0.
This equation is equivalent to a T th order polynomial equation with a unique positive root. The uniqueness follows on rewriting the equation as

T -1 t y it + γ γ + η i exp (β x it ) = 1.
With λ it = η i exp(β x it ), the profile log-likelihood and score for unit i are

l i = T -1 t log Γ (γ + y it ) -log Γ (γ) + γ log γ + y it log λ it -(γ + y it ) log(γ + λ it ) , s iβ = T -1 t γ y it -η i exp (β x it ) γ + η i exp (β x it ) x it , s iγ = T -1 t (ψ (γ + y it ) -ψ (γ) + log γ -log(γ + η i exp (β x it ))
) .

[10]

The expectations of s iβ and s iγ are

Es iβ = ∞ y i1 =0 ... ∞ y iT =0 s iβ t f (y it |x it ; β, γ, η i ), Es iγ = ∞ y i1 =0 ... ∞ y iT =0 s iγ t f (y it |x it ; β, γ, η i ),
but they are difficult to write in a more accessible form. We computed Es iβ and Es iγ for T = 2 and one-dimensional x it . In this case, η i is the largest root of a quadratic equation and the expectations involve only double sums, so they can be evaluated fast and accurately. In general, Es iβ and Es iγ are non-zero and depend on η i . While Es iγ is large, Es iβ is very small. We also computed -(Eh i ) -1 Es i , where h i = ∇ θ s i , as an approximation to the bias of the ( β, γ) , the MLE. When Es iβ is small and, as turns out to be the case, Eh iββ is not too large and Eh iβγ is very small, this approximation to the bias of β is accurate. We computed the approximate bias of β over a range of values of γ and x i1 < x i2 , with η i and β held fixed at one. For each γ, we found that the maximum (approximate) bias occurs as x i1 ↑ 0 and x i2 ↓ 0. Table 1 gives this maximum bias for γ corresponding to moderate to very high levels of overdispersion. Except for very small γ (very large overdispersion), the maximum bias is small. For other values of (b 1 , b 2 ), the bias is typically much smaller than the maximum bias. Table 2 shows Es and E(s -Es) when η i = β = γ = 1 and x i1 = 0, x i1 = log 2, so that the means of y i1 , y i2 are 1, 2 and the variances 2, 6. For β the bias of the profile score is reduced by a factor 4, for γ by a factor 5. Example 8 (Two-period logit) Consider a pair y i = (y i1 , y i2 ) of independent variables y it with mean F (η i + βx it ), where F (z) = (1 + e -z ) -1 is the logistic distribution at z and (x i1 , x i2 ) = (0, 1) (see, e.g., [START_REF] Chamberlain | Analysis of covariance with qualitative data[END_REF]. Here, θ = β, the log-odds ratio, is the parameter [11] of interest. When y i is (0, 0) or [START_REF]η i ∼ N (0, 1[END_REF][START_REF]η i ∼ N (0, 1[END_REF] the MLE of η i for any given β is infinite in absolute value and l i = s i = 0. For the movers, i.e., those units that have y i equal to (0, 1) or (1, 0), the MLE is η i = -β/2. Therefore, the profile loglikelihood and score for unit i are

l i = -d i01 log(1 + e -β/2 ) -d i01 log(1 + e β/2 ), s i = 1 2 d i01 1 + e β/2 - d i10 1 + e -β/2 ,
where d i01 is a binary indicator for y i = (0, 1) and similarly for d i10 . Using

π i01 = Ed i01 = 1 (1 + e η i ) (1 + e -η i -β ) , π i10 = Ed i10 = 1 (1 + e -η i ) (1 + e η i +β ) = π i01 e -β ,
it follows that

Es i = 1 2 π i01 1 + e β/2 - π i01 e -β 1 + e -β/2 = 1 2 1 -e -β/2 1 + e β/2 π i01 .
Hence Es i depends on η i because π i01 does. Therefore, s-Es is not unbiased (see also [START_REF] Mccullagh | A simple method for the adjustment of profile likelihoods[END_REF]. However, a slight modification to s before addressing its bias essentially leads to case (ii), and to the conditional maximum-likelihood estimator [START_REF] Andersen | Asymptotic properties of conditional maximum-likelihood estimators[END_REF]. Write s as s = 1 2N

N 01 1 + e β/2 -N 10 1 + e -β/2 , where N 01 = i d i01 and N 10 is defined similarly. Now consider q = 1 2(N 01 + N 10 )

N 01 1 + e β/2 - N 10 1 + e -β/2
instead of s. Clearly, s = 0 and q = 0 have the same root-the MLE of β. If π 01 = lim N →∞ N -1 i π i01 exists, then π 10 = lim N →∞ N -1 i π i10 = π 01 e -β and plim N →∞ q = 1 2(π 01 + π 10 )

π 01 1 + e β/2 - π 10 1 + e -β/2 = 1 -e -β/2 2(1 + e -β ) (1 + e β/2 ) = q ∞ (say),
where q ∞ is free of the sequence of η i 's. Therefore, q -q ∞ = 0 is a fixed-T unbiased estimating equation. Its solution coincides with the conditional maximum-likelihood estimator, which is known to be efficient [START_REF] Hahn | A note on the efficient semiparametric estimation of some exponential panel models[END_REF]. While in this model there is nothing new in using q -q ∞ = 0, the illustration shows that normalizing the score function by the number of movers (or, more generally, by the number of informative units) can be helpful in models where conditioning is not possible.

[12]

The upper panel in Figure 1 plots the profile score and various adjustments to it for the twoperiod logit model from Example 8. The plot was generated with N = 1, 000, 000, θ = 1, and η i ∼ N (0, 1). The adjusted profile scores, shown up to third order, were obtained by means of one draw at each level of depth. The plot verifies the well-known results that, for this design, plim N →∞ θ = 2θ. The root of the first-order adjusted profile score is already much closer to the true parameter value. The higher-order corrections further reduce the inconsistency and the adjusted profile scores can be observed to converge to the score of the conditional-likelihood estimator as the order of correction increases. Also plotted is the (first-order) adjusted profile score using only units for which y i1 + y i2 = 1, that is, units that contribute to the profile likelihood. Its root is the true parameter value. As higher-order corrections would leave the location of this curve unchanged, these are not plotted. The bottom panel in Figure 1 contains the scores for the probit variant of the two-period model. It may again be observed that a first-order adjustment greatly reduces the inconsistency of the MLE, and that iterating the corrections further centers the profile score. While, again, using only movers improves the situation, it does not lead to an adjusted profile score whose root is at the correct parameter value, although the inconsistency is small.

Figure 2 provides the profile score and the various adjustments for dynamic binary-choice models. They were obtained in an analogous fashion to before, only now with T = 3, which is the shortest panel length for which the MLE is finite. The start-up values for the N time series were drawn from their respective stationary distributions. The plots verify that dynamics tend to increase the magnitude of the inconsistency of the MLE. The re-centering effect of the score corrections is similar to before. The analogy of movers in the static model are units that alternate. More precisly, only units that switch status in each time period contribute to the profile likelihood. Here, using only such sequences does not lead to consistency of the root of the adjusted profile score in the logit model.

A small Monte Carlo experiment was performed to evaluate how much bias can be eliminated in small samples. The models considered are logit and probit variants of the binary-choice model

y it = 1(η i + x it θ ≥ ε it ),
with both η i and x it scalar i.i.d. standard-normal variates. N was set to 100 and θ was fixed at unity throughout. Tables 3-4 contain the mean and median bias of the MLE ( θ) and the root of the first-and second-order corrected profile scores ( θ 1 and θ 2 ), along with their standard deviation (STD) and interquartile range (IQR). The bias of the profile score was computed through simulation, with 10 runs for the outermost iteration step, and a single run in the inner iteration. For logit a useful benchmark is the conditional-likelihood estimator ( θ), and so results for this estimator are also included. For probit, no fixed-T estimator is available.

[13] Score and adjusted score functions for the two-period logit (upper panel) and probit (lower panel) model with a time trend. N = 1, 000, 000; η i ∼ N (0, 1). In both panels, s(θ) (solid) is plotted along with its first-(dashed), second-(dashed; marked +), and third-order (dashed; marked * ) adjustments, together with their scaled counterparts (dashed-dotted). For the logit model, the score of the conditional likelihood (dotted) is also plotted.

[14] Score and adjusted score functions for the three-period AR(1) logit (upper panel) and probit (lower panel) model. In both panels, s(θ) (solid) is plotted along with its first-(dashed), second-(dashed; marked +), and third-order (dashed; marked * ) adjustments, together with their scaled counterparts (dashed-dotted).

The tables show that the MLE suffers from a substantial upward bias. The results also suggest that our large-sample arguments to bias correction tend to give a reasonable approximation. Moreover, solving adjusted profile scores yields estimates with much smaller bias than the MLE. Iterating the correction further reduces the bias. It is apparent from the results that the bias is virtually fully eliminated when T > 2. Interestingly, the bias-corrected estimators are also less variable than is the MLE. It is known from [START_REF] Neyman | Consistent estimates based on partially consistent observations[END_REF] that, with incidental parameters, the MLE need not be asymptotically efficient, even if it is fixed-T consistent. Notice, finally, that bias-corrected estimation of the logit model does not perform better than conditionallikelihood estimation.

Appendix

Weibull durations Write log y it as log y it = κ -1 log(y k it λ it ) -log λ it = κ -1 (log e t -log λ it ) , where the e t are i.i.d. unit exponentials given the x it . Then we have Gamma durations Here, write log y it as log y it = log (y it λ it ) -log η i -β x it = log g t -log η i -β x it where, given the x it , the g t are i.i.d. gamma variates with shape κ and scale one. Then, write s iκ as s iκ = T -1 t -ψ(κ) + log (κT ) + log g t -log t g t .

s iκ = T -1 t κ -1 + κ -1 (log e t -log λ it ) -κ -1 t (
Using E log g t = ψ(κ) and the property that t g t is gamma distributed with shape κT and scale one, it follows that Es iκ = log (κT ) -ψ(κT ).

[17]

  + O(T -1 )) Es i and, on averaging over i, E Es = (1 + O(T -1 )) Es. Therefore, moving from the profile score s to McCullagh and Tibshirani's adjusted profile score s -Es reduces the bias from Es = O(T -1 ) to

Figure 2 :

 2 Figure 2: Plots for dynamic three-period binary-choice models

  γ is Euler's gamma. We used Mathematica to calculate the integral. Setting T = 1 gives E log e t = -γ. On collecting results, Es iκ = (κT ) -1 .

Table 1

 1

		: Maximum approximate bias of β
	γ	.5	1	2	5	10	20
	maxbias( β) .1108 .0532 .0211 .0046 .0012 .0003

Table 2 :

 2 Bias before and after adjustment

		β	γ
	Es	.0021 .133
	E(s -Es) .0005 .025

  Figure 1: Plots for static two-period binary-choice models

						Logit
	0.12
	0.1
	0.08
	0.06
	0.04
	0.02
	0
	-0.02
	-0.04
	-0.06
	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
						!
						Probit
	0.25
	0.2
	0.15
	0.1
	0.05
	0
	-0.05
	-0.1
		0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
						!

  log e t -log λ it ) e t t e

t , with expectation Es iκ = κ -1 1 + E log e t -T E e t log e t t e t . The sum in the denominator is e t +A where A is independent of e t and has the Erlang distribution with density A T -2 exp(-A)/ (T -2)!. Therefore, E e t log e t t e t