Amemiya A Ronald Takeshi

John F Gallant

Cheng Geweke

Peter Hsiao

Arnold Robinson

Zellner

Aït-Sahalia Yacine

H Badi

M W Baltagi

Marcus J Brandt

Songnian Chambers

Manfred Chen

Miguel A Deistler

Jean-Marie Delgado

Sylvia Dufour

Eric Fruhwirth- -Schnatter

John C Ghysels

Javier Ham

Han Hidalgo

Yongmiao Hong

Bo E Hong

Maxwell L Honoré

Yuichi King

G M Kitamura

Chung-Ming Koop

Naoto Kuan

Kajal Kunitomo

Q Lahiri

Tong Li

Oliver Li

James G Linton

Robert Mackinnon

Rosa L Mcculloch

Franz C Matzkin

Dale J Palm

Nicholas Poirier

B M Polson

Ingmar Pötscher

Peter C Prucha

Eric Reiss

Frank Renault

Robin Schorfheide

Fallaw Sickles

G J Sowell

Herman Van Den Berg

Dijk Van

H Quang

Edward Vuong

Tom Vytlacil

Andrew Weiss Wansbeek

Tao

Jan Siem

André Koopman

André Lucas

Monteiro

Koopman

André Monteiro

Jaap Abbring

Luc Bauwens

Arnaud De Servigny

Pieter Klaassen

Michel Mouchart

Olivier Renault

The multi-state latent factor intensity model for credit rating transitions

Keywords: unobserved components, credit cycles, duration model, generator matrix, Monte Carlo likelihood. JEL classification codes: C15, C33, C41, C43, G11, G21

A new empirical reduced-form model for credit rating transitions is introduced. It is a parametric intensity-based duration model with multiple states and driven by exogenous covariates and latent dynamic factors. The model has a generalized semi-Markov structure designed to accommodate many of the stylized facts of credit rating migrations. Parameter estimation is based on Monte Carlo maximum likelihood methods for which the details are discussed in this paper. A simulation experiment is carried out to show the effectiveness of the estimation procedure. An empirical application is presented for transitions in a 7 grade rating system. The model includes a common dynamic component that can be interpreted as the credit cycle. Asymmetric effects of this cycle across rating grades and additional semi-Markov dynamics are found to be statistically significant. Finally, we investigate whether the common factor model suffices to capture systematic risk in rating transition data by introducing multiple factors in the model.

Ratings play a prominent role in the credit industry. Their key purpose is to provide a simple qualitative classification of the solidity, solvency and prospects of a debt issuer. The importance of credit ratings has increased significantly with the introduction of the new regulatory framework known as Basel II (BCBS, 2004). In this framework, ratings can be used directly to determine the size of a bank's capital buffer. As capital constitutes a relatively costly source of funding for a bank, ratings and rating changes directly affect the banks' willingness to grant credit to individual firms. Moreover, if ratings and thus capital requirements co-vary with the business cycle, economic fluctuations may be exacerbated by capital becoming increasingly scarce in adverse economic conditions, precisely when it is needed most. It is clear that a good understanding of the dynamic behavior of ratings and rating changes is therefore important from both a regulatory and financial industry perspective.

In this paper we introduce a new model for rating transitions. The main novelty of our model is that rating transitions are modeled continuously in event time rather than calendar time and are subject to common dynamic latent factors. Although the model is relatively complex, we

show that it can be estimated efficiently using modern importance sampling techniques for non-Gaussian models in state space form.

The literature on modeling credit events such as defaults and rating changes has grown rapidly over the past 10 years. Wilson (1997a,b) modeled default rates using logistic regressions with macroeconomic explanatory variables. [START_REF] Nickell | Stability of rating transitions[END_REF] and [START_REF] Bangia | Ratings migration and the business cycle, cith applications to credit portfolio stress testing[END_REF] show that upgrade, downgrade, and default probabilities differ over different economic regimes, whether characterized by NBER business cycle classifications or by GDP growth rates. Default and downgrade intensities are higher during recessions. In the same spirit, [START_REF] Kavvathas | Estimating credit rating transition probabilities for corporate bonds[END_REF], [START_REF] Carling | Capital charges under Basel II: Corporate credit risk modeling and the macro economy[END_REF], [START_REF] Couderc | Time-to-default: Life cycle, global and industry cycle impacts[END_REF], and Duffie, Saita, and Wang (2006a) use a duration approach conditional on observed macroeconomic and firm characteristics and show that average times-to-default decrease if economic activity decreases. Koopman and Lucas (2005) and [START_REF] Koopman | Empirical credit cycles and capital buffer formation[END_REF] have adopted a direct time series approach and identified the time-varying cyclical nature of default rates over a long historical period. Also Fledelius, [START_REF] Fledelius | Non-parametric analysis of rating transition and default data[END_REF] corroborate the existence of time-fluctuations for credit rating migration rates.

Whereas some of the contributions in the literature introduce observed macro-variables to capture co-variation in default intensities between firms and industries, an alternative approach is to estimate the common components of default risk directly from the data. An advantage of such an approach is that one is less prone to misspecification caused by the use of an incorrect A c c e p t e d m a n u s c r i p t macroeconomic proxy for the credit cycle. [START_REF] Couderc | Time-to-default: Life cycle, global and industry cycle impacts[END_REF] tested a large number of macroeconomic variables for their predictive ability and found five significant factors. Still, a large part of the fluctuations in systematic default probabilities could not be accounted for.

Second, by estimating the default dynamics directly from the data, one obtains an integrated framework for capital determination and risk management, see [START_REF] Koopman | Empirical credit cycles and capital buffer formation[END_REF]. By contrast, if observed macroeconomic variables are used, one needs an auxiliary forecasting model for such variables, see for example Duffie et al. (2006a,b).

Suggestions for dynamic models with latent components are [START_REF] Gagliardini | Stochastic migration models with application to corporate risk[END_REF], [START_REF] Mcneil | Bayesian Inference for Generalized Linear Mixed Models of Portfolio Credit Risk[END_REF], and [START_REF] Koopman | A non-Gaussian panel time series model for estimating and decomposing default risk[END_REF]. These models, however, are all set in a calendar time framework: rating transitions are observed empirically over discrete time slots, e.g., years or quarters. The observed frequencies are subsequently modeled by non-Gaussian time series processes. By contrast, in this paper we use a duration model with unobserved components. The duration (continuous time) approach is the more natural approach in the current context, where durations to transitions are endogenous rather than exogenous. In this way, we are able to use all the information in the data-set [START_REF] Lando | Analyzing rating transitions and rating drift with continuous observations[END_REF] provide a further detailed discussion of the advantages of the continuous-time approach. Our model can be regarded as a multi-state extension of the Latent Factor Intensity (LFI) model of [START_REF] Bauwens | Dynamic latent factor models for intensity processes[END_REF]. The LFI model is a point process model for stock transactions in tick-time. Durations in the LFI model are the time to the next trade. By contrast, in our model it is not only the time to the next rating event that is unknown, but also the type of event that is going to occur, e.g., upgrade, downgrade, or default. In that sense, our model is set in the so-called competing risks framework. Given a firm's initial rating, there are multiple states for the firm's next rating. Each of these states has its own duration process and we observe only the minimum of those. This leads to a more complicated likelihood structure than considered by [START_REF] Bauwens | Dynamic latent factor models for intensity processes[END_REF].

The likelihood function of our model contains a high dimensional integral involving the latent common risk factor. In this way, our parameter driven model differs from well-known observation driven counterparts like the Autoregressive Conditional Duration model (ACD) of [START_REF] Engle | Autoregressive conditional duration: a new model for irregularly spaced transaction data[END_REF], or the Autoregressive Conditional Intensity model (ACI) of [START_REF] Russell | Econometric Modelling of Multivariate Irregularly spaced Highfrequency data[END_REF]. We evaluate the likelihood using a multivariate extension of the Monte Carlo techniques that are developed by Durbin andKoopman (1997, 2001). We demonstrate the effectiveness of the method by means of a simulation experiment.

The model is estimated for the CreditPro7.0 data set from Standard & Poor's, containing all issuer ratings over the period 1981 -2005. We classify firms into 7 standard rating categories and specify a dynamic model for upgrades, downgrades, and defaults using all available

A c c e p t e d m a n u s c r i p t

data. This yields a data set including almost 7000 firms and almost 14000 informative rating events (more than 25000 if we also count sample extension and attrition). We obtain some interesting empirical findings. First, there is significant evidence of a persistent common component in rating transitions. We further show that the impact of this common component with respect to downgrade and upgrade probabilities is asymmetric. Upgrades are idiosyncratic to a large extent, whereas downgrades and defaults tend to cluster together in time. Further, experiments with multiple (latent) factor models cast doubt on the suitability of (widespread) single factor models for risk management and capital buffer determination. Finally, we find statistically significant evidence of semi-Markov effects in transitions and defaults. Investment grade transition intensities tend to increase over the first few years and become more stable thereafter. Sub-investment grade companies, by contrast, show increased transition activity at short durations, probably due to momentum effects, and at longer durations, possibly due to debt roll-over.

The paper is organized as follows. In Section 2, the model is presented. In Section 3 we develop the estimation methodology for this model. Section 4 discusses how to obtain default probabilities over finite time periods from the event time specification. Section 5 contains the results of a Monte Carlo study. Section 6 presents the results of our empirical study. We conclude in Section 7.

The Multi-State Latent Factor Intensity model

The multi-state latent factor intensity (MLFI) model is a multi-state generalization for multivariate point processes of the latent factor intensity (LFI) model of [START_REF] Bauwens | Dynamic latent factor models for intensity processes[END_REF]. Consider a set of K units (or firms) whose event-histories can be adequately described by the history of transitions between a finite set of states. The states in our empirical application will be the set of credit ratings for issuers as assigned by Standard and Poor's (S&P).

