Areski Cousin
email: areski.cousin@univ-lyon1.fr

Mohamed Reda Kheliouen
email: mohamed-reda.kheliouen@etu.univ-lyon1.fr

A comparative study on the estimation of factor migration models

Keywords: Factor migration models, ordered Probit model, multi-state latent intensity model, mobility index, Kalman filter

, this study shows that short migrations of investment grade firms are significantly correlated to the business cycle whereas, because of lack of observations, it is not possible to state any relation between long migrations (more than two grades) and the business cycle. Concerning non investment grade firms, downgrade migrations are negatively related to business cycle whatever the amplitude of the migration.

Introduction

Analyzing the effect of business cycle on rating transition probabilities has been a subject of great interest these last fifteen years, particularly due to the increasing pressure coming from regulators for stress testing.

We consider the statistical estimation of credit migration models where, for each firm, the transition matrix is stochastic and may depend on a pool of common dynamic factors. The underlying factors aim at representing the evolution of the business cycle. The estimation procedure of such models differs given that the factors are considered to be observable or unobservable. In the first approach, one selects observable factors (macroeconomic variables for example) via economic analysis and then estimate the transition probabilities' sensitivities with respect to these covariates. [START_REF] Kavvathas | Estimating credit rating transition probabilities for corporate bonds[END_REF] calibrates a multi-state extension of a Cox proportional hazard model with respect to the 3-months and 10-years interest rates, the equity return and the equity return volatility. He finds that an increase in the interest rates, a lower equity return and a higher equity return volatility are associated with higher downgrade intensities. [START_REF] Nickell | Stability of rating transitions[END_REF] employ an ordered Probit model and prove the dependence of transition probabilities on the industry of the obligor, the business country of the obligor and the stage of the business cycle. [START_REF] Bangia | Ratings migration and the business cycle, with application to credit portfolio stress testing[END_REF] separate the economy into two regimes, expansion and contraction, estimate an ordered Probit model and show that the loss distribution of credit portfolios can differ greatly, as can the concomitant level of economic capital to be assigned. The second approach has emerged in response to criticisms made against the first approach. As [START_REF] Gourieroux | Risque de crédit : une approche avancée[END_REF] point out, the risk in selecting covariates lies in excluding other ones which would be more relevant. Specifications with latent covariates appear in articles such as [START_REF] Gagliardini | Stochastic migration models[END_REF] who consider an ordered Probit model with three unobservable factors and perform its estimation using a Kalman filter on ratings data of french corporates. They find that the two first factors are related to the change in french GDP (considered as a proxy of the business cycle). [START_REF] Koopman | The multi-state latent factor intensity model for credit rating transitions[END_REF] proceed on a parametric intensity model by conditioning the migration intensity on both observable factors and latent dynamic factors. The estimation shows the existence of a common risk factor for all migrations. 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.

The aim of this paper is to assess and compare two alternative stochastic migration models on their ability to link the transition probabilities to either observable or unobservable dynamic risk factors. In this respect, we use the same data set to compare the the multi-state latent factor intensity model used in [START_REF] Koopman | The multi-state latent factor intensity model for credit rating transitions[END_REF] and the ordered Probit model as described for instance in [START_REF] Gagliardini | Stochastic migration models[END_REF]. We use the S&P credit ratings history [01/2006 -01/2014] of a diversified portfolio composed of 2875 obligors and distributed across several regions and sectors. When the underlying factors are assumed to be unobservable, we adapt a method introduced in Gagliardini and [START_REF] Gagliardini | Stochastic migration models[END_REF] to represent the considered factor migration model as a linear Gaussian models. The unobservable factors are then filtered by a standard Kalman filter.

The paper is organized as follows. In Section 2, we present the class of factor migration models. Then, we introduce intensity models and structural models as particular factor migration models. Section 3 describes the estimation procedure for a multi-state factor intensity model in the two cases where the factor is assumed to be observable or unobservable. In Section 4, we present the estimation procedure associated with a structural ordered Probit model and we also carry the case of observable and unobserbable factors. In Section 5, we perform the estimation procedures on S&P credit ratings historical data and provide our main results and findings.

Factor migration models

In this paper, we consider Markovian model consisting of a multivariate factor process X and a vector R representing the rating migration process in a pool of n obligors. More specifically, for any time t ≥ 0, R t = (R 1 t , . . . , R n t ) is a vector in {1, • • • , d} n , where d corresponds to the default state and 1 corresponds to the state with the best credit quality. The l-th entry of R describes the rating migration dynamics of obligor l in the set S = {1, • • • , d}. The transition probabilities are driven by a factor process X. We consider that all sources of risk are defined with respect to the probability space (Ω, F, P) endowed with some reference filtration F = (F u ) 0≤u≤t satisfying the usual conditions. For any l = 1, . . . , n, let H l be the filtration associated with the migration process R l and let H be the global filtration such that

H l u = H 1 u ∨ • • • ∨ H n u , u ≥ 0.
Let G be the filtration associated with the factor process X. The credit migration models we consider in this paper are in the class of factor migration models.

Definition 2.1. A factor credit migration model is a Markov process (X, R) such that

• The factor X is a Markov process (in its own filtration).

• Given the history of X, i.e. given G ∞ , the marginal rating migration processes R 1 , . . . , R n are independent (time-inhomogeneous) Markov chains with the same transition matrices.

Contrary to the definition of stochastic migration models in [START_REF] Gagliardini | Stochastic migration models[END_REF], here the factor process X is not necessarily identified as a common stochastic transition matrix. Note that the factor process X plays two important roles : it introduces dynamic dependence between the obligor rating migration processes and it allows for non-Markovian serial dependence in the migration dynamics: the rating migration process R is not a Markov process. Note also that the conditional Markov chain (R 1 , . . . , R n ) has d n transition states. The conditional independence property allows to significantly reduce the problem dimension when it turns to compute the joint distribution of rating migration events. In that sense, conditionally to the filtration G ∞ , the R 1 t , . . . , R n t are iid for every date t. Moreover, there is no contagion mechanism in this framework. A migration event in the pool has no impact on the obligors migration probabilities since the latter are driven by a process X which is Markov in its own filtration. In [START_REF] Koopman | The multi-state latent factor intensity model for credit rating transitions[END_REF], the factor process X can contain obligor-specific informations (microeconomic variables) as well as common observable factors (macroeconomic variables) or common unobserved latent factors. The model presented in this paper does not distinguish between obligors as the process X stands only for common observable or latent factors.

Intensity models

In this section, we define multi-state factor intensity models as particular factor migration models. Given the history of X (given G ∞ ), the rating migration processes (R l t ) t≥0 , l = 1, . . . , n are conditionally independent continuous-time Markov chains. Moreover, for any time t, they are assumed to have a common generator matrix Λ X (t) defined by

Λ X (t) =         -λ 1 (X t ) λ 12 (X t ) • λ 1,d-1 (X t ) λ 1,d (X t ) λ 21 (X t ) -λ 2 (X t ) • λ 2,d-1 (X t ) λ 2,d (X t ) • • • • • • • • • • λ d-1,1 (X t ) λ d-1,2 (X t ) • -λ d-1 (X t ) λ d-1,d (X t ) 0 0 • 0 0        
, where, for any j = i, λ ij are positive functions and

λ i := j =i λ ij , for i = 1, • • • , d -1.
For j = i, the product λ ij (X t )dt corresponds to the conditional probability of going to rating state j in the small time intervalle (t, t + dt], given G t and the fact that the obligor is in rating state i at time t. For small time length dt, the following first-order approximation holds

λ ij (X t ) dt ≈ P R l t+dt = j | R l t = i, G t , (1)
for any obligor l = 1, . . . , n. In turns out that λ i (X t ) dt is the conditional probability to depart from state i in the small interval (t, t + dt] given G t .

A classical specification of migration intensities is an exponential-affine transformation of the common factor. More specifically, for any i = j,

λ ij (X t ) = exp (α ij + β ij , X t ) , (2)
where α ij is a constant parameter and β ij accounts for the sensitivity of migration intensity λ ij to the common factor X. This specification garanties the positivity of migration intensities. This specification corresponds to a multi-state extension of the Cox proportional hazard model.

In this framework, the conditional probability transition matrix P (given G ∞ ), can be computed by solving the forward Kolmogorov equation

∂P (t, s) ∂s = P (t, s)Λ X (s), P (t, t) = Id. ( 3
)
The unconditional transition matrix is then given by

Π(t, s) = E [P (t, s)] , (4)
where the expectation is taken over the distribution of (X u ) t<u≤s . In practice, the transition matrix Π(t, s) can be approximated by Monte Carlo simulations, based on independent simulations of the path of X between t and s. For each realization of (X u ) t<u≤s , a simple numerical scheme can be used to solve (3).

If the factor process X only changes at discrete times t 1 < . . . < t N , the conditional transition matrix can be expressed as a product involving the following matrix exponential terms

P (t k , t k+1 ) = e Λ X (t k )(t k+1 -t k ) .
(5)

Then, under this assumption, the forward Kolmogorov equation ( 3) has an explicit solution. Note that, if the generator matrix Λ X is time-inhomogeneous, it is generally not true that one can extend (5) to get

P (t, s) = e s t Λ X (u) du . (6)
However, this relation holds for some specification of Λ X for which the matrix Λ X (u) and Λ X (u ) for all u, u in [t, s].

Structural models

We consider a discrete-time structural models where any firm l jumps to a new rating category when a quantitative latent process S l crosses some pre-specified levels or barriers. In the classical structural Merton model, S l is defined as the ratio of asset value and liabilities. The rating of name l at time t is given by the position of the latent variable S l t inside a pre-specified partition of the real line

-∞ = C d+1 < C d < • • • < C i+1 < C i < • • • < C 1 = +∞. More formally, R l t = d i=1 i1 {Ci+1≤S l t <Ci} . (7)
Several models exist in literature for the specification of S l (see, e.g., [START_REF] Nickell | Stability of rating transitions[END_REF], [START_REF] Bangia | Ratings migration and the business cycle, with application to credit portfolio stress testing[END_REF], [START_REF] Albanese | Credit barrier models[END_REF], [START_REF] Feng | The ordered qualitative model for credit rating transition[END_REF]). In this study, we choose [START_REF] Gagliardini | Stochastic migration models[END_REF]'s approach since it is a generalization of previously cited models, it is also investigated in [START_REF] Gourieroux | Risque de crédit : une approche avancée[END_REF] and [START_REF] Feng | The ordered qualitative model for credit rating transition[END_REF]. For any obligor l = 1, . . . , n and any time t, the latent process S l t is expressed as a deterministic affine transformation of a common factor X t and of an independent idiosyncratic factor ε l t . The characteristics of this affine transformation may depend on the rating state at the preceding date t -1. The latent process S l is then described by