The data set has a clear panel structure and consists of the exact dates and the corresponding type of the rating changes recorded for each firm in the sample. In order to account for unobserved dependence between the transition histories in a parsimonious way, we introduce a common factor ψ(t). We assume that conditional on ψ(t), rating events are independent across firms (i.e., along the cross section dimension). This assumption is standard in the credit risk literature and is used to prevent the model's corresponding joint state-space becoming quickly unmanageable due to its size. [START_REF] Gagliardini | Stochastic migration models with application to corporate risk[END_REF] processes N k (t) and N (t). The process N (t) makes a jump of unit size at each time there is a rating event for one of the K units.1 Similarly, N k (t) jumps at the times there is a credit event for unit k such that

N (t) = K k=1 N k (t).
These point processes are marked because at each event time we also observe the transition type of the unit, i.e., the specific type of upgrade or downgrade. In fact, the counting process N k (t) can be expressed as the sum of S counting processes N sk (t) that keep track of the total number of transitions of type s for firm k. It follows that

N k (t) = S s=1 N sk (t), N (t) = K k=1 N k (t) = S s=1 K k=1 N sk (t).
Corresponding to each of these point processes N sk (t) we assume there is a finite stochastic intensity λ sk (t). In practical terms this intensity describes the instantaneous probability of unit k experiencing a type s rating transition at time t conditional on the information available just before time t. Naturally, such transition intensities are only defined at time t if the unit actually is 'at risk' for transition type s at t -< t, where t -t -is arbitrarily small. For example, the downgrade intensity from AAA to AA for firm k at time t is only defined if firm k actually has an AAA rating just prior to t. The intensity for each point process2 λ sk (t) can be (informally) defined by

λ sk (t) = lim ∆↓0 P [N sk ((t + ∆) -) -N sk (t -) > 0 | F t -] ∆ ,
see for example Andersen et al. (1993, p. 51). The conditional information up to (but not including) time t is represented by F t -= ∪ τ <t F τ for an appropriate filtration F τ .

Define R sk (t) as a dummy variable that takes the value one if unit k is 'at risk' for transition type s ∈ S at time t -, and zero otherwise. Note that unit k can be at risk for multiple transition types at the same time. For example, both the AAA to AA and the AAA to A transitions may be at risk simultaneously. Obvious reasons for a transition type not to be at risk for firm k at We adopt a proportional hazards specification. The model specification for the conditional hazards is given by

λ sk (t) = exp [η s + γ ′ s w k (t) + α s ψ(t)] • H sk (t), (1)
with s = 1, . . . , S and k = 1, . . . , K, where (i) scalar η s , m × 1 vector γ s , and scalar α s are fixed unknown coefficients, (ii) m × 1 vector w k (t) contains explanatory variables (covariates), (iii) scalar ψ(t) represents a latent dynamic factor, and (iv) scalar function H sk (t) represents the generalized baseline hazard function, which can be used to model duration dependence of the multivariate type. This specification encompasses for example the homogeneous continuoustime Markov chain model that is frequently used in the empirical credit risk literature, see, for example, [START_REF] Kavvathas | Estimating credit rating transition probabilities for corporate bonds[END_REF] and [START_REF] Lando | Analyzing rating transitions and rating drift with continuous observations[END_REF]. Note we assume λ sk (t) to be only defined if the corresponding R sk (t) = 1.

A more detailed discussion of the intensity specification (1) follows below. The parameter η s represents the reference-level log-intensity of transition type s. It is independent of time and common across all units k = 1, . . . , K. The parameter vector γ s and scalar α s measure the sensitivity of unit k's log-intensity for transition type s with respect to observed covariates w k (t) and the unobserved process ψ(t), respectively. The m-dimensional vector of covariates w k (t) can contain unit-specific information such as leverage and profitability ratios, industry dummies, stock volatilities or statistics depending on the rating process.3 Further, w k (t) can include macroeconomic information such as economic growth rates, interest rate levels and term structure variables. In this case subscript k can be dropped from the notation. Note that phenomena like rating momentum can also be included in w k (t) such that past downgrades and upgrades make subsequent downgrades and upgrades more likely, respectively.

The coefficients α s depend on the transition type s ∈ S. This implies that α s can depend on both the origin and the destination state. In the empirical literature it is common practice to let α s parameters depend on the origin state, i.e., the initial rating, only. Here, however, the impact of the common risk factor ψ(t) depends on the type of transition, and therefore on the destination state as well. For example, upgrades might be less subject to common risk factors than downgrades, see [START_REF] Gagliardini | Migration correlation: Definition and efficient estimation[END_REF]. Restrictions on α s can be tested explicitly using the maximum likelihood based procedure of Section 3.

Following the empirical work in the credit risk literature, we assume all intensities are subject

A c c e p t e d m a n u s c r i p t

to the same unobserved dynamic common factor ψ(t). Relaxing this assumption by making ψ(t), for example, rating or industry specific is conceptually straightforward in our modeling framework. The latent process might even be unit specific as in [START_REF] Bauwens | Dynamic latent factor models for intensity processes[END_REF].

In the case of rating transition data, however, specifying unit-specific processes is not really feasible. The number of rating events for an individual firm is usually too small, even over a prolonged period of time. This is a direct consequence of the rating agencies' policy to provide stable ratings to the investment community.

Assume that ψ(t) only changes at observed event times t i for i = 1, . . . , N (T ) where T denotes the time index of the last observation (right-censoring of type I). The specification of ψ(t) as a stochastic process with piecewise constant (left-continuous) sample paths is intuitive since the intensity of the pooled process (pooled over firms and transition types) is not identified between two consecutive events. Moreover, in the context of credit rating transitions, ψ(t) is intended to capture low-frequency co-movements in the vector of migration intensities. In the empirical illustration of Section 6, the average duration of the pooled process is 1.2 days.

Therefore, no serious bias will arise from disregarding possible changes in the macroeconomic variables over the almost bi-daily spells of the pooled process.

Let ψ i = ψ(t i ) denote the value of the common risk factor ψ(t) over the interval t ∈ (t i-1 , t i ].

In order to capture serial correlation in the intensity of the pooled process, the dynamic process for ψ i can be specified, for example, by a first order autoregressive (AR) equation

ψ i = ρ t i -t i-1 s ψ i-1 + σ i ε i , (2)
where ε i is a set of i.i.d. N(0, 1) innovations, the AR parameter ρ ∈ [-1, 1], and where σ i may depend on the length of the spell (t i -t i-1 ). The specification in (2) has small changes in ψ(t) over short spells. This is in line with our interpretation of ψ(t) as an economy wide risk factor, which one would not expect to vary wildly at high frequencies. As not all α s parameters from

(1) and σ i can be identified simultaneously, we normalize σ i to

σ 2 i =    (1 -ρ 2(t i -t i-1 )/D )/(1 -ρ 2/D ) for -1 < ρ < 1, t i -t i-1 for ρ = 1, (3)
in our empirical work, with t i measured in days of a D = 260 day business year. This brings ψ(t)

close to an Ornstein-Uhlenbeck (for |ρ| < 1) or Brownian Motion (for ρ = 1) process observed at the event days t i . More general dynamic specifications for ψ i can be easily incorporated in the state space framework of the next subsection. For example, different AR processes can be considered for specific rating transitions. See also the empirical application in Section 6.

The baseline hazard H sk (t) is specified by the deterministic function The inclusion of H sk (t) introduces duration dependence into the model and, therefore, relaxes the Markov assumption. More precisely, if H sk (t) is allowed to depend only on t -t N k (t),k , then each unit follows a semi-Markov process. In the general case a generalized semi-Markov process is obtained, see [START_REF] Glynn | A GSMP formalism for discrete-event systems[END_REF]. Possible choices for H sk (t) include the hazard function of a multivariate Weibull distribution, given by

H sk (t) = H s (t -t 0k , t -t 1k , . . . , t -t N k (t),k ), (4)
H s (x 0 , . . . , x N ) = N i=0 a si x bs-1 i , (5)
with x i ≥ 0 and fixed coefficients a si > 0 and b s > 0 for i = 0, 1, . . . , N . Another valid alternative is the self-excitation mechanism introduced by [START_REF] Hawkes | Spectra of some self-exciting and mutually exciting point processes[END_REF] and also considered for the LFI model by [START_REF] Bauwens | Dynamic latent factor models for intensity processes[END_REF].

We note that k's observed duration or spell

t N k (t),k -t N k (t)-1,k is the minimum of s R sk (t)
latent durations corresponding to the set of feasible transitions 'at risk' for unit k at time t.

We adopt the standard practice of assuming that the latent duration processes are mutually independent conditional on the common factor ψ(t).4 See van den Berg ( 2001) for a detailed discussion on identification problems in this setting.

To complete the model specification, an additional set of identifying assumptions for the parameters is required. The global identification of intensity specification (1) requires a sign restriction for α s . Changing the sign simultaneously for all α s 's and for the complete path of ψ(t) clearly yields the same path for intensity λ sk (t). We therefore set α s < 0 for s = S.

For a vector of unknown parameters θ, the likelihood function conditional on the initial ratings, pre-sample event histories,5 and on the complete path of the unobserved process, as defined by Ψ N (T ) = {ψ i } N (T ) i=0 , can be written as

L θ | F T , Ψ N (T ) = N (T ) i=1 K k=1 S s=1 exp Y sk (t i ) ln{λ sk (t i )} -R sk (t i ) t i t i-1 λ sk (t)dt , (6)
where dummy variable Y sk (t) is one if unit k at time t experiences a rating event of type s, and zero otherwise, and F T denotes the relevant observable filtration. The likelihood function (6) has an intuitive interpretation. Unit k only contributes to the (conditional) likelihood if it is at risk, that is if R sk (t i ) = 1. In this case, the likelihood contains the probability of survival of unit k in its current state over each spell of the pooled point process if there was no rating A c c e p t e d m a n u s c r i p t event for this unit at risk. When rating event i takes place at the end of the spell of the pooled process for firm k, that is if Y sk (t i ) = 1, the survival probability is multiplied by the hazard rate to yield the probability density of the rating event.

The likelihood in (6) can be decomposed in a likelihood of a spell length in the pooled point process, and a likelihood of the associated mark. The spell length follows from the cumulated (over all s and k) intensities at risk, R sk (t)λ sk (t). The mark then follows from the multinomial distribution with probabilities R sk (t)λ sk (t)/ s,k R sk (t)λ sk (t) for firm k experiencing rating event type s. See also Section 5. This is the most general specification. 6 More restrictive specifications can of course be accommodated as well. For example, one might take the ordered nature of ratings into account by restricting the α s parameters to depend on the initial/input rating only, and not on the output rating. In this paper, however, we stick to the more general specification.

In order to estimate the parameter vector θ, the conditional likelihood function must be integrated with respect to the complete path Ψ N (T ) of the unobserved process ψ(t). The maximum likelihood problem becomes

max θ L(θ | F T ), (7)
where

L(θ | F T ) = L θ | F T , Ψ N (T ) p(Ψ N (T ) )dΨ N (T ) , (8)
and p(Ψ N (T ) ) denotes the density function of Ψ N (T ) .

Monte Carlo Maximum Likelihood Estimation

The main difficulty with maximum likelihood estimation in (8) is the computation of the highdimensional integral. In a typical application such as the one in the next section, this integral is much more than 4000 dimensional. [START_REF] Mcneil | Bayesian Inference for Generalized Linear Mixed Models of Portfolio Credit Risk[END_REF] address a similar problem by adopting a Bayesian perspective, albeit in a lower dimensional space (around 50). [START_REF] Bauwens | Dynamic latent factor models for intensity processes[END_REF] adopt the simulated maximum likelihood method of [START_REF] Liesenfeld | Univariate and multivariate stochastic volatility models: Estimation and diagnostics[END_REF] for the estimation of a single-state LFI model.

6 Given the close resemblence of the process for the marks with qualitative response models, well-known problems for those classes of models like (in)dependence of irrelevant alternatives (IIR), will be an issue here as well, at least for the general multinomial specification. Though not the issue of this paper, the current methodology is also applicable in more general classes of point process models that are less subject to such problems. Moreover, in our specific application to credit risk modeling within a 7 grade ratings system, IIR is not a key issue.

A c c e p t e d m a n u s c r i p t

By contrast, in this paper, we adopt the general method of Monte Carlo maximum likelihood for a multi-state LFI (MLFI) model. To overcome the inefficiency problem of direct Monte Carlo estimation of the high-dimensional integral in (8) we use a combination of importance sampling and the Kalman filter and smoother as described in Durbin and Koopman (Part II, 2001). It is shown that the methodology can be made applicable for high-dimensional problems. In this section, the model is formulated in state space form in Subsection 3.1. The Monte Carlo simulation method for likelihood evaluation is discussed in Subsection 3.2.

Statistical model specification

The MLFI model considers the following three sources of stochastic variation: (i) the duration between events in the pooled process, denoted by τ i = t i -t i-1 ; (ii) the transition types s being at risk at t - i for unit k, denoted by R sk (t i ); (iii) the specific transition type s at time t i for unit k, denoted by Y sk (t i ). These stochastic variables are collected in the vector z i for i = 1, . . . , N with N = N (T ), where z i is defined as

z i = {τ i , R 11 (t i ) , . . . , R SK (t i ) , Y 11 (t i ) , . . . , Y SK (t i )} ′ .
The vector z i can be constructed (or observed) at each event i = 1, . . . , N . The analogue of the observation equation for z i is implied by the non-Gaussian conditional likelihood in (6). In particular, for the ith event time of the pooled process, we have the conditional log-density

ln p(z i |ψ i , F t - i ) = S s=1 K k=1 Y sk (t i ) ln{λ sk (t i )} -R sk (t i ) t i t i-1 λ sk (t)dt, (9)
for i = 1, . . . , N .

The intensity specification (1) can be formulated more generally via vector ν i that contains latent processes and fixed effects. We have

λ sk (t) = exp (Z ski ν i ) • H sk (t), for t i-1 < t ≤ t i , (10)
where Z ski is a fixed and known 'selection' vector, for s = 1, . . . , S, k = 1, . . . , K, and i = 1, . . . , N . In case

H sk (t) = 1, intensity λ sk (t) is constant for t i-1 < t ≤ t i .
To show that specifications (1) and ( 10) can be equivalent, we take

ν i = {η 1 , . . . , η S , γ ′ 1 , . . . , γ ′ S , ψ(t i )} ′ , Z ski = {e ′ s , e ′ s ⊗ w k (t i ) ′ , α s } ,
where e s is the s-th column of I S . It follows that

Z ski ν i = η s + γ ′ s w k (t i ) + α s ψ i .
If another specification for λ sk (t) is considered, the specifications for Z ski and ν i need to be adjusted accordingly.

A c c e p t e d m a n u s c r i p t

The vector ν i can contain both fixed unknown coefficients and dynamic latent processes.

We therefore model ν i by the general Markovian process

ν i = F i ν i-1 + Ri η i , η i ∼ NIID(0, Q i ), i = 1, . . . , N, (11)
with initial condition ν 0 ∼ N(a, P ). The vector a and the matrices F i , Ri , Q i and P are fixed matrices that may depend on the parameter vector θ. If the vector ν i only consists of fixed unknown coefficients, we set a = 0, F i = Ri = I, Q i = 0 and P = κI, where κ is the so-called diffuse prior constant. Usually, κ is set to some large value in numerical software, see Harvey (1989, pp. 367-8). Exact solutions for κ → ∞ are available as well, see Durbin and Koopman (2001, Ch. 4). If the vector ν i only contains the latent autoregressive process [START_REF] Gagliardini | Stochastic migration models with application to corporate risk[END_REF], that is

ν i = ψ i , we set a = 0, F i = ρ t i -t i-1 , Ri = 1, Q i = σ 2 i and P = (1 -ρ 2 ) -1 .
A combination of unknown coefficients and latent time series processes can be incorporated in (11) in a straightforward way. For example, in the case of (1) with w k (t) = 0, we have

ν i = (η 1 , . . . , η S , ψ i ) ′ with a = 0, F i =   I S 0 0 ρ t i -t i-1   , Ri =   0 1   , Q i = σ 2 i , P =   κI S 0 0 (1 -ρ 2 ) -1   .
In case multiple latent factors or higher order dynamics of the latent factors are part of the model, the state vector ν i can be extended in a natural way. The state space matrices need to be adjusted accordingly. A general framework for the MLFI model can be summarized by the observation log-density for z i conditional on the state ν i . This is given by ( 9) where λ sk (t) is

given by (10) for t i-1 < t ≤ t i , and where ν i is modeled by ( 11) with i = 1, . . . , N . This set of equations makes up a nonlinear non-Gaussian state space model as considered by [START_REF] Shephard | Likelihood analysis of non-Gaussian measurement time series[END_REF] and Durbin andKoopman (1997, 2001). Importance sampling methods enable the evaluation of the likelihood function but also the estimation of ν i and the computation of the corresponding standard errors for i = 1, . . . , N . The details of these methods are given in Subsection 3.2.

We notice the high dimension of the state vector ν i since it includes the scalars η 1 , . . . , η S , where S can be as large as 49 = 7 2 in a 7 ratings system. However, our state space analysis can accommodate this aspect of the model in a feasible way since it relies on computationally efficient methods such as the Kalman filter and associated algorithms. Furthermore, it leads to a significant reduction of the size of the parameter vector θ. Since θ needs to be estimated via the numerical optimization of the likelihood, computation time is also reduced as a result.

The remaining parameters in θ are ρ and α s s. These coefficients can be placed in the state vector ν i as well although in this case the model becomes nonlinear in the state equation.

The linear Markovian process (11) for the state vector ν i is not applicable and we cannot rely

Monte Carlo likelihood evaluation

Given the statistical model specification of the previous subsection, the likelihood function ( 8) can be reformulated by

L(θ | F T ) = N i=1 p (z i |ν i , F i-1 ) p(ν | F T )dν, ( 12
)
where p (z i |ν i , F i-1 ) is given by ( 9) and the model for

ν = (ν ′ 1 , . . . , ν ′ N ) ′ is implied by (11). Both p (z i |ν i , F i-1
) and p(ν|F T ) depend on the parameter vector θ for i = 1, . . . , N . An analytical expression for (12) does not exist and therefore we rely on numerical techniques for the evaluation of ( 12). For this purpose we explore the technique of Monte Carlo integration using the method of importance sampling. The basic idea is simple. First, we simulate M paths of ν from p(ν|F T ) denoted by ν 1 , . . . , ν M where M is a large number. Second, we compute the Monte Carlo estimator of ( 12) given by

L(θ | F T ) = M -1 M m=1 N i=1 p (z i |ν m i , F i-1 ) (13)
where ν m i is the ith element from ν m . The estimator ( 13) is poor since ν m is simulated 'unconditionally' and is therefore likely to make little contribution to the likelihood. A more efficient approach is to simulate from p(ν|z, F T ), but this is not feasible since no analytical expression exists for this density. The idea of importance sampling is to replace p(ν|z, F T ) by the more convenient Gaussian density p G (ν|z, F T ) for simulating ν's. The basic algorithm is then adjusted as follows. First, simulate M paths of ν from p G (ν|z, F T ) denoted by ν 1 , . . . , ν M where M is a large number. Second, compute the Monte Carlo estimator of (12) as given by

L(θ | F T ) = M -1 M m=1 N i=1 p (z i |ν m i , F i-1 ) p(ν m |F T ) p G (ν m |z, F T ) = p G (z|F T )M -1 M m=1 N i=1 p (z i |ν m i , F i-1 ) 1 p G (z|ν m , F T ) , ( 14
) since p G (ν|F T ) = p(ν|F T ) and p G (ν|z, F T ) = p G (z|ν, F T ) p G (ν|F T )/p G (z|F T ).
We refer to this estimator as the Monte Carlo likelihood. The construction of p G (ν|z, F T ) and the evaluation of the different densities is described in detail below.

A c c e p t e d m a n u s c r i p t

Step 1: Simulate paths of ν from p G (ν|z, F T )