S l t = d i=1 α i + β i , X t + σ i ε l t 1 {R l t-1 =i} , (8)
where α i is the level of S l at rating i, β i represents the sensitivity in rating i of the latent factor S l to the common factor X, σ i corresponds to the volatility of residuals and ε l t are iid random variables independent of X. When the idiosyncratic residual processes ε l , l = 1, . . . , n are not specified, this corresponds to an ordered polytomous model. The most common version is the ordered Probit model where ε l t are independent standard Gaussian variables (see [START_REF] Gagliardini | Stochastic migration models[END_REF] or [START_REF] Feng | The ordered qualitative model for credit rating transition[END_REF] for more details). In what follows, we consider the ordered Probit model and we denote by θ = (C i+1 , α i , β i , σ i ) i=1,...,d-1 the set of unknown parameters.

In this framework, the conditional transition probabilities p ij are given by1 ,

p ij = P R l t = j | R l t-1 = i, G ∞ , = P C j+1 ≤ S l t < C j | R l t-1 = i, G ∞ , = P C j+1 ≤ α i + β i , X t + σ i ε l t < C j | R l t-1 = i, X t , = P C j+1 -α i -β i , X t σ i ≤ ε l t < C j -α i -β i , X t σ i | X t .
Then, if Φ(•) is the cumulative distribution function of a standard Gaussian variable, we obtain, for

any i = 1, • • • , d -1, p ij = Φ C j -α i -β i , X t σ i -Φ C j+1 -α i -β i , X t σ i , j = 2, • • • , d -1, (9)
p i1 = 1 -Φ C 2 -α i -β i , X t σ i , (10)
p id = Φ C d -α i -β i , X t σ i . ( 11
)
3 Statistical estimation of intensity models

In this section, we consider a multi-state intensity model as described in Section 2.1 and we explain how to estimate the parameters of the generator matrix given the sample history of the obligors rating migrations. In a general setting, we first give the expression of the conditional likelihood function. For the ease of the presentation, we distinguish the simple case where migration intensities are constant over time and the case where, as in a (multi-state) proportional hazard model, the intensities are given as exponential-affine transformation of an observable factor X. In this setting, the estimation of model parameters is made by maximizing the conditional likelihood function. When the factor X is assumed to be unobservable, the estimation requires the computation of the unconditional likelihood function. This task may be computationally intensive as explained in [START_REF] Koopman | The multi-state latent factor intensity model for credit rating transitions[END_REF]. In this paper, we choose to adapt the approach of [START_REF] Gagliardini | Stochastic migration models[END_REF] (estimation of an ordered Probit model with unobervable factors) to the multi-state factor intensity model.

We assume that, for any obligor l = 1, • • • , n, the observed number of ratings visited during the period [0, t) is denoted by N l (N l ≥ 1). For any k = 1, • • • , N l , the time intervals t l k-1 , t l k correspond to the visiting of new state r l k where t l 0 = 0 and t l N l = t. Then, the observed path of ratings of obligor l during the period [0, t) is described by

r l u = N l k=1 r l k 1 {t l k-1 ≤u<t l k } , (12)
There are N l -1 migration events observed for obligor l in the time interval [0, t), each of them has took place at time

t l k for k = 1, • • • , N l -1.

Conditional likelihood function

Proposition 3.1. Let θ be the set of parameters which characterizes the functional link between the generator matrix Λ X and the factor process X. The conditional likelihood function given the observed path of rating migration histories r l u 0≤u<t and the path of the risk factor (X u ) 0≤u≤t can be expressed as

L(θ | G t ) = n l=1 N l k=1 λ r l k ,r l k+1 (X t l k ) e - t l k t l k-1 λ r l k (Xu) du , (13)
with the convention λ r l N l

,r l N l +1 = 1. Proof. See proof in appendix A.

Note that the conditional likelihood function can be expressed as a product of (conditional) marginal likelihood function, that is

L(θ | G t ) = i =j L ij (θ | G t ) (14)
where the product is taken over all transition types and L ij (θ | G t ) is given in the following proposition.

Proposition 3.2. The marginal likelihood function L ij (θ | G t ) associated with migration from state i to state j (i = j) is given by,

L ij (θ | G t ) = n l=1 N l k=1 exp Y l ij (t l k ) log λ ij (X t l k ) -S l i (t l k ) t l k t l k-1 λ ij (X u ) du (15)
where

Y l ij (t l k ) = 1 if i = r l t l k and j = r l t l k+1 , 0 else, (16)
for k = 1, . . . , N l -1, Y l ij (t l k ) = 0 for k = N l and S l i (t l k ) = 1 if i = r l t l k 0 else, (17)
for k = 1, . . . , N l .

According to (12), the quantity Y l ij (t l k ) is equal to one (and zero otherwise) if the couple of indices (i, j) corresponds to the observed migration event of obligor l at time t l k , i.e., if i = r l t l k and j = r l t l k+1 .

Moreover, Y l ij (t l k ) = 0 for k = N l since, by construction, no migration occurs after time t l N l -1 . The quantity S l i (t l k ) is equal to one (and zero otherwise) if index i corresponds to the rating state visited by obligor l just before migration date t l k , i.e., if i = r l t l k .

Proof. Proposition (3.2) is a direct consequence of Proposition (3.1).

If each migration intensity function λ ij is described by specific parameters, the maximization of the conditional likelihood function (13) can be done by maximizing each marginal likelihood independently with respect to its own set of parameters. Numerically, it is more convenient to work with the loglikelihood function, we then transform (15) to

log (L ij (θ | G t )) = n l=1 N l k=1 Y l ij (t l k ) log λ ij (X t l k ) - n l=1 N l k=1 S l i (t l k ) t l k t l k-1 λ ij (X u ) du, (18)
Likelihood function in the multi-state intensity model is classical from the point of view of point processes [START_REF] Hougaard | Analysis of Multivaiate Survival Data[END_REF]). We compute here the likelihood function for the conditional Markov chain assumption. The link between the two approaches is clear as soon as we recall that Markov chains are particular marked point processes. The likelihood expression ( 18) is the same one obtained in [START_REF] Koopman | The multi-state latent factor intensity model for credit rating transitions[END_REF] except that in [START_REF] Koopman | The multi-state latent factor intensity model for credit rating transitions[END_REF] the time is discrete and divided according to the observed dates of change in value of the common process X.

Time-homogeneous intensities

Let us assume that the transition intensity functions λ ij are constant, so that, λ ij (X t ) = λ ij,0 for any time t. As explained above, the max-likelihood estimation of λ ij,0 , i = j is obtained by maximizing each corresponding marginal likelihood functions.

Proposition 3.3. The maximum likelihood estimate of the transition intensities is given by,

λij,0 = n l=1 N l k=1 Y l ij (t l k ) n l=1 N l k=1 S l i (t l k ) t l k -t l k-1 . ( 19
)
Note that the estimator λij,0 is expressed as the ratio of the total number of observed transitions from state i to state j over the cumulated time spent by the obligors in state i.

Proof. Given (18), λij,0 is the solution of ∂ log (L ij (θ | G t ))
∂λ ij,0 = 0. And, we have

∂ log (L ij (θ | G t )) ∂λ ij,0 = n l=1 N l k=1 Y l ij (t l k ) ∂ log (λ ij ) ∂λ ij,0 - n l=1 N l k=1 S l i (t l k ) ∂λ ij ∂λ ij,0 t l k -t l k-1 = 1 λ ij,0 n l=1 N l k=1 Y l ij - n l=1 N l k=1 S l i (t l k ) t l k -t l k-1 = 0
Moreover, the marginal likelihood function is concave since it is easily seen that

∂ 2 log (L ij (θ | G t )) ∂ 2 λ ij,0 = n l=1 N l k=1 -Y l ij (t l k ) (λ ij,0 ) 2 < 0.
Under this specification, the dynamics of rating migrations does not depend on the business cycle. In the numerical part (see Section 5), this estimation procedure will be used to construct Through The Cycle2 (TTC) transition matrix.

Intensities depending on observable factors

We now assume that migration intensities depend on a K-dimensional observable factor X = (X 1 , . . . , X K ) through an exponential-affine relation:

λ ij (X t ) = λ ij,0 exp ( β ij , X t ) , i = j, i < d. ( 20
)
for all. The vector

β ij = β 1 ij , • • • , β K ij contains
the sensitivities of λ ij with respect to each component of factor X. The term β ij , X t denotes the inner product between vectors β ij and X t . The baseline intensity λ ij,0 is assumed to be constant. This specification corresponds to a multi-state version of the so-called Cox proportional hazards regression model, where here the baseline function is assumed to be constant. This setting has been investigated by among others [START_REF] Kavvathas | Estimating credit rating transition probabilities for corporate bonds[END_REF], [START_REF] Lando | Analyzing rating transitions and rating drift with continuous observations[END_REF], [START_REF] Koopman | The multi-state latent factor intensity model for credit rating transitions[END_REF], [START_REF] Naldi | Quantitative credit quarterly : Transition matrix estimation[END_REF].

Proposition 3.4. Under specification (20), for any transition type (i, j), the maximum likelihood

estimate βij = β1 ij , • • • , βK ij
is solution to the following non-linear system:

n l=1 N l k=1 Y l ij (t l k ) X s t l k n l=1 N l k=1 Y l ij (t l k ) = n l=1 N l k=1 S l i (t l k ) t l k t l k-1 X s,u e βij ,Xu du n l=1 N l k=1 S l i (t l k ) t l k t l k-1 e βij ,Xu du , s = 1, • • • , K. (21)
The maximum likelihood estimate λij,0 is given by

λij,0 = n l=1 N l k=1 Y l ij (t l k ) n l=1 N l k=1 S l i (t l k ) t l k t l k-1 e βij ,Xu du . ( 22
)
where βij is the solution of system (21).

Proof. See proof in appendix B.

We notice from ( 22) that, when βij is fixed to 0, the maximum likelihood estimate of the baseline intensity falls back to the one of the homogeneous case (see ( 19)).

Intensities depending on unobservable factors

In the previous setting, we implicitly assume that the dynamics of migration intensities are fully explained by pre-identified observable factors. However, it may be the case that the dynamics of the intensities are driven by other (un-identified) factors or even by totally unobservable factors. We consider here that the migration intensities follow specification (20) but the vector of factors X is now assumed to be unobservable. [START_REF] Koopman | The multi-state latent factor intensity model for credit rating transitions[END_REF] use a similar latent factor specification. We assume that the K-dimensional dynamic factor X only changes at times t 1 , • • • , t N and that its dynamics follows an auto-regressive (AR) process

X t = AX t-1 + ζ t , t = t 1 , . . . , t N (23)
where the matrix A characterizes the auto-regression coefficients and ζ t are iid and standard Gaussian variables. The maximum likelihood estimation of λ ij,0 and β ij involves the computation of the unconditional likelihood function L(θ) given by,