To build a device for simulating from the conditional Gaussian density p G (ν|z, F T ), an approximating linear Gaussian model needs to be formulated that represents the joint density p G (ν, z|F T ). This density for the linear Gaussian model ideally resembles the true density p(ν, z|F T ) as close as possible because samples generated from the conditional Gaussian density p G (ν|z, F T ) may then be similar to samples from the conditional density p(ν|z, F T ). An appropriate linear Gaussian model can be obtained using the method described in Durbin and Koopman (2001, Part II) and is based on the linearization of the observational log-density using a second-order Taylor expansion.

In the context of the model described in Subsection 3.1, the basic idea is to construct a linear Gaussian state space model for the series of rating event indicators at event i as given by

{Y 11i , . . . , Y S1i , Y 12i , . . . , Y SKi } ,
where Y ski ≡ Y sk (t i ) and Y sk (t i ) is one or zero to indicate whether a rating event of type s has taken place for unit k at time t i . Such a rating event is triggered by the signal

Z ski ν i = η s + γ ′ s w k (t i ) + α s ψ i which determines the intensity λ sk (t) for t i-1 < t ≤ t i , see Subsection 3.1.
To establish an approximating Gaussian model that relates the signal Z ski ν i to Y ski , we consider the linear Gaussian observation equation The implementation of this procedure is relatively simple. An initial guess for the mode ν needs to be found that is denoted by ν (0) . The linear Gaussian model ( 15) is constructed for

Y ski = c ski + Z ski ν i + ξ ski , ξ ski ∼ NIID (0, C ski ) , (15)
c ski = Y ski -Z ski ν (j) i -C ski Z ski ∇ ln p(z|ν, F T ) i , C ski = -[Z ski ∇ 2 ln p(z|ν, F T ) i Z ′ ski ] -1 , (16)
where

∇ ln p(z|ν, F T ) i = ∂ ln p(z|ν, F T ) ∂ν i ν= ν (j) , ∇ 2 ln p(z|ν, F T ) i = ∂ 2 ln p(z|ν, F T ) ∂ν i ∂ν ′ i ν= ν (j)
.

A new guess of the mode for ν is obtained by estimating the conditional mean of ν conditional on Y for the approximating linear Gaussian state space model ( 15) and ( 11). j+1) . New guesses for the mode are obtained by the KFS based on ( 16) for j = 1, 2, . . . until convergence is reached according to some metric.

Usually convergence takes place after 5 to 10 iterations.

The approximating linear Gaussian model consists of ( 11) and ( 15), with joint density p G (ν, z|F T ) and where ( 16) is evaluated at ν = ν with ν as the estimated mode. We adopt this model to generate conditional samples for ν from p G (ν|z, F T ). Direct sampling from such a high-dimensional Gaussian density requires many high-dimensional matrix operations. These numerical problems can be overcome because the model is formulated as a linear Gaussian state space model. Therefore, the simulation smoothing algorithms of de Jong and [START_REF] De Jong | The simulation smoother for time series models[END_REF] or [START_REF] Durbin | A simple and efficient simulation smoother for state space time series models[END_REF] can be used to generate conditional samples for ν, denoted as

ν m for m = 1, . . . , M .
Step 2: Compute the Monte Carlo likelihood ( 14)

Given a set of simulated samples from p G (ν|z, F T ) ≡ p G (ν|Y, F T ) and denoted by ν m , the computation of the Monte Carlo likelihood ( 14) is relatively simple. The Gaussian density

p G (z|ν, F T ) ≡ p G (Y |ν, F T
) is conditional on ν and its expression is well-known for the linear model (15). Further, the observation density of interest p(z i |ν i , F i-1 ) is given by equation ( 9) and can also be computed straightforwardly.

The Monte Carlo likelihood is then maximized with respect to θ for a particular choice of M . The maximization can be carried out by a numerical optimization procedure. For example, a quasi-Newton method can be used for this purpose. To ensure a likelihood surface that is continuous (or smooth) in θ, the same random numbers are used for the sampling in Step 1

A c c e p t e d m a n u s c r i p t

of the M signals from p G (ν|z, F T ). The optimization procedure requires an initial estimate of θ that is chosen ideally in the neighborhood of its final estimate. However, even when the initialisation is poor, several modifications in numerical optimization procedures exist that make global convergence attainable.

Step 3: Smoothed estimates of the state vector

The state vector ν i contains fixed unknown coefficients and dynamic latent processes. Estimating the state vector for each i leads to estimates of regression parameters and latent processes such as ψ i . A straightforward estimate of the state vector, given the data, is obtained by weighting each simulated state vector ν m i by its contribution to the likelihood function, that is

νi|N = M m=1 w m × ν m i M m=1 w m , (17)
where

w m = N i=1 p (z i |ν m i , F i-1 ) p G (z|ν m , F T ). ( 18
)
Standard errors for νi|N are obtained by taking the square root of

M m=1 w m × (ν m i ) 2 M m=1 w m -(ν i|N ) 2 . ( 19
)

Implied Transition Matrices

We For the MLFI model, a parametric form for G k (t) conditional on observed regressors and an unobserved factor is assumed. In [START_REF] Aalen | An empirical transition matrix for non-homogeneous Markov chains based on censored observations[END_REF], by contrast, G k (t) is left completely unspecified under the assumption that duration and self-excitation effects are absent.

We therefore use the Aalen Johansen estimator in our empirical section as a benchmark for evaluating model adequacy. We assume that the elements of G k (t) are adapted to the observable filtration F t . In this situation P k (T, T +∆) becomes a random variable, and we want to compute its expectation conditional on F T . This expectation can be interpreted as the transition matrix over the interval [T, T + ∆],

Pk (T, T + ∆) = E [P k (T, T + ∆) | F T ] = E   T +∆ T (I + G k (t)dt) F T   . ( 21
)
To estimate the conditional expectation in (21), we suggest two approximating schemes.

In the first scheme, we generate a large number of high frequency (say daily) paths from the process ψ(t), compare ( 2) and (3) in Section 2. Using these paths, a consistent estimator for Pk (T, T + ∆) is given by

Pk (T, T + ∆) = 1 M M m=1 T +∆ T (I + G m k (t)dt) , (22)
where G m k (t) denotes unit k's realized matrix of intensities for replication m. This approximation skips the fact that in our empirical model specification, ψ(t) only jumps at event times. This last assumption, however, is not material in the model's specification, and jumps at higher frequencies (such as every day) are easily allowed for. The second scheme to obtain estimates of annual transition probabilities uses the bootstrap. Here, we build on the empirical model's assumption that the common factor ψ(t) only jumps at event times. We start with the estimates of the unknown model parameters and the smoothed estimates of the latent process, E[ψ i |F T ] for i = 1, . . . , N (T ). Next, we simulate a large number M of possible future sample paths over the [T, T + ∆] interval for the full panel of K firms as well as for the unobserved risk factor ψ i . 9

8 See [START_REF] Gill | Product integration[END_REF] for an exposition on product integration. The product integral is the continuous counterpart of the standard, discrete product operator, just as the integral is the continuous counterpart of the summation operator. Informally, the product integral of a function f (t) over the interval [T, T + ∆] is

T +∆ T (1 + df (t)) = lim n→∞ n i=1 (1 + f (t i ) -f (t i-1
)) for a partition T = t 0 < t 1 < . . . < t n = T + ∆. 9 If weakly exogenous covariates were included in equation ( 1), then an auxiliary model is needed to forecast the future path of these covariates (as mentioned in the introduction). One resulting possibility is the estimation of scenario forecasts.

A c c e p t e d m a n u s c r i p t

The pooled process over [T, T + ∆] for replication m provides a partition T = t m 0 < t m 1 < . . . < t m n = T + ∆, over which the product integral can be factored, see Andersen et al. (1993, p. 91). For the empirical model in Section 6 these calculations become particularly manageable.

The estimates of Pk [T, T + ∆] can be used to compute several interesting risk measures. For example, one can compute the average transition probabilities over a specific time interval for a portfolio of firms,

P [T, T + ∆] = K -1 K k=1 Pk [T, T + ∆].
One can also compute (non-linear) functions of the default probabilities in Pk [T, T + ∆] to obtain direct estimates of capital requirements according to the official Basel II regulations. This is especially interesting if one does not average over simulations of ψ(t), but considers quantiles instead. See the forecasting exercise in Section 6. This also allows one to consider stress scenarios in terms of extreme ψ(t) realizations.

So far, we have discussed how to obtain estimates of transition matrices for Markov processes.

If we move on to semi-Markov processes, the equations become more involved. In the current case of semi-Markov behavior and a common risk factor ψ(t), the relevant equations are worked out in [START_REF] Monteiro | Nonparametric Estimation for Non-Homogeneous Semi-Markov Processes: An Application to Credit Risk[END_REF]. If we follow the first approximation scheme above, we then have to solve a system of Volterra integral equations for a fixed path of ψ(t) to obtain the transition matrix. This matrix subsequently has to be averaged over different simulated paths of ψ. Alternatively, we could follow the second approach and simulate the complete panel of firms to obtain a realized path of ψ that only jumps at event times. This could then be used in the Volterra equations to obtain a transition matrix. The resulting matrix again would have to be averaged over many replications. Working out the finite sample properties of these different approaches is beyond the scope of the current paper. Here we first focus on estimating the empirical dynamics of systematic credit risk. We leave the implication of these dynamics for one-year default probabilities and risk measures for future research.

Simulation Results

To assess the performance of the Monte Carlo maximum likelihood method in a controlled environment, a simulation experiment is carried out. The modeling framework resembles closely one of the model specifications in the empirical study of Section 6. We consider 7 rating classes.

The states can be interpreted as the familiar grades AAA, AA, A, BBB, . . . , CCC, and default.

Default is modeled as an absorbing state. In the simulation section, we restrict the intensities

(t) = R sk (t) • exp [η s + α s ψ(t)] ,
where ψ(t) is a step function that jumps at the endogenous event times t i as in [START_REF] Gagliardini | Stochastic migration models with application to corporate risk[END_REF]. The firm heterogeneity in this specification enters through the different parameters η s for the different transition types s = 1, . . . , S. Another source of heterogeneity is the latent process ψ(t) that can be interpreted as the (unobserved) macroeconomic effect. The benchmark model in this simulation exercise abstracts from duration dependence by setting H sk (t) ≡ 1. This assumption is relaxed in the empirical section. Further parsimony is introduced by setting α s = α down < 0 for downgrades, and α s = α up > 0 for upgrades. The parameter values used for the simulation can be found in the first column of Table 1.

<INSERT TABLE 1 ABOUT HERE>

Note that the number of parameters, even in this simple model specification, is large. For 7 rating classes, we have 49 possible rating transitions. Given the underlying data generating process (dgp) in this simulation study, some of the transitions are extremely unlikely, e.g., from AAA to default. As a result, in a particular simulation run there may be no transitions of this type. In such cases, the corresponding η s parameter is not estimated. This means that not all η s parameters can be estimated for every simulation. Table 1 only reports the simulation results for those η s parameters for which a reasonable (50 out of 500) number of simulations exists for at least one of the four data generating processes (DGPs) presented in Table 1. The model settings we consider are K = 70 and K = 700 combined with autoregressive parameters ρ = 1 and ρ = 0.9.

For both K = 70 and K = 700, a panel of firms and rating transitions is generated as follows.

At time t 0 = 0, the sample contains an equal number of firms in each rating category. The unobserved process ψ(t) is initialized at zero. Given the parameters, this completely specifies the intensities up to the event date t 1 . For the time interval (t i-1 , t i ], the intensity of the pooled process is defined by

λ * (t i ) = K k=1 S s=1 λ sk (t i ), (23)
with λ * (t 1 ) applicable over the first spell (t 0 , t 1 ]. The length of any spell in the pooled process can therefore be drawn from the exponential distribution with intensity parameter λ * (t i ). Given the durations of the spells (t i-1 , t i ] for i = 1, . . . , N (T ), the firm experiencing a rating event is drawn from the univariate Multinomial{π 1 (t i ), . . . , π K (t i )} distribution where the probability

A c c e p t e d m a n u s c r i p t

of drawing unit k is given by

π k (t i ) = [λ * (t i )] -1 S s=1 λ sk (t i ), k = 1, . . . , K. (24)
Next, the type of rating event for unit k is drawn from the multinomial distribution with the probability of state s being drawn for unit k given by

π sk (t i ) = S s=1 λ sk (t i ) -1 λ sk (t i ), (25)
for s = 1, . . . , S and k = 1, . . . , K. If the event is a default, the dummy variable R ik (t) jumps to zero. Finally, the unobserved common risk factor ψ i = ψ(t i ) is updated using ( 2) with ρ = 1 or ρ = 0.9 and where the disturbances ε i , i = 1, . . . , N (T ), are drawn from a standard normal distribution. This process is repeated until all units have entered the absorbing default state, or until the event time t i exceeds the maximum period of 25 years. For each panel size, we performed 500 replications of the simulations. All calculations in this paper were performed using the Ox matrix programming language of [START_REF] Doornik | Object-oriented matrix programming using Ox[END_REF] and the estimation and smoothing routines in the package SsfPack of [START_REF] Koopman | Statistical algorithms for models in state space form using SsfPack 2.2[END_REF].

The simulation results for the Monte Carlo maximum likelihood procedure discussed in Section 3 are shown in Table 1. We first concentrate on the model with a random walk factor (ρ = 1). We see that many of the parameters are estimated accurately for a panel with 70 firms.

The parameters that are estimated less accurately, correspond to larger rating transitions (e.g, AAA to BB or AA to CCC). As the larger rating transitions are much less likely, the Monte Carlo averages of the corresponding η s parameters are based on less replications and, therefore, less accurate themselves. By contrast, the presented averages for the smaller rating transitions are all very close to the true parameter values. The accuracy of the η s parameters increases further if the number of units is increased to K = 700. This is due to the fact that for a larger panel, we will observe more types of transitions in the same period of 25 years. The increase in precision also holds for the α s parameters. By considering the average estimate, it appears that the estimator for α s is somewhat biased towards zero. This implies the magnitude of the common risk factor is underestimated. Consequently, it is more difficult to find significant evidence of such a factor in the empirical section later on. The strength of the common factor actually found in the empirical section might thus be a lower bound on its true value.

If we consider a stationary specification for ψ with ρ = 0.9, the results are similar. Parameters for small rating transitions are estimated accurately. Parameters for the larger rating transitions are more difficult to estimate for the smaller panel due to the limited number of observations. For the large panel, the bias in the estimates of α s appears smaller than for a

A c c e p t e d m a n u s c r i p t

non-stationary ψ(t). For the stationary model, we estimate the additional AR(1) parameter ρ.

For the small panel of K = 70, the average estimate is near its true value of 0.9. The slight negative finite-sample bias is well-known from the linear model context. If the panel size increases to K = 700, the estimate of ρ remains stable, while its Monte Carlo standard deviation decreases somewhat. The mild improvement is largely due to the parameterization chosen. We increase the cross-sectional dimension of the panel while keeping T at 25 years. This increases the number of event times, without increasing the calendar time available. Given the scaling of the parameters of the autoregressive process with the length of the spell intervals, there is limited additional information on the long-term dynamics of the AR process. There is only a moderate precision gain caused by a more precise estimate of the signal given the increased number of events.

As explained in Section 3 (Step 3), we can use smoothing techniques to obtain an estimate of the unobserved ψ(t) factor. Figure 1 illustrates the result for a single 'representative' simulation. For this simulation, we plotted the true value of ψ(t) against its estimated (smoothed) counterpart using our importance sampling scheme. The 95% confidence bounds are also provided. The algorithm clearly performs adequately in recovering the characteristics of the true, unobserved ψ(t) process from the observed data. As expected, the true ψ(t) is much more volatile at high frequencies than its smoothed counterpart. Local and global peaks and troughs of the series, however, appear correctly positioned in calendar time. The true process also falls inside the 95% confidence interval most of the time.

<INSERT FIGURE 1 ABOUT HERE>

6 Empirical Results

Data

The data consist of rating transitions obtained from Standard & Poor's. The rating histories of all issuers are recorded in the CreditPro 7.0 database. The sample period is from the end of 1980 (the left-censoring time point) until June 2005. We express the durations of the pooled process as a fraction of the business year. Note that there may be multiple rating events on a single day. This is captured by the variables Y sk (t i ) in ( 6). The rating histories in the data set distinguish between more than 18 different rating classes. To illustrate our methodology, we consider only seven broad classes, namely AAA, AA, A, BBB, BB, B, CCC (and lower), and default. This produces 49 possible rating transitions. Out of these 49, 42 are observed in

A c c e p t e d m a n u s c r i p t

the current sample. In terms of model ( 9), we therefore have S = 42, K is almost 7000, and N > 25000 (or 14000 if the ψ(t) process does not jump at transitions involving the non-rated class). Even with these sample sizes, the model can be implemented on a standard desktop PC.

New firms enter the sample when they receive a rating for the first time. Firms leave the sample when they enter the default state or when their rating is withdrawn. However, S&P continues to track firms whose ratings are withdrawn. It is notified in the database when such firms default at a later stage. This should substantially mitigate any biases caused by strategic behavior of firms in maintaining a rating at S&P. If a firm first enters the non-rated class and later defaults, we skip the transition to the non-rating class. In total, we observe 7000 firms, though many of these firms are only observed over part of the sample. The number of rated firms increases over time. Accounting for attrition (defaults and transitions to non-rated) and sample extension, the time series average of the number of firms available at any moment is around 2200.

Given our interest in the clustering effect of defaults, we first clean the database for alternative forms of clustering. We use two filters. First, we check the data for firms that have long histories of coincident rating increases and decreases. Typical examples of these are firms that have merged during the time of the database. In such cases, we leave out one of the two firms from the time of the merger onwards. Second, there appeared to be some policy of the rating agencies in clustering re-ratings, e.g., centered around meeting times of the committee.

To account for this, we Winsorized the 'number of events per day' to 3 by replacing Y sk (t i ) by

Y sk (t i ) • min(3, Y s (t i ))/Y s (t i ), with Y s (t i ) = k Y sk (t i ).
Experiments with Winsorizing values between 1 and 5 yielded similar results. The Winsorizing procedure has the largest effect on the estimation of the dynamics of the latent component, summarized by the AR parameter ρ.

Erratic clustering due to rating agencies' policies, e.g., around committee meetings, corrupts the dynamics of ψ and causes a downward bias in the estimate of ρ. This effect is mitigated considerably by the Winsorizing procedure. In future research, these anomalous clustering effects can be dealt with differently, e.g., by trying to model them explicitly.

Some descriptive features of the data are as follows. The pooled process has a high intensity of migrations, resulting in an average duration between transitions, births, or withdrawals of 1.2 days. There is a large number of downgrades and upgrades. The counts of transition events is given in Table 2. It is clear that most transitions take place to adjacent rating categories.

Some of the transitions are very rare, e.g., large down-grades or up-grades. In order to check the sensitivity of our results to these rare events, we also performed a robustness check where we only incorporated transition types with more than 20 events. The results are not sensitive The data are visualized in Figure 2. In order to keep the number of graphs manageable, we temporarily clustered the ratings further into investment grade (AAA-BBB) and subinvestment grade (BB-CCC). The upper two plots in Figure 2 show the number of investment grade downgrades and defaults on a daily basis since December 31, 1980, respectively. We see that downgrades and defaults tend to cluster in time. This follows from the concentration of vertical lines, which originate dark and bright areas along the horizontal (time) axis. The lower two plots in Figure 2 contain the number of sub-investment grade upgrades and downgrades on any given day, respectively. Interestingly, the plots complement each other. When downgrades and defaults cluster, upgrades are more sparse, and vise versa. This suggests that the model specification used with a single common risk factor ψ(t) might already captures the most salient features of the data. We also experiment, however, with multiple factor models at a later stage in the analysis.

<INSERT FIGURE 2 ABOUT HERE>