L(θ) = E [L (θ | G t )] , (24)
where expectation in ( 24) is taken over the joint distribution of (X t1 , • • • , X t N ). Consequently, maximizing the likelihood function is computationally very intensive3 (see [START_REF] Koopman | The multi-state latent factor intensity model for credit rating transitions[END_REF]) as soon as N is larger than a few units. We choose to follow another route by adapting the approach of [START_REF] Gagliardini | Stochastic migration models[END_REF] (initially proposed for a structural model) to the multi-state factor intensity model. The idea is to construct an approximation of model ( 20) as a linear Gaussian model which can be dealt with a Kalman filter. As usual, the estimation of model parameters is made by maximizing the likelihood of the filtered model.

Representation as a linear Gaussian model

Let us first remark that equation ( 20) is equivalent to the following linear relation

y ij,t = α ij + β ij , X t , t = t 1 , . . . , t N (25)
where y ij,t = log (λ ij,t ) and α ij = log (λ ij,0 ). Note that, even if the transition intensities λ ij,t , t = t 1 , . . . , t N (and then y ij,t ) are not directly observed, they can be estimated from the panel data. We assume that, for any time indices k = 1, . . . , N and for any transition type (i, j), λ ij,t k can be estimated from the migration dynamics observed in time interval (t k-1 , t k ) using the max-likelihood estimate (19). We then construct a time series of estimated migration intensities λij,t and log-intensities ŷij,t := log( λij,t ) t = t 1 , . . . , t N . The asymptotic normality of the maximum likelihood estimate ( 19)

writes √ n λij,t -λ ij,t → N (0, σ 2 ij,t
) as the number n of obligors goes to infinity. The asymptotic variance σ 2 ij,t can be approximated easily4 (see [START_REF] Hougaard | Analysis of Multivaiate Survival Data[END_REF]). Following the Delta method, we also know that

√ n (ŷ ij,t -y ij,t ) → N 0, χ 2 ij,t where χ 2 ij,t = σ 2 ij,t (λ ij,t ) 2 .
As a result, if the panel data is sufficiently large, by using (25), we obtain the following linear Gaussian system

ŷij,t α ij + β ij , X t + χij,t κ ij,t , ∀i = j, t = t 1 , . . . , t N X t = AX t-1 + ζ t , (26)
where for any i = j, χij,t := χij,t √ n , κ ij,t and ζ t are independent error terms, distributed as standard Gaussian variables. For any t = t 1 , . . . , t N , we denote by ŷt the column vector composed of the elements ŷij,t , i = j, i < d. Let p be the number of rows of ŷt and Y t = (ŷ t1 , . . . , ŷt ) be the information collected up to time t. Using the Kalman filter, the unobserved factor X can be filtered recursively given the estimated process ŷ. For any t = t 1 , . . . , t N , let us denote by Xt := E(X t | Y t ) the filtered version of the factor process and by Xt|t-1 := E(X t | Y t-1 ) the best prediction of X t given information up to time t -1. In the Kalman filter terminology, the first equation in ( 26) corresponds to the measurement equation, the second to the transition equation.

Maximum likelihood of the filtered model

Let θ = (α ij , β ij ) i =j,i<d be the vector of unknown model parameters. The log-likelihood function associated with the filtered Gaussian model is given by (see, e.g., [START_REF] Durbin | Time Series Analysis by State Space Methods, Second Edition[END_REF]),

log (L(θ, A)) = - N p 2 log (2π) - 1 2 N t=1 log |F t | + e T t F -1 t e t , (27)
For any t = t 1 , . . . , t N , the column vector e t is the prediction error and is equal to

e t := ŷt -E[ŷ t | Y t-1 ].
The matrix F t is the conditional variance of e t given Y t-1 , i.e. F t = Var(e t | Y t-1 ), and |F t | is its determinant. Given the vector of "observations" Y t N , the prediction error e t and the matrix F t can be obtained at each time t as outputs of the Kalman routine.

The parameters αij , βij are chosen as the ones that maximize the log-likelihood function ( 27). We then define the estimated migration intensities as

λij ( Xt ) = exp αij + βij , Xt , (28)
where X the filtered factor process associated with max-likelihood parameters θ, Â. The relation (28) will be used in the application on S&P data (see Subsection 5.4).

Statistical estimation of structural models

In this section, we consider the structural model described in Section 2.2 and we explain how to estimate the parameters of the conditional transition probability p ij (see (9-11)) given the sample history of the obligors rating migrations. As for the intensity model, we distinguish the case where the underlying factor is observable and the case where it is unobservable. When the factor X is observable, the estimation of model parameters is made by maximizing the conditional likelihood function. When the factor X is assumed to be unobservable, the estimation requires the computation of the unconditional likelihood function. This task can be computationally very intensive. As for the intensity model, we consider as an alternative route the estimation approach used in [START_REF] Gagliardini | Stochastic migration models[END_REF].

Observable factors

Assuming that the conditional transition probability p ij depends on a K-dimensional observable factor X = (X 1 , . . . , X K ) and considering θ = (C i+1 , α i , β i , σ i ) i=1,••• ,d-1 as the set of parameters which characterize the functional link between p ij and the factor process X. For any k = 1, • • • , N , the time intervals [t k-1 , t k ) correspond to the visiting of the state r l k where t 0 = 0 and t N = t. The observed path of ratings of obligor l during the period [0, t) is described by

r l t = N k=0 r l k 1 {t=t k } , (29)
Proposition 4.1. The conditional likelihood function associated with the observed path of rating migration histories r l t , t = t 0 , . . . , t N , l = 1, . . . , n and given the initial ratings and the path of the risk factor (X t ), t = t 1 , . . . , t N can be expressed as,

L(θ | G t N ∨ H t0 ) = n l=1 N k=1 i =j (p ij (X t k )) Y l ij (t k ) , (30)
where, for k = 1, . . . , N ,

Y l ij (t k ) = 1 if i = r l t k-1 and j = r l t k , 0 else. Proof. L(θ | G t N ∨ H t0 ) = n l=1 P R l u = r l u , u = t 1 , . . . , t N | G t N , R l t0 = r l t0 = n l=1 N k=1 P R l t k = r l t k | R l t k-1 = r l t k-1 , G t N p ij (X t k
) is defined by ( 9), ( 10) and ( 11)

and θ = (C i+1 , α i , β i , σ i ) i=1,••• ,d-1 .
We deduce that,

θ = arg max θ    n l=1 N k=1 i =j Y l ij (t k ) log (p ij (X t k ))    , (31)
solving ( 31) is done numerically by considering constraints on the positivity of σ i > 0 and the or-

der of threshold C i (-∞ = C d+1 < C d < • • • < C 2 < C 1 = +∞) (see appendix C to have more details).

Unobservable factors

We now consider that the conditional transition probability p ij depends on a K-dimensional unobservable factor X = (X 1 , . . . , X K ). The dynamic factor X is assumed to change at times t 1 , • • • , t N with the following auto-regressive dynamics,

X t = AX t-1 + η t , (32)
where the matrix A characterizes the covariates dynamics and η t are iid and standard Gaussian.

The maximum likelihood estimation of θ = (C i+1 , α i , β i , σ i ) i=1,••• ,d-1 involves the computation of the unconditional likelihood function L(θ) given by,

L(θ) = E   n l=1 N k=1 j =i (p ij (X t k )) Y l ij (t k )   . ( 33
)
Estimating θ using (33) raises the issue of the potentially high number of integrals to compute. Indeed, the distribution must be integrated with respect to covariates values

X t1 , • • • , X t N , which represents [K × N ] integrals.
Similarly to the intensity model, we use the method presented in [START_REF] Gagliardini | Stochastic migration models[END_REF] which consists in transforming the probability of migrating towards a rating J or below, denoted p * iJ,t (see ( 34)), as a linear function of latent covariates (see (36) below) and thus construct a linear Gaussian model.

Representation as a linear Gaussian model

First, we formally define p * iJ,t for t = t 1 , . . . , t N as,

p * iJ,t k = J j=d p ij,t k = P R l t k < J|R l t k-1 = i, X t k , (34)
Note that p * i1,t k = 1 and thus the associated threshold C 1 = +∞. The thresholds to be estimated are

C J for J = 2, • • • , d, or equivalently C i+1 for i = 1, • • • , d -1. We use the formulation C i+1 instead of C J in the parameters vector θ = (C i+1 , α i , β i , σ i ) i=1,••• ,d-1
for the ease of notation.

Assuming an ordered Probit model, (34) becomes,

p * iJ,t = Φ C J -α i -β i , X t σ i , ⇒ Φ -1 p * iJ,t = C J -α i σ i - 1 σ i β i , X t , (35)
or equivalently,

ψ iJ,t = C J -α i σ i - 1 σ i β i , X t , (36)
Note that p * iJ,t and ψ iJ,t , t = t 1 , . . . , t N are not observed but can be estimated from the panel data. We assume that, for any time indices k = 1, . . . , N and for any arrival state J, p * iJ,t (and then ψ iJ,t ) can be estimated from the migration dynamics observed in time interval (t k-1 , t k ) using the max-likelihood estimate (19) and the exponential matrix term expressed in (5). We then construct a time series of estimated p * iJ,t and ψiJ,t , t = t 1 , . . . , t N . The asymptotic normality of the estimator p * iJ,t writes √ n p * iJ,t -p * iJ,t → N (0, Ω 2 iJ,t ) as the number n of obligors goes to infinity5 . Following the Delta method, we also know that √ n ψiJ,t -ψ iJ,t → N 0, Ξ 2 iJ,t where Ξ iJ,t = Φ -1 (p * iJ,t )Ω iJ,t . As a result, if the panel data is sufficiently large, by using (36), we obtain the following linear Gaussian representation

   ψiJ,t C J -α i σ i - 1 σ i β i , X t + ΩiJ,t υ iJ,t , ∀i, J, t = t 1 , . . . , t N , X t = AX t-1 + η t , (37)
where ΩiJ,t :=

Ω iJ,t √
n , υ iJ,t and η t are independent error terms, distributed as standard Gaussian variables. For any t = t 1 , . . . , t N , we denote by ψt the column vector composed of the elements ψiJ,t , i < d. Let p be the number of rows of ψt and Ψ t = ( ψt1 , . . . , ψt ) be the information collected up to time t. Using the Kalman filter, the unobserved factor X can be filtered recursively given the estimated process ψt . For any t = t 1 , . . . , t N , let us denote Xt := E(X t | Ψ t ) the filtered version of the factor process and by Xt|t-1 := E(X t | Ψ t-1 ) the best prediction of X t given information up to time t -1. In the Kalman filter terminology, the first equation in (37) corresponds to the measurement equation, the second to the transition equation.

Maximum likelihood of the filtered model

Let θ = (C i+1 , α i , β i , σ i ) i=1,••• ,d-1
be the vector of unknown model parameters. The log-likelihood function associated with the filtered Gaussian model is given by [START_REF] Durbin | Time Series Analysis by State Space Methods, Second Edition[END_REF] as,