Homogeneous continuous-time Markov chain model

ηs = ln N i=1 K k=1 Y ski -ln N i=1 K k=1 t i R ski . ( 26
)
Table 3 presents the parameter estimates based on (26) as well as the estimates obtained by using the Monte Carlo methods of Section 3. As expected, the parameter estimates are almost identical for all transition types for which the number of observations is sufficiently large. But even for the transition types with only one observation, the importance sampling estimates (SML) never differ more than 10% from their closed form counterparts. This is well within a bound of two standard errors. As stated before, we also performed a robustness check only including the transition types with 20 observations or more. Here the differences between the sampling and analytic approach appear negligible. Also note that the estimate of the likelihood differs by less than 0.01% from its analytic counterpart. Again, this difference shrinks to zero if only transition types with more than 20 observations are used.

A c c e p t e d m a n u s c r i p t

<INSERT TABLE 3 ABOUT HERE>

Estimation results for the MLFI model

We start our empirical exercise with the introduction of a single random walk component ψ i capturing systematic credit risk. This implies the AR parameter ρ is set to unity. The estimation results are in Table 4, model B. Due to the restriction ρ = 1, ψ(t 1 ) and η s for s = 1, . . . , S(= 42) are not jointly identified. Therefore we start the latent process ψ(t) at ψ 1 = ψ(t 1 ) = 0. This means that ψ(t) can be interpreted as a relative credit index compared to its starting level in December 31, 1980.

For the intensity specification (1) in model B, we set α s = -ᾱ < 0 for downgrades and α s = ᾱ > 0 for upgrades. Since ψ i is interpreted as the (unobserved) credit cycle, these sign restrictions on α s imply an increase in the probability of downgrades and defaults if ψ i is negative, and a simultaneous decrease in the probability of rating upgrades. Conversely, if ψ i is positive, it leads to an increase in the probability of firms being upgraded. The sign restrictions are relaxed later on.

It is worth mentioning here that all models reported in Table 4 have an additional 42 (unreported) η s parameters that need to be estimated. This is the number of non-zeros in Table 2. Given space constraints, we do not report them here. They are however a part of the estimation problem as sketched in Section 3.

<INSERT TABLE 4 ABOUT HERE>

The likelihood increase from model A (no common factor) to model B (a single common random walk ψ(t)) is large: almost 350 points upon adding one parameter. We conclude that there is ample evidence of a systematic credit risk component in the data at hand. This also appears from the standard error of the estimate of α s .

The above result, however, hinges on the assumption of a single sensitivity parameter for all transition types. This may not hold empirically. For example, upgrades can be less sensitive to common risk factors than downgrades and defaults, see the discussion in [START_REF] Kavvathas | Estimating credit rating transition probabilities for corporate bonds[END_REF]. As a first interesting refinement of model B, we allow in model C the single systematic credit risk factor ψ(t) to have a different impact on downgrade intensities compared to upgrade intensities.

We again restrict the α s coefficient of downgrades to be negative. For the upgrade α s , no sign restriction is imposed. We find that the improvement of model C compared to B is again statistically significant: an almost 18 point increase with one additional parameter. The sign of the upgrade α s is intuitive. The ψ(t) factor has the opposite effect on downgrade versus upgrade

A c c e p t e d m a n u s c r i p t

intensities. Moreover, the systematic risk in downgrade intensities appears to be significantly higher than in upgrades. The α s for upgrades (1.61) is about half of that for downgrades (3.13). This result is in line with the results of [START_REF] Das | Correlated default risk[END_REF] and [START_REF] Kavvathas | Estimating credit rating transition probabilities for corporate bonds[END_REF] based on observed macroeconomic variables as proxies for the common risk factor.

Given the significant difference between systematic risk in upgrades and downgrades, we make a further distinction. We allow for a different systematic risk sensitivity for investment grade (AAA-BBB) companies that (i) upgrade, (ii) downgrade to investment grade, and (iii) downgrade to sub-investment grade. Similarly, for sub-investment grade companies we allow for a different sensitivity for downgrades, upgrades within the sub-investment grade class, and upgrades to the investment grade class. The results are presented as model D in Table 4. The likelihood increase from model C to D is again statistically significant: a 7.9 points increase upon adding 4 parameters. This gives a likelihood ratio test statistic of 15.8 with p-value 0.003.

Again, we see the largest systematic credit risk sensitivity for sub-investment grade downgrades and defaults, followed by small and large investment grade downgrades. The smallest systematic risk exposure appears to be for investment grade upgrades. Though still significant, it is less than 50% of any of the other sensitivities.

The left-hand panel in Figure 2 presents the smoothed estimate of the factor ψ(t) for model D. The factor is obviously low (so credit risk is high) in the mid 80s, early 90s, and early 2000s. Given the 95% confidence bands, at each of these three episodes the credit risk factor is significantly negative. Given the signs of the α s coefficients, this implies that downgrade and default intensities were significantly higher in those periods. The converse holds for the upgrade intensities.

<INSERT FIGURE 3 ABOUT HERE>

In order to investigate whether a single risk factor suffices to model systematic credit risk, we allow for a different factor ψ s (t) for each of the 6 blocks used in model D. The results are presented as model E. The likelihood increases by about 110 points from model D to E. This cannot be a formal test, however, as the models are non-nested. The results are, however, highly suggestive of different systematic risk factors being relevant for different transition types. This also appears from the smoothed estimates of the common risk factors to be presented later.

By looking at the sensitivity parameters α s , we can get a first idea of which transition types might show a deviating systematic risk exposure. For investment grade companies, we see large increases in the (absolute) magnitude of the α s coefficients compared to model D for rating migrations within the investment grade class. This signals that some clustering may take place, but that the clustering need not coincide with the systematic default risk component.

A c c e p t e d m a n u s c r i p t

The same holds for limited sub-investment grade upgrades, i.e., within the sub-investment grade class. The coefficient almost doubles compared to model D. In short, the results are suggestive of a difference in clustering effects between defaults and large rating movements on the one hand, and limited re-ratings on the other hand. The difference might be attributable to the through-the-cycle policy of rating agencies. Here large rating movements and especially defaults might be attributable to economic up and down-turns, whereas limited re-ratings might be more company, industry, or rating policy specific, thus causing a different type of clustering.

So far, we imposed the restriction of a unit root on the latent component(s) ψ(t). Given the short average spell length of 1.2 days, this should give a good in-sample description of the systematic risk factor. Out-of-sample, however, the restriction of a unit root may be less realistic. It implies that even when the credit risk factor is very low, it can still move up or down with equal probability. This may result in over-estimates of credit Value-at-Risks (VaRs) for risk management purposes. To check the sensitivity of the results for the unit root assumption, we re-estimate all models using stationary AR(1) models. We use the specification in ( 2) and

(3) in Section 2. The autoregressive parameter and the innovation variance are now tuned automatically to the spell length. Short durations have a high persistence parameter and a small innovation variance. By using this specification and the data at hand, the ρ parameter can readily be interpreted as an autoregressive parameter at an 'annual' (260 business days) frequency.

The results for the single factor model with stationary latent component and common sensitivity parameter α s are presented as model F. The likelihood increases by 4.5 points compared to model B upon adding one parameter. We cannot say, however, whether this is statistically significant, because this is a unit root test in a class of non-linear models. Working out the statistical properties of this test is well beyond the scope of the current paper. A similar increase of 4.1 points is seen for the case with unrestricted α s s, compare model G to D.

The estimated value of the AR parameter in models F (0.38) and G (0.37) is very similar.

The estimates are much smaller than those in Koopman et al. (2005a). The latter, however, are not based on an intensity framework and use annual default counts. The estimates in [START_REF] Mcneil | Bayesian Inference for Generalized Linear Mixed Models of Portfolio Credit Risk[END_REF] are based on quarterly counts and are also somewhat higher, around 0.6 to 0.7. Comparing models G and D, we see that the reduced persistence of ψ by the lower autoregressive parameter causes the sensitivity coefficients α s to increase. This is due to the fact that the model also tries to match the unconditional in-sample variance in default intensities.

A graph of the smoothed AR(1) factor for model G is presented in the right-hand panel of and troughs are roughly similar. There are three main differences: the location, scale, and standard error bands. It must be kept in mind, however, that ψ(t) enters the intensity process as η s + α s ψ(t). The location and scale effect is thus partly off-set by the differences in estimates of α s and η s between the models. On average, the η s in model D is about 0.4 lower than in model G for downgrades and defaults. For upgrades, the average difference is -0.2. Similarly, the increase in the (absolute) magnitude of the α s s from model D to G off-sets the smaller scale for the stationary ψ(t) factor in model G. As a result, given the sign differences in α s between upgrades and downgrades, the difference in impact of the systematic risk factor between models D and G appears negligible.

In model H we again look at a multiple factor model. The six ψ factors now all follow separate (independent) stationary autoregressive processes with parameters ρ s . Again we see the familiar effect that lower persistence in ψ is matched by a higher α s . More interestingly, however, are the differences in the estimated autoregressive parameters between transition types. The persistence of the sub-investment grade downgrade factor (0.40) is similar as in model G. The persistence in within investment grade downgrades and within sub-investment grade upgrades, however, is virtually absent. The common factor here does not appear to be sticky at all. Though not apparent immediately, this implies that there is not really a systematic credit risk factor for these types of transitions. If ρ s is close to zero, the ψ factor essentially collapses to another source of firm (and time) specific error given its frequency of updating (i.e., at every event). Consequently, the durations in the pooled process are completely idiosyncratic, though no longer (conditionally) exponential. By contrast, we see that the persistence of the common factor for large up or downgrades is high (0.93 and 0.92). This again underlines our earlier finding that the common risk factors may be different across the different transition types. This is also visualized in Figure 4. The default and downgrade factors ψ(t) for IS down and SS down clearly reveal a different pattern from the upgrades, II up and SI up .

To conclude our set of model specifications, we also consider a semi-Markov model. In particular, we allow the intensities to depend on the time spent in the current rating. We capture this effect through the baseline hazard functions H sk in (1). We allow for a different baseline hazard for investment grade versus sub-investment grade companies. The baseline hazards are flexibly parameterized by cubic splines with knots at 0, 1, 3, and 6 years. The spline passes through (0,0) in order to identify both the spline parameters and the η s coefficients. The results are presented as model I in Table 4. The increase in likelihood from model H to I is again statistically significant: 76.9 points upon adding 6 parameters. The persistence parameters ρ and sensitivity parameters α s are very similar between the two models. This is probably due to the rating momentum effect: firms that have been downgraded recently are more likely to downgrade in the near future because of the smooth rating revision policy by rating agencies. A second peak emerges after 3 to 5 years. Here, default typically occurs because initial debt taken by these companies has to be rolled over. This causes a peak in defaults because of some sub-investment quality firms getting into financial distress at that time.

Forecasting transition probabilities

To illustrate the impact of fluctuations in the estimated latent risk component for the implied 1-year transition/default probabilities, a recursive (partially) out-of-sample forecasting exercise is carried out. The details of the forecasting study are as follows. We consider the data windows Dec 1980-June 1992, Dec 1980-June 1993, . . . , Dec 1980-June 2005. For each data window, we estimate the parameters η s , the smoothed risk factor estimate ψ and its covariance matrix.

The estimates of the loadings α s and ρ s are kept fixed at their estimates from Table 4 in this illustration. This reduces the computations needed. Note that in the forecasting exercise, we do not Winsorize the realizations in the forecasting period as we did during the estimation stage.

Given the estimates in a specific year, we simulate daily paths of the common risk factor(s) ψ s (t) over the subsequent year, starting with a draw from the smoothed distribution. For each path ψ, the product integral in (20) is computed. This gives us a conditional (on ψ) annual transition matrix. Finally, we compute the 99% quantile of the default probabilities from these simulated transition matrices for each initial rating category. The 99% quantile is taken over the different simulated paths of ψ. We do so for models A, D, E, G, and H from Table 4. As mentioned in Section 4, if a semi-Markov specification were used like model I, computing the annual transition matrix would be much more involved and dependent on the current portfolio composition. See [START_REF] Monteiro | Nonparametric Estimation for Non-Homogeneous Semi-Markov Processes: An Application to Credit Risk[END_REF] for more details on the semi-Markov case.

To benchmark all model forecasts, we also compute the nonparametric Aalen Johansen (AJ)

estimator of the transition probabilities over each year. This empirical transition matrix can be seen in the current context as a proxy for the 'true' or realized 10 transition probabilities, see 10 The term realized is chosen deliberately here. Its well-known analogue is the nonparametric computation A c c e p t e d m a n u s c r i p t [START_REF] Lando | Analyzing rating transitions and rating drift with continuous observations[END_REF]. Figures 6 and7 presents the results.

<INSERT FIGURES 6 AND 7 ABOUT HERE>

The forecasts from the MLFI model move much more in accordance with the nonparametric ex-post AJ estimates than the HCTMC model (HMC). Note that though the HMC does not include a dynamic component,the forecasts still adapt over the years due to the recursive nature of the estimation procedure for η s . The changes are nevertheless insufficient to keep up with the changes in the true default rates (see AJ curve).

The 99% quantiles for the MLFI models with 1 or 6 AR risk factors appear reasonable upper bounds for all rating categories except single A. By contrast, the forecasted default probabilities for the non-stationary models are much larger. This is due to the random walk nature of the factors. Even if the common risk factor is already very low (so default probabilities are high), it is still equally likely that the factor goes further down or up. Excessively large forecasts for default probabilities are the result and these subsequently lead to inefficient (though highly prudent) capital buffers at most times. For the stationary risk factors, low factors tend to mean revert, resulting in less extreme estimates of future default probabilities. This comes at the cost of an increased probability of violating the 99% quantile bounds around early 2000s

for some initial rating categories.

An interesting feature of the MLFI specifications is its increase in predicted default probabilities before the actual crash. This holds in many cases, especially around the record default years in the early 2000s. The HMC model with recursively estimated parameters, however, obviously lags the dynamic pattern in defaults. This suggests that for risk management purposes, models that account for the dynamic pattern in systematic risk factors are very important.

Our forecasting procedure produces much more than only the default rates. Also all rerating probabilities are available. To get a (parsimonious) impression of how well the different models perform with respect to all the remaining transition probabilities, a summary statistic of the transition matrix is needed. [START_REF] Jafry | Measurement, estimation and comparison of credit migration matrices[END_REF] It is clear that the current forecasting exercise can be extended further. Different specifications for the MLFI could be tried. Also, the semi-Markov property could be included at the cost of a significant increase in the required computations in line with [START_REF] Monteiro | Nonparametric Estimation for Non-Homogeneous Semi-Markov Processes: An Application to Credit Risk[END_REF].

Moreover, dollar based metrics could be easily developed. For example, given the conditional forecasts of transition matrices from the product integral expression, one can construct (conditional) forecasts of capital requirements in line with the New Basle Capital Accord, see also [START_REF] Jafry | Measurement, estimation and comparison of credit migration matrices[END_REF]. We leave such extensions to future applications.

Conclusion

In this paper we have motivated and introduced a multi-state latent factor intensity (MLFI) model for credit rating transitions. The model can be regarded as a generalization of the latent factor intensity point process introduced in [START_REF] Bauwens | Dynamic latent factor models for intensity processes[END_REF] to a situation with multiple origin and destination states. However, the econometric issues related to this generalization are intricate and the computational consequences are severe. We have discussed the details for the estimation of the MLFI model using a Monte Carlo maximum likelihood procedure that consists of a combination of importance sampling techniques and state space methods, as outlined in Durbin andKoopman (1997, 2001). It is shown that this procedure can be extended successfully to a multivariate class of non-Gaussian models. A simulation study was carried out to show that the estimation procedure works well in recovering the parameters of the MLFI model. Finally, we have applied the model to a real-world dataset of credit rating migrations. A significant common risk factor in credit rating migrations is found.

The impact of this risk factor is higher for downgrades than for upgrades. This empirical result suggests that upgrades are more subject to idiosyncratic shocks than downgrades. This finding is consistent with the conclusions in the earlier studies of [START_REF] Kavvathas | Estimating credit rating transition probabilities for corporate bonds[END_REF] and [START_REF] Das | Correlated default risk[END_REF]. Furthermore, our results pointed towards important potential differences in the systematic risk of different transition types. This can have important implications for single factor models for risk management, as they are currently widely employed in the industry.

Generalizations with respect to the current empirical specification are easily incorporated in the structure of the MLFI model. For example, the general specification allows for the inclusion of observed firm-specific and economic variables as in Duffie et al. (2006b) and [START_REF] Koopman | Credit Cycles and Macro Fundamentals[END_REF], self-exciting processes, and additional dynamic components. A further interesting feature of the model is that it produces a high-frequency credit cycle index estimated directly from default and rating migration data. Finally, within the model specification proposed in ). The transition types are sorted on their number of observations. The first two columns provide the rating transition type. The third column indicates the number of steps taken for this transition type in the 7 grade system. The fourth column contains the number of observations for this transition type. The column headed ML-η s contains the closed form maximum likelihood (ML) estimates from (26). Its ML asymptotic standard error is in parentheses in column six. The Monte Carlo maximum likelihood estimates SML-η s and their simulated standard errors (using 100 samples) follow in columns seven and eight. The column ∆η s and ∆s.e. give the difference between the ML and SML estimates, and between their standard errors, respectively. ψ s,i-1 + σ si ε si , with ε i i.i.d. standard normal random variables, t i for i = 1, . . . , N the event times of the pooled process measured in fractions of the 260 day business year (instead of in days as in (3)), and σ 2 si = (1 -ρ ti-ti-1 s

)/(1 -ρ 2/260 s

). For identification, ψ s0 = 0. The log baseline hazard function h s depends on the time spent in the current rating and follows a cubic spline function with knots at 0, 1, 3, and 10 years. Given the presence of the η s factor, we restrict h s (0) = 0. Spline values at spell lengths exceeding 10 years are set to the spline value at 10 years. The η s are allowed to be different for all transition types (yielding almost 49 parameter estimates that are not presented here). The log baseline hazard h s only differs between investment grade (AAA-BBB) and sub-investment (BB-CCC) grade firms (and has 3 parameters per spline). We allow for a different α s for within investment grade upgrades (II up ) and downgrades (II down ), within sub-investment grade upgrades (SS up ) and downgrades (SS down ), and across investment-sub-investment grade upgrades (SI up ) and downgrades (IS down ). A similar block structure is imposed on the ψ s and ρ s if we allow for multi-factor models. For models A-E, the latent component(s) ψ s (t) are random walks, i.e., ρ s ≡ 1. Standard errors (in parentheses) and log-likelihood values are based on 1,000 importance samples. The baseline model and the simulation set-up are the same as explained in the note to Table 1. The thick, solid curve is the smoothed estimate of ψ(t). The thin solid curve is the true ψ(t) path in a typical simulation. The dotted lines give the 95% confidence band of the smoothed ψ(t). 4 are models D (1 common random walk (RW) factor), G (1 common autoregressive (AR) factor), E (6 independent RW factors), and H (6 independent AR factors). The Aalen-Johansen (AJ) estimates for each year are also plotted as a proxy for the observed transition rates. There are 4 plots for the rating categories AAA, AA, A, and BBB, respectively. See the note to Figure 6 for further explanation. There are 3 plots for the rating categories BB, B, and CCC, respectively. The lower right-hand plot provides the average singular value measure for transition matrices as proposed in [START_REF] Jafry | Measurement, estimation and comparison of credit migration matrices[END_REF] as a general measure for the forecasting accuracy of transition matrices.

Model