log (L(θ, A)) = - N p 2 log (2π) - 1 2 N t=0 log |F t | + e T t F -1 t e t , (38)
For any t = t 1 , . . . , t N , the column vector e t is the prediction error and is equal to

e t := ψt -E[ ψt | Ψ t-1 ].
The matrix F t is the conditional variance of e t given Ψ t-1 , i.e. F t = Var(e t | Ψ t-1 ), and |F t | is its determinant. Given the vector of "observations" Ψ t N , the prediction error e t and the matrix F t can be obtained at each time t as outputs of the Kalman routine.

The parameters θ = Ĉi+1 , αi , βi , σi

i=1,••• ,d-1
are chosen as the one that maximize the loglikelihood function ( 38) of the filtered model. We then define the estimated probability of migrating towards a rating J or below as

ψ iJ,t ( Xt ) = C J -α i σ i - 1 σ i β i , Xt , (39)
where X the filtered factor process associated with max-likelihood parameters θ, Â. The relation (39) will be used in the application on S&P data (see Subsection 5.4).

Empirical study

In order to perform the estimation procedures developed in sections 3 and 4, we have built a worldwide bonds portfolio, composed with a mixture of investment and non investment grade positions. The diversification of the portfolio in terms of sectors and regions will allow to compare between the two estimation approaches; estimation with observable factors and estimation with latent factors.

In this section, we describe the data and explain how we proceed to compute the transition and generator matrices in time homogeneous Markov chain (no covariates). In the second part, i.e. nonhomogeneous Markov chain (with observable or latent covariates), we perform the parameters estimation on both the intensity and the structural approaches. The third part of the application deals with the comparison between models and assesses the ability of both models to adjust the migration dynamics in terms of mobility (see section 5.5.2 for more details on mobility).

Data description

Our data stands for S&P credit ratings history, it covers the period January 2006 to January 2014 with a monthly frequency observation on a portfolio of 2875 obligors. The database has a total of 16168 obligor years excluding withdrawn ratings. The overall shares of the most dominant regions in the dataset, i.e. North America, Western Europe and Asia are respectively 48.4%, 24.5% and 17.2% (see Table (1)). The split by sector is less concentrated as the top five dominant sectors are Finance, Energy, Industry, Utilities and Telecom with respectively 17.5%, 9.1%, 7%, 6% and 5.1% of the total shares (see Table (2)). We notice also that the portfolio is dynamic, as the composition is changing through time but the overall shares per rating remain stable (see Table (3)) .

The S&P ratings contain 8 classes of risk, AAA (lowest risk), AA, A, BBB, BB, B, CCC and D which stands for Default. For convenience we replace these qualitative notations with numeric equivalent, as 1 denote the AAA, 2 the AA,... until 8 for D. .9 8.0 28.5 31.9 13.4 12.5 1.9 0.9 100 2011 3.4 7.5 28.2 32.5 13.8 13.3 1.0 0.3 100 2012 2.4 6.5 26.4 32.6 16.5 13.8 1.4 0.4 100 2013 2.5 6.2 25.7 34.0 15.9 14.1 1.4 0.2 100 Table 3: Percentage of obligors per ratings and per years

Time-homogeneous intensity model (no covariates)