{1, 2 ,

 2 provide a short discussion of this curse of dimensionality problem.The multi-state feature of the model is represented as a set S of transition types, S = . . . , S}. For example, in the case of three rating classes (AAA,AA,A), s = 1 denotes a downgrade from AAA to AA, s = 2 from AAA to A, s = 3 an upgrade from AA to AAA, . . ., up to s = S = 6 an upgrade from A to AA. Next, define the right-continuous counting

  time t are that unit k has the incorrect current initial rating, has defaulted, or dropped out of the sample earlier for other reasons.

  where t -t ik denotes the backward-recurrence time of unit k with respect to its past ith transition moment. The function H sk (•) can be any non-negative function of its arguments.

  efficient Kalman filter methods as we do below. The treatment of nonlinear state processes is more involved and computationally more demanding. The details of the estimation procedures for ν i and θ in the current framework are presented in the next subsection.

  for s = 1, . . . , S, k = 1, . . . , K and i = 1, . . . , N , where scalar constant c ski and scalar variance C ski are considered as auxiliary and unknown variables that need to be constructed in a consistent fashion as is shown below. The observation Y ski is linear in vector ν i and modeled by the linear Gaussian process (11). Therefore, observation equation (15) and the dynamic latent process (11) make up a standard linear Gaussian state space model, seeDurbin and Koopman (2001, part I) for a detailed discussion on this class of models.The constant c ski and variance C ski of the observation equation (15) are constructed in such a way that the conditional density of the model of interest p(z|ν, F T ) and the conditional density of the approximating model p G (Y |ν, F T ) have the same mode for ν, where Y = (Y 111 , . . . , Y SKN ) ′ . The joint solution for c ski and C ski to obtain the mode denoted by ν can be obtained recursively, see the treatment in Durbin and Koopman (2001, Chapter 11).

  now turn to the issue of estimating the transition matrix given the Monte Carlo maximum likelihood estimates of the parameters. Typical examples include 1-year transition matrices as the ones published by Standard & Poor's and Moody's. We start by recalling the connection for unit k between the infinitesimal generator matrix G k (t) and the implied matrix P k of transition probabilities for a continuous-time finite-state Markov process. The case of semi-Markov processes is dealt with later on.7 The matrix G k (t) contains the hazard rates for each origin and destination state combination. In particular, the (i, j)th element of G k (t) equals λ sk (t) for s corresponding to a transition from origin state i to destination state j. The diagonal elements of G k (t) are such that the rows of G k (t) sum to zero. Consider an interval [T, T + ∆]. Then the matrix of transition probabilities over the interval [T, T + ∆] is given by + G k (t)dt) .(20)

  To get a first impression of the adequacy of the empirical model specification, we consider the MLFI model without any latent dynamics. In this case, the model has λ sk (t) = exp(η s ) and the MLFI model reduces to a simple homogeneous continuous-time Markov chain (HCTMC) model. The maximum likelihood (ML) estimator of η s for the HCTMC model has a closed-form expression and is given by

Figure 3 .

 3 Figure 3. Comparing this to the smoothed estimate for the random walk model in the lefthand panel (model D), we notice that the characteristics of the two paths in terms of peaks

Figure 5

 5 Figure5presents the log baseline hazard functions (cubic splines) capturing the semi-Markov

  discuss a number of metrics to compare the forecasts of transition matrices. In particular, they propose the average of singular values of the transition matrix minus the unit matrix. We computed this metric over time for the different model specifications used. The results are in the lower-right panel of Figure7.Again, we see a much larger resemblance between the dynamics of the AJ realizations with the MLFI specifications compared to the HMC model. It appears that the random walk models are the only ones able to explain the large increase in singular values in the early 2000s. As of realized volatility in the empirical finance literature, see[START_REF] Andersen | Modelling and forecasting volatility[END_REF]. The AJ estimator is its counterpart in the context of transition models.

  however, this comes at the cost of very prudent capital standards in other years.

  is possible to estimate and test formally for the number of latent factors driving the default and rating migration intensities. The implementation of the model for credit risk simulations is straightforward given the integrated structure of the model. The dynamics of the common risk factor are estimated simultaneously with the development of default events conditional on this common risk factor. Therefore, they can also be easily integrated in a forecasting exercise as the one presented in Section 6. . 1997a, Portfolio credit risk part I.Risk, 111-117. Wilson, T. 1997b, Portfolio credit risk, part II. Risk, 56-61.

  (t))], for k = 1, . . . , K with K the number of firms, s ∈ S = {AAA → AA, AAA → A, . . . , CCC → B, CCC → Df lt}, the 49 different possible rating transitions for the 7 rating classes AAA, . . . , CCC. The factor ψ is = ψ s (t i ) follows an AR(1) process ψ si = ρ ti-ti-1 s

Figure 1 :

 1 Figure 1: True versus smoothed estimate of ψ(t)

Figure 3 :

 3 Figure 3: Smoothed credit cycle ψ i with 95% confidence band for model D and G The left-hand panel contains the smoothed estimate of ψ i for a model with a single common random walk component (model D). The right-hand panel gives the smoothed estimate for a model with a single common AR(1) credit risk component (model G).

Figure 4 :

 4 Figure4: Smoothed credit cycles ψ s (t i ) with 95% confidence band for model I This figure contains the smoothed estimates of ψ s (t i ) for a model with 6 independent AR(1) systematic credit risk factors. The 6 factors describe the systematic credit risk in upgrades from investment grade to investment grade (upper-left: II up ), downgrades from investment grade to investment grade (upper-middle: II down ), downgrades from investment grade to sub-investment grade and default (upper-right: IS down ), upgrades from sub-investment grade to investment grade (lower-left: SI up ), upgrades from sub-investment grade to subinvestment grade (lower-middle: SS up ), and downgrades from sub-investment grade to sub-investment grade and default (lower-right: SS down ).

Figure 6 :

 6 Figure 6: Forecasted Investment Grade default probabilities Using a recursive state estimation procedure with extending data window from Dec 1980-June 1992 to Dec 1980-June 2005. Transition probabilities are estimated using the methodology of Section 4 for the MLFI model and the homogeneous continuous-time Markov chain (HMC) model. The MLFI specifications used from Table4are models D (1 common random walk (RW) factor), G (1 common autoregressive (AR) factor), E (6 independent RW factors), and H (6 independent AR factors). The Aalen-Johansen (AJ) estimates for each year are also plotted as a proxy for the observed transition rates. There are 4 plots for the rating categories AAA, AA, A, and BBB, respectively.

Figure 7 :

 7 Figure 7: Forecasted Sub-Investment Grade default probabilities

  The conditional mean of ν can be computed by the Kalman filter and smoothing (KFS) algorithm. More formally,

the KFS method computes E G (ν|Y ) where E G (•) is with respect to the approximating linear Gaussian model. It is well-known that the mode and the mean are equivalent in a Gaussian model. The new estimate of ν is denoted by ν (

Table 2 :

 2 Event counts

	This table contains the number of events over the sample period
	from (row) a specific rating to a new (column) rating. D denotes
	default.
	From				To
		AAA	AA	A BBB	BB	B CCC	D
	AAA		95	8	1	2	0	0	0
	AA	34		475	25	2	5	1	0
	A	10	219		714	40	19	0	9
	BBB	2	27	457		605	53	6	24
	BB	4	8	36	467		769	48	65
	B	0	6	16	34	524		855 349
	CCC	1	0	5	8	13	135		678

Table 3 :

 3 Parameter estimates of the HCTMC model This table presents estimates of a homogeneous continuous-time Markov chain (HCTMC) model with intensities λ sk (t) = exp(η s

Table 4 :

 4 Parameter estimates of the MLFI modelThis table contains the parameter estimates ρ s and αs = 100α s of the MLFI model, λ sk (t) = R sk (t) exp[η s + α s ψ s (t) + h s (u sk

We assume there are no simultaneous rating transitions. In practice the S&P's database is recorded at a daily frequency. This means multiple rating actions can be observed on a single day (for distinct firms). Our likelihood specification in Section 3 incorporates this phenomenon.

We assume N sk (t) to be a conditionally orderly process, i.e., it satisfiesP [N sk ((t + ∆) -) -N sk (t -) > 1 | F t -] = o (∆) P [N sk ((t + ∆) -) -N sk (t -) = 1 | F t -], such that we can discard the probability of a jump larger than 1 in N sk (t).

The possible endogenous nature of a selection of (time-varying) covariates leads to an inference procedure that can no longer be interpreted as full (conditional) maximum likelihood. Instead, we then have a partial likelihood inference framework, see[START_REF] Lancaster | The econometric analysis of transition data[END_REF].

If no exogenous covariates are included, as in the empirical illustration of Section 6, this is an innocuous assumption, seeTsiatis (1975, Theorem 2).

A discussion of the initial conditions problem in event-history models is provided by[START_REF] Van Den Berg | Duration Models: Specification, identification, and multiple durations[END_REF].

For a Markov chain, the entries of the generator matrix are either constants or (deterministic) functions of time. However, for generalized semi-Markov processes the entries of the generator matrix are, in general, stochastic processes.

for helpful comments. Financial support from Fundação para a Ciência e a Tecnologia (Portuguese Foundation for Science and Technology) and the Asset Management Program at Vrije Universiteit Amsterdam is gratefully acknowledged. The rating transition data for this research were generously supplied by Standard and Poor's.

A c c e p t e d m a n u s c r i p t Table 1: Monte Carlo Results

This table contains parameter estimates for the MLFI model, λ sk (t) = R sk (t) • exp[η s + α s ψ(t)], for k = 1, . . . , K with K the number of units, s = 1, . . . , 49. The true parameters are taken in accordance with the empirical results in Section 6, with different loading α s for up and down-grades, model C in Table 4. There are 7 rating classes. Initial ratings are distributed evenly over these classes. We estimate the model for K = 70 and K = 700, both with a random walk common factor ψ and with a stationary AR, parameter ρ = 0.9. The maximum time T is set to 25 years, unless the complete sample has defaulted at an earlier stage. We performed 500 replications for each parameter combination. Monte-Carlo averages and standard errors (in parentheses) are presented for those parameters that have a sufficient number of occurrences over all simulations.

True K=70 K=700 K=70 K=700 ρ true = 1 ρ true = 1 ρ true = 0.9 ρ true = 0.9