In the case of homogeneous Markov chain, the transition probabilities are time-invariant. This means that both the cohort approach (see [START_REF] Löffler | Credit Risk Modelling using VBA and Excel[END_REF]) and the intensity approach with a constant generator matrix (see (19) in section 3.2) can be used to compute the transition matrix. We have investigated both of them. The cohort approach turns out to produce a high number of null transition probabilities, which makes the calculation of quantities like Φ -1 p * ij impossible 6 . This is not the case of the intensity approach where almost all transition probabilities were not null 7 . For instance, using a constant generator matrix on time interval [01/2013 -01/2014] leads to the following results on the 2013 one year transition matrix 8 , see Table (4 7 This difference is due to the fact that the cohort approach does not make full use of the available data. Specifically, the estimates of the cohort approach are not affected by the timing and sequencing of transitions within the period where the intensity approach captures within-period transitions.

8 The probabilities lower than 0.001 are reported as null in the matrix.

Table 4: One year transition matrix on 2013

The transition matrix in Table ( 4) is typical of the migration dynamic behaviour. The probabilities on the diagonal are close to one which means that the ratings are quite stable during a short time period (which is the case for 1 year). The probabilities around the main diagonal (above and below) are also significant, while the other ones tend generally to zero, indeed rating migrations are usually by one grade. However, for the non investment grades, i.e. ratings 5, 6 and 7 it is common to have more heterogeneity than the investment grades. This is due in one hand to the financial situation of the firms standing at these ratings and in the other hand to the correction effect of agencies. When a default which has not been predicted by the rating agencies occurs, several adjustments are performed to correct for this prediction error, some downgrades are then applied to the firms which are in the same situation as the defaulted one.

If we consider a constant matrix generator on the overall study period, i.e. [01/2006 -01/2014] we get the following transition matrix that we can consider as a Through The Cycle transition matrix(see Table (5 Comparing to the 2013 matrix, the TTC matrix is less concentrated on the diagonal and almost all the transition probabilities are not null. This is due to computation with respect to all the duration of data history, all the migrations are likely to happen during this time.

Estimation with observable factors

Many studies have addressed the question of selecting the appropriate covariates to explain time variation in the behaviour of ratings. This question is a large topic of investigation, one can refer to [START_REF] Kavvathas | Estimating credit rating transition probabilities for corporate bonds[END_REF], [START_REF] Bangia | Ratings migration and the business cycle, with application to credit portfolio stress testing[END_REF], [START_REF] Couderc | Time-to-default: Life cycle, global and industry cycle impacts[END_REF] and [START_REF] Koopman | Credit cycles and macro fundamentals[END_REF] to have more details. In this study, we have followed the idea of [START_REF] Couderc | Time-to-default: Life cycle, global and industry cycle impacts[END_REF] by selecting variables according to three covariates "families", i.e. macroeconomic information, financial and commodity markets information and credit market information. As our portfolio is mainly composed of US and Euro positions we focused on US and Euro indicators. The entire list of covariates is given below9 .

Macroeconomic variables : US and Euro Real Gross Domestic Product, US and Euro Consumer Price Indices, US and Euro Civilian Unemployment rates, US Effective Federal Funds rate (short term interest rate), Euro short term interest rates, US 10 years treasury constant maturity rate, longterm Government bond yields 10-years, US and Euro Industrial Production Indices, US and Euro Purchasing Manager Indices, Case-Shiller Index (real estate price Index for US).

Financial and commodity markets variables : Crude Oil price, Gold price, NASDAQ 100 Index, Dow Jones Index, Euro STOXX Index and Chicago Board Options Exchange Market Volatility Index.

Credit market variables : Bank of America Merrill Lynch US Corp BBB total return Index, Bank of America Merrill Lynch US Corp AA total return Index, Markit iBoxx AA and Markit iBoxx BBB.

The data time series go from 01/2006 to 01/2014 which encompasses our study period. Each variable consists in a year to year growth without overlapping between the sub-periods, this gives 8 sub-periods. For example, the US Real GDP growth for the date 01/2014 (last sub-period) is calculated as the relative variation during the period [31/12/2012 -31/12/2013] The first factor explains almost 44% of the total variance, where the rest of factors (2 until 8) explain 50% cumulated variance. As there is not a strict rule to choose the number of covariates, we have investigated the estimation procedure for different number of covariates. Above five factors the log-likelihood is no more significantly improved, we have then chosen the five first factors with a cumulated explained variance of 86% to realise the estimation procedure. The number of selected factors can be challenged in the light of a deeper investigation, but it is not the purpose of this paper.

Parameters estimation: we recall that for the intensity model, the parameters to be estimated are the baseline intensity Λ 0 and the matrices of factor sensitivities β s , s = 1, • • • , 5. We use proposition 3.4 to estimate these parameters from the panel data composed of 8 sub-periods. For the structural model, the parameters are θ

= (C i+1 , α i , β i , σ i ) i=1,••• ,7 (see (9 -11)).
The estimated parameters are reported in the following tables. See Tables (7), ( 8) below and in the appendix for the intensity model and Tables ( 9), ( 10) for the structural model. The migrations are expected to be correlated with the business cycle, negatively when the migration is a downgrade and positively when the migration is an upgrade. According the Table (8), where the first factor can be considered as a proxy of the business cycle, this behaviour is confirmed for the short migrations (i.e. one grade migration, for instance β12 = -0.06, β23 = -0.11, β32 = 0.22, β21 = 0.34). When the migration is farther (i.e. more than two grades), due to the lack of observations, it is no more possible to state any significant relation with the business cycle. When we focus on the non investment grade obligors (ratings 5, 6 and 7), we can observe that, all downgrades are negatively correlated to the business cycle, whatever the amplitude of the migration.

Ĉ

α σ Ĉ2 = 5.16 α1 = 5.08 σ1 = 0.17 Ĉ3 = 4.35 α2 =3.83 σ2 = 0.38 Ĉ4 = 3.47 α3 = 3.14 σ3 = 0.74 Ĉ5 = 2.16 α4 = 2.47 σ4 = 1.18 Ĉ6 = 1.05 α5 = 1.52 σ5 = 1.17 Ĉ7 = -0.19 α6 =0.63 σ6 = 1.17 Ĉ8 = -1.55 α7 =-1.01 σ7 = 1.72 The results of the Tables ( 9) and ( 10) are consistent with those obtained by [START_REF] Gagliardini | Stochastic migration models[END_REF] as the estimated barrier levels Ĉ increase with respect to rating quality. The same conclusion holds for the intercepts αi which confirms that the downgrade risk is higher for low quality credit. The estimated rating volatilities are consistent with the fact that the non investment grade firms are likely to move in a larger set of ratings. Indeed, σ7 = 1.72 and σ6 = 1.17, whereas σ1 = 0.17 and σ2 = 0.38. One should also notice that the β coefficients for the non investment grade ratings 6 and 7 are significantly higher than for the investment grade firms. This confirms that the weakest firms are more subject to the uncertainty of the economy.

Estimation with unobservable latent factors

The use of latent factors to explain time variation of the credit ratings requires to compute a relatively high number of empirical transition intensities λij,t , i ∈ {1, • • • , 7} , j ∈ {1, • • • , 8}, i = j. Indeed, λij,t are used as inputs for both models, they are transformed into ψiJ,t for the structural model (see next section) and kept as such for the intensity model. In order to increase the number of observations, we have considered an overlapping year to year sub-period with one month slope. With 8 years of monthly frequency data, this allows to have 86 observations instead of 8 observations if the sub-periods were not overlapping10 . For instance, λij,t1 is calculated on the sub-period [01/2006 -01/2007], λij,t2 is calculated on the sub-period [02/2006 -02/2007], etc. We recall that λij,t are computed using the proposition 3.3.

Determining the number of unobservable factors

Since the covariates are unobserved we proceed by using the method proposed in [START_REF] Gagliardini | Stochastic migration models[END_REF] to find the number of covariates in both structural and intensity models. The idea is to build time series of λij,t and ψiJ,t for t = t 1 , • • • , t 86 , to consider them as the realisations of random variables and then carry on a PCA to find the number of common factors.

For the intensity model, we count 49 variables λ12,t , λ13,t , • • • , λ78,t , each one going from t 1 to t 86 which makes 49 time series. The PCA is then applied on this set to compute the common factors. Concerning the structural model, the λij,t are transformed into empirical transition probabilities pij,t using (5), then the series ψiJ,t such as ψiJ,t = Φ -1 p * For the λij,t time series, the first two factors explain 98% of the total variance. The variance of ψiJ,t time series are explained at 72% by the two first factors. We have applied the Kalman filter for one common factor and two common factors, we found that the filtered process was the same in both cases. We then consider that the first factor represent the business cycle. Once the number of covariates known, we have done the parameters estimation using the two steps procedure, filtering of the process Xt and maximisation of the log-likelihood (see sections 3.3 and 4.2). We present and discuss the results for the one factor models below.

Results and discussion

We report in the Figures (1) and (2) the filtered processes 2), ( 3) and ( 4)). Parameters estimation: the estimated parameters for the one factor models are reported in the following tables. See Tables ( 12) and ( 13) for the one factor intensity parameters and The parameters obtained for both models are consistent with those obtained in the case of observable factors. However, the values β1 = 0.43 and σ1 = 1.71 are not in line with the previous results. This is more a numerical "issue" than a interpretable result.

E(X t | Y t ) and E(X t | Ψ t ) respectively 11 , for t = t 1 , • • • , t 86 . Recall that E(X t | Y t ) and E(X t | Ψ t ) correspond

Comparison between models

As stated previously, the aim of this paper is to assess and compare the migration models on their ability to link the transition probabilities to dynamic risk factors. In order to achieve this objective, we proceed in two steps. First, we compare the one year default probabilities implied by each model during the study period with the empirical default probabilities. Second, we compare transition matrices of each model with respect to a metric called SVD mobility index and introduced in [START_REF] Jafry | Measurement, estimation and comparison of credit migration matrices[END_REF]. The mobility index is a function that measures the ability of a transition matrix to generate migration events (the mobility index is detailed in the section 5.5.2). The model parameters used for the computation of one-year default probabilities and mobility indices are those obtained from the estimation procedures of the two previous sections (observable and latent covariates). The comparison is made with respect to a credit migration model with no covariate (called the empirical model) where, on each sub-period, the transition matrix is estimated using the procedure described in proposition 3.3 (see also section 5.2). Recall that we have 8 non overlapping sub-periods for the observable factors and 86 overlapping sub-periods for the unobservable factors.

Comparison using the implied one year default probabilities

We compute the one year matrices implied on each sub-period by both intensity and structural models. We present graphically the evolution of the implied default probabilities of the ratings 4, 5, 6 and 7 (the default probabilities of the ratings 1, 2, and 3 are almost null). Concerning the one year PD for the ratings 6 and 7 (see Figure ( 6)), the structural model shows a quasi non sensitivity to the business cycle. Indeed, the curves are high and almost flat comparing to the empirical one year PD. The intensity model shows a better adjustment in the sense that it is more reactive than the structural model to the business cycle but it is still unable to reproduce the empirical PD during the trough period.

Comparison using the SVD mobility index

In literature, comparing two transition matrices is usually done using three methods, the euclidean distance (see [START_REF] Israel | Finding generators for Markov chains via empirical transition matrices, with applications to credit ratings[END_REF]), statistic testing, t-test or χ 2 (see [START_REF] Nickell | Stability of rating transitions[END_REF], [START_REF] Foulcher | Term structure of defaults and ratings[END_REF]) or with a mobility index (see [START_REF] Geweke | Mobility indices in continuous time Markov chains[END_REF], [START_REF] Jafry | Measurement, estimation and comparison of credit migration matrices[END_REF]). The mobility index was introduced to compare two transition matrices on their ability to generate migration events. It is defined as a function M (Π) of R d×d -→ R where by convention M (I) = 0, with I the identity matrix. A good mobility index must be assessed on three criteria, persistence, convergence and temporal aggregation (see the well documented papers of [START_REF] Shorrocks | The measurement of mobility[END_REF], [START_REF] Geweke | Mobility indices in continuous time Markov chains[END_REF] and [START_REF] Jafry | Measurement, estimation and comparison of credit migration matrices[END_REF] to have a complete overview of this topic). The interest of the mobility index, is that instead of comparing the transition probabilities individually it compares the matrices according to their overall dynamic. We present hereafter the SVD mobility index introduced in Jafry and Schuermann ( 2004), we will use it as a metric to compare the matrices implied by the intensity and the structural models.

For a d-dimensional square matrix Π, the SVD mobility index is given by,

M SV D (Π) = d i=1 λ i Π Π d , (40)
where Π = Π -I and Π is its transpose. λ i Π Π is the i th eigenvalue of Π Π, sorted in decreasing

order, i.e. λ 1 ( Π Π) > • • • > λ d ( Π Π).
We have applied the SVD mobility index on each matrix calculated according to the intensity and structural models. Figures (7) As a conclusion for the section 5.5 we can state that the intensity model with observable factors is the one that best fits the empirical probabilities. This is mainly due to following points:

• for the intensity model with observable factors, the parameters λ ij,0 and β ij are estimated through the maximization of the marginal likelihood L ij . This means that among the set of explaining factors, it is most likely to have at least one factor that explains correctly the intensity migration λ ij . The estimation of the structural model parameters is done through the maximization of the global likelihood function, this is less tractable than the intensity model in the sens that it is not possible to explain separately the transition probabilities.

• for the intensity model with observable factors, the maximum likelihood problem is solved using a semi-analytical solution which is not the case of the other models. Indeed the baseline intensity estimator λij,0 is given analytically when the regression coefficients β are obtained through the resolution of the gradient of the log-likelihood function which represents a simple non-linear equations system. The parameters of the structural model with observable factors are also estimated through the resolution of the gradient of the log-likelihood function but the non-linear equations system is far more difficult to resolve.

• for the statistical estimation using unobservable factors, the solution of [START_REF] Gagliardini | Stochastic migration models[END_REF] to skip the issue of high dimension integrals requires to build a linear Gaussian model to apply the Kalman filter. The assumption of an infinite number of obligors n and the asymptotic normality of the maximum likelihood estimator at each time t and for every migration (i,j) is not satisfied. Indeed, even with a very high value of n, some migrations are unlikely to happen frequently like the long migrations (three grades or more). The intensity model with observable factors use directly the observations without making some approximations based on very restrictive assumptions.

• the latent factors models have more parameters to estimate (the filtered process Xt and the autoregression matrix A). This constraint added to the limited quantity of data plays an important role in the final results.

Conclusion

The aim of this paper is to assess and compare two alternative stochastic migration models on their ability to link the transition probabilities to either observable or unobservable dynamic risk factors.

In this respect, we use the same data set (S&P credit ratings on the period [01/2006 -01/2014]) to estimate model parameters in the multi-state latent factor intensity model and in the ordered Probit model. The estimation procedure is detailed in the two cases where the underlying dynamic factors are assumed to be either observable or unobervable. When the underlying factors are unobservable, we adapt a method given in [START_REF] Gagliardini | Stochastic migration models[END_REF] to represent the considered factor migration model as a linear Gaussian model. In that case, we identify the business cycle as the factor process filtered by a standard Kalman filter. The estimation methods are compared on their ability to fit the one-year empirical default probabilities and their ability to reproduce the empirical SVD mobility index. We conclude that the intensity model with observable factors was the one which best fits the S&P rating history. In line with [START_REF] Kavvathas | Estimating credit rating transition probabilities for corporate bonds[END_REF], this study shows that short migrations of investment grade firms are significantly correlated to the business cycle whereas, because of lack of observations, it is not possible to state any relation between long migrations (more than two grades) and the business cycle. Concerning non investment grade firms, downgrade migrations are negatively related to business cycle whatever the amplitude of the migration.

Stepping back recursively on H l u leads to get, B Proof of Proposition 3.4

Proof. Contrary to the homogeneous case, the maximum likelihood estimator is not only related to λ ij,0 but to λ ij,0 and the vector β ij , we denote our vector parameters θ = (λ ij,0 , β 1,ij , β 2,ij , • • • , β K,ij ) and look to estimate θ such as,

∂ log (L ij (θ | G t )) ∂θ = 0 ⇒ ∂ log (L ij (θ | G t )) ∂λ ij,0 , ∂ log (L ij (θ | G t )) ∂β ij = 0. ∂ log (L ij (θ | G t )) ∂λ ij,0 = n l=1 N l k=1 Y l ij (t l k ) ∂ log λ ij (X t l k ) ∂λ ij,0 - n l=1 N l k=1 S l i (t l k ) t l k t l k-1 ∂λ ij (X u ) ∂λ ij,0 du.
We recall that, Getting vector βij goes through solving equation system ∂ log (L ij (θ | G t )) ∂β ij = 0, every partial derivative will result in a relation between vector β ij and past realised covariates values (i.e. X s,t ) as well as migration observations. This non linear equation system has no analytical solution, it can be solved by numerical optimisation (multi-dimensional Newton-Raphson for instance). The solution provides the values of βij .

λ ij (X t ) = λ ij,0 e βij ,Xt
∂ log (L ij (θ | G t )) ∂β ij = 0 ⇔                  ∂ log (L ij (θ | G t )) ∂β 1,ij = 0 . . . ∂ log (L ij (θ | G t )) ∂β 2,ij = 0 . . . . . . ∂ log (L ij (θ | G t )) ∂β K,ij = 0                  , (45)
Once we have βij we go back to λ ij,0 to get its estimator λij,0 .

C Maximisation of equation (31) using gradient

One can notice that the gradient of (31) with respect to α i , β i and σ i can easily be derived. This means that it is possible to transform the maximization problem to a system of non-linear equations and then to resolve this system.

∂ log (L(θ | G t N ∨ H t0 )) ∂θ = 0 ⇒ ∂ log (L(θ | G t N ∨ H t0 )) ∂α i , ∂ log (L(θ | G t N ∨ H t0 )) ∂β i , ∂ log (L(θ | G t N ∨ H t0 )) ∂σ i = 0. For i ∈ {1, • • • , d -1}, with a fixed values of threshold C j (-∞ = C d+1 < C d < • • • < C 2 < C 1 =
+∞), every vector θi = αi , βi , σi is obtained through the resolution of the following non-linear equations system,

                                                 n l=1 N k=1 j =i Y l ij,t k p ij (X t k ) φ C j -α i -β i , X t k σ i -φ C j+1 -α i -β i , X t k σ i = 0 n l=1 N k=1 i =j Y l ij,t k p ij (X t k ) (C j -α i -β i , X t k ) φ C j -α i -β i , X t k σ i - (C j+1 -α i -β i , X t k ) φ C j+1 -α i -β i , X t k σ i = 0 n l=1 N k=1 i =j Y l ij,t k X t k ,s p ij (X t k ) φ C j -α i -β i , X t k σ i -φ C j+1 -α i -β i , X t k σ i = 0 for s = 1, • • • , K (46)
Where, φ(•) is the density function of the standard normal distribution. By convention, we have,

                         φ C 1 -α i -β i , X t k σ i = 0 φ C d+1 -α i -β i , X t k σ i = 0 (C 1 -α i -β i , X t k ) φ C 1 -α i -β i , X t k σ i = 0 (C d+1 -α i -β i , X t k ) φ C d+1 -α i -β i , X t k σ i = 0.
(47)

  to the filtered factor process associated with max-likelihood models (see (3.4) and (4.2)). The figures below clearly show periods of peak, normal times and trough which are properties of the business cycle. Our period of study encompasses what it is known as the Liquidity Crisis and Sovereign Debt Crisis, with some famous events like the Lehman Brothers bankruptcy on September 2008 and the falls of market indices like CAC40, DAX and FTSE100 during the summer 2011. These crisis are perfectly captured by the structural model while the intensity model fails to capture the sovereign debt crisis. The study of the default probability of non investment grade ratings shows two peaks during the periods [01/2009 -01/2010] and [01/2012 -01/2013] which correspond to the trough periods of the business cycle (see Figures (1), (

Figure 1 :

 1 Figure 1: Filtered common factor in the intensity model

  56 σ3 = 0.69 β3 = -0.06 Ĉ5 = 1.84 α4 = 1.83 σ4 = 0.69 β4 = -0.05 Ĉ6 = 0.96 α5 = 1.23 σ5 = 0.75 β5 = 0.03 Ĉ7 = 0.001 α6 = 0.12 σ6 = 0.87 β6 = 0.03 Ĉ8 = -0.43 α7 = -0.56 σ7 = 1.09 β7 = -0.03Table 14: Parameters for 1 factor structural model

  Observable factors: Figures (3) and (4) represent the estimated one year default probabilities in the empirical model (solid line), in the intensity model (dashed line) and in the structural model (dotted line). One can see on the Figure (3) that for ratings 4 and 5 the intensity model adjusts perfectly the empirical default probabilities, the curves are almost identical. The structural model over-estimates the default probabilities for the rating 5, when the curve of the default probabilities for the rating 4 is almost flat.

Figure 3 :

 3 Figure 3: Ratings 4 and 5: implied one year default probabilities. Model estimation performed in the case of observable factors. Empirical model in solid line, intensity model in dashed line, structural model in dotted line.

Figure 4 :

 4 Figure 4: Ratings 6 and 7: implied one year default probabilities. Model estimation performed in the case of observable factors. Empirical model in solid line, intensity model in dashed line, structural model in dotted line.

  Unobservable factors: Figure (5) shows that both models are not able to reproduce the one year empirical PD during the trough period, see the period [01/2009 -07/2010]. Concerning the rest of the study period, the structural model over-estimates the PD for the rating 5 where it shows an almost flat curve for the PD rating 4. The PD curves of the intensity model (for ratings 4 and 5) are most of the time under the empirical PD.

Figure 5 :

 5 Figure 5: Ratings 4 and 5: implied one year default probabilities. Model estimation performed in the case of unobservable factors. Empirical model in solid line, intensity model in dashed line, structural model in dotted line.

Figure 6 :

 6 Figure 6: Ratings 6 and 7: implied one year default probabilities. Model estimation performed in the case of unobservable factors. Empirical model in solid line, intensity model in dashed line, structural model in dotted line.

  and (8) show the evolution of the SVD mobility index of the empirical model (black curve), of the on intensity model (red curve) and of the structural model (blue curve).

  Figure (7) corresponds to estimations with observable factors whereas Figure (8) corresponds to estimations with unobservable factors. The two figures confirm the statements made previously : the intensity model with observable factors fits very well the empirical transition matrices whereas the other models, i.e., the structural with observable factors, the intensity with latent covariates and the structural with latent covariates, over-estimates, under-estimates and over-estimates respectively the empirical transition matrices.

Figure 7 :

 7 Figure 7: Implied SVD mobility index implied when model estimation is performed in the case of observable factors. Empirical model in black, intensity model in red, structural model in blue.

Figure 8 :

 8 Figure 8: Implied SVD mobility index implied when model estimation is performed in the case of unobservable factors. Empirical model in black, intensity model in red, structural model in blue.

  for k = N l .

  and ∂λ ij (X t ) ∂λ ij,0 = e βij ,Xt . ∂ log (L ij (θ | G t )) ∂λ ij

Table 5 :

 5 One year Through The Cycle transition matrix

	)).
		1	2	3	4	5	6	7	8
	1 92.69 6.47	0.48	0.03	0.31	0.02	0	0
	2 0.64 89.93 8.96	0.46	0.01	0	0	0
	3 0.02	1.47 92.37 5.91	0.20	0.03	0	0
	4	0	0.08	3.12 92.13 4.16	0.41	0.06	0.04
	5 0.06	0.09	0.25	6.89 83.99 8.16	0.41	0.15
	6	0	0	0.13	0.52	7.96 82.77 6.10	2.52
	7	0	0	0.03	0.61	1.24 22.92 43.27 31.93
	8	0	0	0	0	0	0	0	100

Table 6 :

 6 , which gives, As the number of covariates is important (25 series), we have performed a Principal Component Analysis (PCA) in order to reduce the number of explicative factors.Table (6) summarizes the 8 first eigenvalues with cumulated variance explained in percentage. Observable factors eigenvalues

	Real GDP Growth 01/2014 =	Real GDP 31/12/2013 -Real GDP 31/12/2012 Real GDP 31/12/2012
	=	15916 b$ -15433 b$ 15433 b$	= 3.13%
	Eigenvalues	Percentage	Cumulated
		of variance	percentage
	8.79	43.99	43.99
	3.39	16.98	60.98
	2.43	12.19	73.18
	1.61	8.09	81.27
	0.95	4.77	86.04
	0.72	3.64	89.69
	0.56	2.83	92.53
	0.36	1.83	94.36

Table 7 :

 7 Baseline intensity generator for observable factors intensity model

		1	2	3	4	5	6	7	8
	1	0	-0.06	0	0	0	0	0	0
	2 0.34	0	-0.11 1.58	0	0	0	0
	3 -4.93 0.22	0	-0.03 0.28	2.77	0	0
	4	0	12.38 0.02	0	-0.06 0.07	0.80 -2.47
	5 -2.42 -2.42 -0.17 0.04	0	-0.04 -2.92 -0.17
	6	0	0	-7.67 1.63	0.05	0	-0.04 -0.08
	7	0	0	0	3.52	0	0.09	0	-0.11
	8	0	0	0	0	0	0	0	0

Table 8 :

 8 Beta for the first factor intensity model

Table 9 :

 9 Parameters for observable factors structural model

	β1	β2	β3	β3	β5
	β11 =0.01 β21 = 0.03	β31 =0.05	β41 =-0.06	β51 =0.27
	β12 =0.03	β22 =0.02	β32 =-0.02	β42 =0.09	β52 =-0.25
	β13 =0.06 β23 = 0.06 β33 =-0.06	β43 =0.08	β53 =-0.17
	β14 =0.06	β24 =0.04	β34 =0.00	β44 =-0.01	β54 =0.00
	β15 =0.04	β25 =0.05	β35 =-0.06 β45 =-0.03	β55 =0.01
	β16 =0.07	β26 =0.13	β36 =0.14	β46 =-0.19	β56 =0.40
	β17 =0.14	β27 =0.23	β37 =0.22	β47 =0.07	β57 =0.38

Table 10 :

 10 Beta for observable factors structural model

  iJ,t for t = t 1 , • • • , t 86 are computed (see (35) in 4.2 to get the definition of p * iJ,t ). We get also 49 times series on which we apply the PCA. The corresponding eigenvalues are given in decreasing order in the following table.

	intensity	Cumulated	structural	Cumulated
	eigenvalues	percentage	eigenvalues	percentage
	1.83	79.50	1853.46	48.89
	0.43	98.28	735.23	72.29
	0.02	99.27	515.81	81.52
	0.007	99.58	311.41	89.04
	0.003	99.71	306.01	93.81
	0.002	99.82	247.94	96.50
		Table 11: Eigenvalues

Table (

 (

									14) for
	the one factor structural parameters.
		1	2	3	4	5	6	7	8
	1 -0.03 0.03	0	0	0	0	0	0
	2	0	-0.08 0.08	0	0	0	0	0
	3	0	0	-0.07 0.07	0	0	0	0
	4	0	0	0.03 -0.07 0.04	0	0	0
	5	0	0		0.07 -0.17	0.1	0	0
	6	0	0	0	0	0.07 -0.15 0.08	0.
	7	0	0	0	0	0	0.27 -0.67 0.42
	8	0	0	0	0	0	0	0	0
		Table 12: Baseline intensity for one factor intensity model
		1	2	3	4	5	6	7	8
	1	0	-0.25	0	0	0	0	0	0
	2 0.43	0	-0.12 0.20	0	0	0	0
	3 0.02	0.19	0	-0.05 0.16	0.05	0	0
	4	0	0.07	0.03	0	-0.06 0.06	0.05 -0.23
	5 -0.16 -0.27 -0.23 0.03	0	-0.07 -0.46 -0.18
	6	0	0	-0.25 0.12	0.12	0	-0.15 -0.28
	7	0	0	0	-0.19	0	-0.03	0	-0.11
	8	0	0	0	0	0	0	0	0

Table 13 :

 13 Beta for one factor intensity model

  By replacing with λ ij,0 's expression we get the equation system, X s,u e βij ,Xu du

	n l=1	N l k=1 Y l ij (t l k ) X s,t l k l=1 n N l k=1 Y l ij (t l k )	=	n l=1	N l k=1 S l i (t l k ) l=1 n N l k=1 S l i (t l k ) t l k t l k-1 k k-1 t l t l	e βij ,Xu du	.
	for s = 1, • • • , K
		n l=1	N l k=1	Y l ij (t l k ) X s,t l k -λ ij,0	n l=1	N l k=1	S l i (t l k )	t l k t l k-1

∂ log (L ij (θ | G t )) ∂β s,ij = X s,u e βij ,Xu du = 0

See Gagliardini and Gourieroux (2005) for more details.

Through The Cycle transition matrix corresponds to the long term time-invariant transition matrix of the homogeneous Markov chain. The transition probabilities are almost unaffected by the economical conditions.

[START_REF] Koopman | The multi-state latent factor intensity model for credit rating transitions[END_REF] propose to use a suitable important sampling technique to improve the efficiency of the Monte Carlo estimator.

Following Hougaard (2000), the estimate σij,t of σ ij,t is given by σij,t = N ij,t (T i,t ) 2 , where N ij,t is the total number of transitions from i to j (numerator of (19) for t ∈ [t k-1 , t k )), and T i,t is the total time spent in state i (denominator of (19) for t ∈ [t k-1 , t k )).

Ω 2 iJ,t is approximated by Ω2 iJ,t = p * iJ,t (1 -p * iJ,t ). Indeed, if we consider the binomial variable : starting from i, either go to J or below, it is clear that the standard error for p * iJ,t can be calculated as a binomial standard error. see[START_REF] Nickell | Stability of rating transitions[END_REF].

The macroeconomic indicators time series are in free access on the website of the FEDERAL RESERVE BANK OF ST. LOUIS, where the financial, commodity and credit indicators time series can be retrieved from Bloomberg.

We need a high number of observations in order to perform the Kalman filter, we have performed the Kalman filter on 8 observations but we didn't get a satisfying result.

Note that the process E(Xt | Ψt) has been centered in the Figure(2) in order to be easily interpreted.

Appendices

A Proof of Proposition 3.1

We recall the model specification, for each obligor l = 1, • • • , n;

• we denote R l u the random variable representing the rating of obligor l at time u, R l u ∈ {1, • • • , d}, where d is the default state and 1 the best credit quality state;

• we fix a finite horizon date t > 0. The underlying probability space (Ω, F, P) is endowed with some reference filtration F = (F u ) 0≤u≤t and is sufficiently rich to support the stochastic processes R l u , 0 ≤ u < t and {X u , 0 ≤ u ≤ t};

• let H l u(0≤u<t) be the filtration associated to the process R l u , 0 ≤ u < t for every value of l = 1, • • • , n. We consider a global filtration H u(0≤u<t) which is the union of all the filtrations H l u , that is,

• let G u(0≤u≤t) be the filtration associated to the process {X u , 0 ≤ u ≤ t}. We assume that conditional on G t , the processes R l u , 0 ≤ u < t are independent across obligors;

• we assume that intensities are adapted to the filtration G u(0≤u≤t) and independent from the filtration H u(0≤u<t) .

• we denote N l the number of ratings visited during the period [0, t) by the obligor l, with the condition

corresponds to the visiting of the state r l k , we then have the visited states r l 1 , • • • , r l N l . Note that each consecutive states are obligatory different, thus r l k-1 = r l k ;

• there are N l -1 migration events for each obligor, each migration event takes place at the time

• we denote r l u , 0 ≤ u < t the observed path of ratings of obligor l during the period [0, t); Proof. We know that,

The equation ( 42) describes the likelihood of the entire path of the obligor l. The obligor starts in the rating r l k at time t l k-1 , it spends the period t l k-1 , t l k in this state until it jumps to the rating r l k+1 at time t k . The last period t l N l -1 , t l N l is particular because we only observe the rating R l u until time t l- N l (which corresponds to time t -). We can only conclude on the stay of the obligor l in the state r l N l .

The Tower Property is used with respect to G t ∨ H l t N l -1 in order to exhibit the probability of staying in the rating r l N l during the period t l N l -1 , t l N l . The rest of the path is also obtained with a recursive use of the Tower Property by stepping back each time with one period on H l u . We show the case of G t ∨ H l t N l -2 and deduce the rest to get (13).


Alban Vergnaud
email: alban.vergnaud@univ-angers.fr

Thanh Phong Tran

Laetitia Perez

Philippe Lucidarme

Laurent Autrique

Deployment strategies of mobile sensors for monitoring of mobile sources: method and prototype

Keywords: inverse problem, parametric identification, thermal procedure, Conjugate Gradient Method, Experimentation

In this communication, the characterization of mobile sources is studied. A quasi online identification method based on the iterative regularization of ill-posed inverse problem is proposed, it allows monitoring of evolutionary phenomena. Several strategies of choice of sensors (Adaptive selection of relevant sensors within the fixed network or management of a pack of mobile sensors) are discussed. To illustrate the approaches, an experimental prototype has been designed in the context of thermal engineering. The experiment is described: design, characterization, exchanges protocols of signals acquisition, movement of mobile elements (sources and sensors).

I. INTRODUCTION

The distributed parameters identification (continuously depending on space and time) in systems described by nonlinear partial differential equations (PDE) is, in general, difficult. A typical approach is to minimise an output error problem. For example, to minimize a quadratic criterion describing the difference between outputs predicted by the mathematical model and observed outputs collected by using sensors. This approach can be defined as an inverse problem [START_REF] Isakov | Inverse problems for Partial Differential Equations[END_REF].

In more specific context of thermal engineering, inverse problems of heat conduction studied to identify properties, heat flows or initial conditions, are usually ill-posed [START_REF] Alifanov | Inverse heat transfer problems[END_REF]. In fact, the stability condition (in Hadamard's sense [START_REF] Isakov | Inverse problems for Partial Differential Equations[END_REF]) is generally not satisfied because of weak disturbances on measurements generate large variations on the identified parameters. This can be shown in a matrix environment where, in the numerical solution of partial differential equations system, assembled matrices obtained by finite element method [START_REF] Pepper | The finite element methodbasic concepts and applications[END_REF] are ill-conditioned.

Matrices inversion is not relevant in the presence of uncertainties on measurements [START_REF] Khachfe | Résolution numérique de problèmes inverses 2D non linéaires de conduction de la chaleur par la méthode des éléments finis et l'algorithme du gradient conjugué-Validation expérimentaleThesis of the University of Nantes[END_REF][START_REF] Beddiaf | Identification paramétrique de systèmes d'EDP paraboliques non linéaires en géométrie 3D par une méthode de régularisation itérative[END_REF]. In this context, Tikhonov has developed methods allowing the construction of stable solutions [START_REF] Isakov | Inverse problems for Partial Differential Equations[END_REF][START_REF] Alifanov | Inverse heat transfer problems[END_REF][START_REF] Morozov | Methods for solving incorrectly posed problems[END_REF]. These methods are generally called as "regularization method ". The main principle is to formulate a new problem considering an additional parameter (the regularization parameter) such that the new problem is stable. The construction of these new problems is not trivial.

Among various existing approaches, the Conjugate gradient method [START_REF] Minoux | Programmation mathématiquethéorie et algorithmes[END_REF] allowing the iterative minimization of the output error has a regularisation property in context of the inverse problems of heat conduction.

Indeed, this method leads to solve iteratively three problems well posed [START_REF] Jarny | A general optimization method using adjoint equation for solving multidimensional inverse heat conduction[END_REF]. In [START_REF] Alifanov | Inverse heat transfer problems[END_REF], Alifanov says: "such a method of damping the instability when specifying an approximate solution for an ill-posed problem is based on viscous properties of numerical algorithms of optimization". This stabilizing effect is illustrated in [START_REF] Prud'homme | On the iterative regularization of inverse heat conduction problems by conjugate gradient method[END_REF] for an academic situation corresponding to a 1D geometry for which analytical solutions are available. It is highlighted that the main structure of the unknown parameter is identified from early iterations.

The Conjugate gradient algorithm behaves as a filter able to reject disturbances on measurements during the iterative process of minimization. The number of iteration is often considered to be the parameter of regularization of these methods.

This inverse problem resolution is traditionally implemented offline: identification begins when all data are collected [START_REF] Autrique | On the implementation of a finite element method for parameter identification[END_REF][START_REF] Perez | Implementation of a conjugate in a three-dimensional domain[END_REF][12][START_REF] Beddiaf | Parametric identification of a heating mobile source in a three-dimensional geometry[END_REF]. However recent adaptations have shown its potential for a quasi-online use in a numerical context [START_REF] Vergnaud | On-line monitoring of surfacic mobile heating sources[END_REF]. To illustrate the interest of this new approach in a real situation, experimentation was designed for the identification of heat mobile sources (See Fig. 1). This main objective of the prototype is to estimate the trajectory and the heating of one or more surfacic heating sources from temperature measuring.

II. EXPERIMENTAL DEVICE

In this part, the studied thermal situation is described. It is intended to illustrate the quasi online resolution of an inverse problem of heat conduction using the iterative adjustment of the conjugate gradient method. It also allows testing strategies of sensors choice (positioning, deployment).

A. Field of study

In this study, the hypothesis of a geometry 2D (infinitely fine plate) was considered in order to reduce the computing times. This assumption is valid only if we can neglect the heat transfer in the thickness of the plate. To do this, it is necessary to consider a fine plate of a metal having a high thermal conductivity. For financial reasons, aluminum was chosen (gold, silver and copper not being possible). This square plate 3m of the each side posed horizontally on a support providing insulation. This last is composed of Rockwool thermal insulation to high temperatures.

In order to verify that the heat transfers are surfacing on the aluminum plate, two mathematical models were compared.

The first system [START_REF] Isakov | Inverse problems for Partial Differential Equations[END_REF]

                                      (1)
where c  is the volumetric heat capacity in The reference of this material is Rockwool-Rockmur-Kraft.

-3 -1 J.m .K ,  is the thermal conductivity in -1 -1 W.m .K , 0  in K
The mathematical model describes the heat transfer in 3D geometry is the following: 2) are solved by using the finite element method [START_REF] Pepper | The finite element methodbasic concepts and applications[END_REF] implemented by the Comsol code interfaced with Matlab. On figure 1, the evolutions of temperature in various points of the superior face of the plate are submitted to a plate thickness e  2mm.

              0 0 , , ; 0 , , , ;0 , , ; 0 , , ; 1 ,
                                                                  (2)

Fig. 1. Comparison of the 2D and 3D models

Figure 1 shows that the 2D model is very satisfactory for points near the heating source. For points far away from the source, weak errors of model will be present. The following table indicates in various points the mean absolute error between two models, as well as the absolute average deviation of temperature emod between 2 points located on both sides of the plate (information ecart obtained with the model ( 2)). component is only related to controlling the power supply.

Indeed, the delay inherent in the heating of the filament is considered negligible compared with the duration of the experimentation (approximately 30 minutes) and with dynamic delays for the rise in temperature. Thus, it is justified to write for each source

      , ; , heat x y t t f x y   where   t 
is function of the power provided by the power supply and   , f x y is spatial distribution of the heat flux density on the plate. Identification campaigns must be realized to know these parameters. The function   , f x y depends on the geometry of the bulb and the distance to the plate. A reflector placed behind the bulb allows increasing the flux received by the plate. This last depends on the distance between the heat source and the plate. If a constant flux on a delimited surface is desired then a kaleidoscope (device for homogenizing the flux) can be used; see an application in [START_REF] Vergnaud | DARC-EDP: Conception du prototype[END_REF]. This last will require that the distance between the plate and the heat source is fixed (equal to the length of the kaleidoscope). In a general case, it is necessary to identify the spatial distribution   The robots are connected to the management computer via a wireless router (wifi router) through which they can exchange information in the form of frames using the TCP/IP protocol (client/server).

In order to obtain the precise position of each robot, a visible camera (camera AVT StingrayFireWire, F-046C, field of view 61.9° -76.7°) is wired with the management computer. This camera has 780x580 pixels and a frequency acquisition at 55Hz. By using the software SSL-vision [START_REF] Martinez-Gomez | Real Time Vision System for a Small Size League Team[END_REF][START_REF] Zickler | SSL-Vision: The Shared Vision System for the RoboCup Small Size League[END_REF], the displacement of the robots on the plate is registered and the coordinates as well as the orientations of the robots are sent to the management computer to recalibrate. In fact, the robots have various sensors (10 infrared sensors, 5 ultrasonic telemeters) and wheels with incremental encoders which allow knowing their relative positions but these would suffer an important drift. By using the camera, the accuracy of the position of the robots is about the centimeter.

E. In practice…

The temperatures as well as the precise positions of the robots are measured every second. The robots embedding the heating sources can move continuously. The robots embedding the pyrometers move according to the provided indications every second by the management computer (possibly connected to a second computer dedicated to the resolution of the problem of identification).

F. Vision system prototype

The heating sources are embedded on mobile robots which deploy on a plate. The power supply of heat sources is controlled so as to describe various profiles of evolution. On

H. Localization based on the vision

The robots Khepera III have wheels encoders able to calculate the relative displacement carried out at the time of the starting of the motors (odometric measurement). However, frictions between the wheels of the robot and the plate generate bias (in particular at the time of tight turns). In order to correct errors of trajectory, the robot location checking is implemented via image processing software (software SSLvision [START_REF] Martinez-Gomez | Real Time Vision System for a Small Size League Team[END_REF][START_REF] Zickler | SSL-Vision: The Shared Vision System for the RoboCup Small Size League[END_REF]). This software uses the images provided by the camera and detects markers previously recorded (figure 5).

These are used as a unique identifier; by placing a marker on each robot, it is possible to know its position  

III. IDENTIFICATION

A. Direct problem

When all the parameters of the model are known

      2 * 1 1 Ârg min , , 2
N heat i heat i i C t t dt             T (3)
Where N is number of sensors i C and T the time interval during which measures are taken into account.

In following paragraphs, sensor choice (fixed or mobile) strategies are considered and the definition of the intervals of observations T is discussed.

C. The conjugate gradient method (Offline method)

To minimize the cost function defined in (3), the iterative Conjugate Gradient Method (CGM) is implemented [START_REF] Minoux | Programmation mathématiquethéorie et algorithmes[END_REF]. An application to thermal engineering is presented in [START_REF] Jarny | A general optimization method using adjoint equation for solving multidimensional inverse heat conduction[END_REF] and its regularization properties are illustrated in [START_REF] Prud'homme | On the iterative regularization of inverse heat conduction problems by conjugate gradient method[END_REF]. The selected algorithm is solving iteratively three well posed problems on a time interval. At the end of the resolution of these three problems, a new value of



1 k heat   is calculated.           2 ,; , , ; 0 ,
; , ; 0

h x y t c E te x y x y x y t x y t n                                      T T (4)
For the Adjoint problem (4),

 

,;

E x y t depends on the difference between the estimated temperature and measurements at sensor position. Note that this problem is retrograde in time (

f t   and f t is equal to the end of the time interval).                 ,;
, , ;0 0 ,; , ; 0

k x y t c t x y x y x y t x y t n                                      T T (5)
The variation of the heat flux density used in sensitivity problem is: 4) and ( 5), the equations to obtain the gradient (with  

    2 k heat k h e      Problems (
,;

x y t



) and the formulation of the descent depth are explained in [START_REF] Beddiaf | Identification paramétrique de systèmes d'EDP paraboliques non linéaires en géométrie 3D par une méthode de régularisation itérative[END_REF][12][START_REF] Beddiaf | Parametric identification of a heating mobile source in a three-dimensional geometry[END_REF][START_REF] Vergnaud | On-line monitoring of surfacic mobile heating sources[END_REF]. The iterative process continues until the cost function (3) is under a permissible limit [START_REF] Morozov | Methods for solving incorrectly posed problems[END_REF] [2]:

2 1 2 stop Jn   (6)
Where n is number of measurements i C considering on time interval T and  the standard deviation of the noise measurements (assumed as a Gaussian profile with mean value equal to zero). This method is offline when T  T , the identification process begins only when all measurement are collected. It is thus necessary to wait the end of the experiment to identify the trajectory and power of heating sources.

D. Adaptation for a quasi-online method

For identification during the process, when measurements are available, the previous method should be adapted. The choice of these sliding time intervals (start, length, recovery rates), is explained in [START_REF] Vergnaud | On-line monitoring of surfacic mobile heating sources[END_REF] where several strategies are presented. Note in [START_REF] Vergnaud | On-line monitoring of surfacic mobile heating sources[END_REF] that the experimental design and the choices of sensor procedures are not presented.

Techniques for online identification are based on the behavior of the minimization algorithm to handle compromises between quality of the estimate and delay between estimation and measurements.

The algorithm is started only if new observations are not in adequacy with predicted temperature evolution. Lengths of time intervals must be choice considering the dynamics of the experimentation (priori experimental knowledge can be taken into account). Knowledge on the physical system studied is an undeniable asset to 'adjust' the online procedure. This note also applies to the choice of sensors presented in the next paragraph.

In the considered situation, it is shown that offline identification leads to results 70 minutes after the end of the experiment of 10 minutes, the online method provides satisfactory results only 1 minute after the end. This approach also allows estimation of the unknown parameter even before the end of the process (the average difference between the measurements and the results of online method is about 30 seconds).

E. Choice of relevant sensors.

In order to identify unknown trajectories of heating sources, a network of fixed sensors can be used. Indeed, temperature sensors (thermocouples) are not very expensive. However, this type of sensors introduces a non-negligible noise on measurement. On a plate (9m²), it's clear that several sensors are not sensitive to sources displacement. Also the problem concerns the selection of a few "relevant" sensors within a large number of "blind" sensors.

In [19], methods based on the optimal design of experiments are presented. In general, these approaches require a priori knowledge of a nominal value for the unknown parameter.

Below, the proposed method is based on analysis of the iteratively solved sensitivities functions [START_REF] Beddiaf | Identification paramétrique de systèmes d'EDP paraboliques non linéaires en géométrie 3D par une méthode de régularisation itérative[END_REF].

On the time interval T the sensors with sufficient sensitivity are selected (an example inspired by this approach is presented in [20]). One possible strategy would be, for an example, if you are looking 4 relevant sensors, to divide the interval T into 4 subintervals i T of equal length, and on each, selected the sensor with the best sensitivity (norm sense)

      2 2
,; the analysis of the functions of sensitivity (in the direction of descent k  ); see problem (5).

At each iteration of the minimization algorithm, it is easy to estimate the relevant spatial areas in time and select sensors.

The initialization of the 12 sensors amongst 16 (see iteration 1) is made to cover the plate without a priori.

In figure 7 This concept can be used in strategies of mobile sensors deployment.

III. CONCLUSION AND OUTLOOKS

A

For

  quasi online identification, the problem of optimum selection of sensors can be posed. Two strategies are discussed in this paper: the first concerns the selection quasionline of most relevant localized sensors within a network of sensors, the second works on the problem of the intelligent movement of some mobile sensors. In this communication, the complete description of the experimental device is presented. The choice of dimensions and materials are justified. Heat sources are detailed and an identification procedure of heat flux density is proposed. The acquisition of collected data is presented. Heat elements and sensors are embedded on mobile robots (Khepera III) and the procedure is detailed. The methodology for the identification is presented and the innovative aspects (procedure quasionline, strategies of sensors choice) are explained. Finally, conclusions and outlooks are briefly trained.

,Fig. 2 .Fig. 3 .

 23 Fig. 2. The heating sources

Fig. 4 .

 4 Fig. 4. Camera and Khepera-III robot equipped with a pyrometer For each robot equipped with a pyrometer, the temperature is measured and the analogical values (currents) are processed by a digital analogical converter (integrated in I/O card KoreIO connected on the expansion card). The robot then sends a frame to the management computer (by using the protocol TCP/IP client/server) containing the position and the temperature recorded at that position and the time and date of acquisition of the measure.

Fig. 5 . 6 .Fig. 6 .

 566 Fig. 5. Vision system prototype

0 

 0 in a frame of a vector space Oxy . Once the new coordinates calculated, they are compared with the odometric measurements of the robots in order to estimate as well as possible the coordinates of the best robots on the plate.After measure treatments and process of quasi online identification (see paragraph III), each robot receives a frame of data which contains its current coordinates   any). With the help of the SSL-vision software, the positions of the robots are updated every second to correct errors and avoid increasing the inaccuracy of the temperature measurements.



  Adjoint problem (4) written with Lagrangian formulation to compute the gradient of the cost function and the descent direction at iteration k .  Sensitivity problem (5): calculation of the sensitivity function defined as the variation of temperature induced by variation of the heat flux density (in the descent direction).

  of this method is proposed below for identification the heat flux density of a mobile heating source moving in thin plaque (1m²) with an offline CGM. Noisy measurements are available from 16 fixed sensors. In order to reduce the influence of noise measurement on the cost function stop J (6), 4 least relevant sensors are rejected during the iterative process based on norm   2 i L T . In Fig 7, is presented, 12 sensors (at different iterations) used to estimate the unknown heat flux density.

Fig. 7 :

 7 Fig. 7: Choice of sensors with offline method. The choice of the most relevant sensors during the iterative process of minimization is not based on priori information. It is based on the results of the identification and in particular

Fig. 8 . 3 5

 83 Fig. 8. Heat flux density identification at iteration 23.

  complete set (experimentation, definition of the inverse problem, resolution offline and online) is presented in this communication to examine different strategies for mobile sources tracking observation. Where the State is described by a system of parabolic PDE (thermal), the selection of the relevant sensors within a fixed network or the definition online of the trajectories of mobile sensors is discussed. The prospects envisaged at the end of this work consist in the realization of many experimental campaigns and to the study of non-linear systems for which input parameters would depend on the State. This work was partly funded with ANR-12-BS03-008-03. gradient algorithm for thermal diffusivity identification in a moving boundaries system. Journal of physics: Conf. Series, vol. 135, 2008. [12] Beddiaf S., Perez L., Autrique L. et Jolly J.C. Simultaneous determination of time-varying strength and location of a heating source

  , ;

	x y z t	T	c	t
	x y z			x y
	x y z t	lat	T	n
	x y z t	sup	T	n	e
	x y z t	inf	T	p n	R