Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
95,574 characters · 33 sections · 73 citation commands
Composite Likelihood for Stochastic Migration Model with Unobserved Factor
\affil[1]{York University} \affil[2]{Toronto University, Toulouse School of Economics and CREST}\affil[3]{Equitable (EQ) Bank, Toronto} \thispagestyle{empty}
\setcounter{page}{1}
Under the “internal-ratings-based” (IRB) approach advocated in the Basel II and III regulation, banks use their internal risk rating systems to estimate the risk exposures, credit rating migration probabilities and the probability of default (PD) in order to evaluate their regulatory capital requirements [see bas2,bas3, hull2012risk, GG2016]. Under Pillar II, financial institutions must also conduct stress tests to determine the level of capital needed to absorb losses in worsening economic conditions and be protected against systemic risk. For these reasons, banks perform their own credit rating migration analysis in order to monitor the changes in borrowers’ credit quality and to predict borrowers’ potential default in a volatile economic environment. This analysis concerns the “internal” or “in-house” established credit rating histories of borrowers, classified into credit quality categories, which are determined independently of the ratings publicly provided by the rating agencies such as the Moody's.\footnote{Publicly available credit ratings of large obligors are available from the rating agencies such as the Moody's, Standard and Poor's (S&P), and Fitch.} The internal credit rating analysis is applied to the historical probabilities of default and migration probabilities. It differs from the analysis of their risk-neutral counterparts, which underlies the pricing of credit derivatives, such as credit default swaps (CDS), Collaterized Debt Obligations (CDO), or derivatives written on iTraxx [see, e.g. DEHS2009, AGS2018 in continuous time, GMP2006 in discrete time, GMMR2021 for joint historical and risk-neutral analysis]. The internal ratings are used for pricing the portfolios of credits offered to a large number of small and medium-size firms whose assets are not traded on the markets. Even for large firms, the historical and risk-neutral probabilities of default can differ significantly. In our paper, the analysis of internal ratings is consistent with prudential banking supervision and aims at avoiding a pure mark-to-market pricing of risk.\footnote{There is often a confusion about the notions of historical and risk-neutral risks. For example, Moody's Analytics provides “EDF” estimates of the historical probability of default by considering default frequencies of firms with the same distance-to-default (DD). However, the notion of DD is risk-neutral.}
The credit rating migration analysis concerns the changes [i.e. upgrades or downgrades] of borrowers' credit quality over time with respect to their previous ratings [altman1997credit]. These data are available from monthly or quarterly time series of credit migration matrices comprising the qualitative ratings of firms, ranked from the low risk category A to the most risky rating D of default. The ordered Probit model for credit ratings arises as a natural specification, which has been extended to the Asymptotic Single Risk Factor (ASRF) model by Vasicek1991 [see also Vasicek1987, NPV2001]. The ASFR is a stochastic factor probit model of default with an independent and identically distributed common random unobserved factor capturing the systemic risk effect. The factor is assumed to drive the parameters of a latent quantitative score function in the model, which is transformed into qualitative ratings. Due to the presence of the unobserved common factor, the observed rating histories are cross-sectionally dependent, which can explain default correlation. GG2005a, Feng2008, \ extended this setup to multiple credit rating categories with common systemic factors that can be serially correlated, in order to predict the future credit ratings of firms. This extension is strongly recommended under the Basel III regulatory measures: “Interdependence between issuers is frequently modelled in a similar way to the regulatory framework, using a combination of an idiosyncratic (i.e. individual) and one or more systemic risk factors” [EBA2012, article 12 on Systemic Risk Factor]. Moreover, the dynamic ordered probit model takes into account the heterogeneity of issuers and satisfies the requirement that “Separate transition matrices may be applied for specific groups of issuers and geographical areas” [EBA2012]. It also reproduces other stylized facts such as the rating momentum [see e.g. AK1992].
The estimation of the ordered Probit model with a latent common factor is challenging. In order to derive the joint density of observed ratings, the history of the latent factor has to be integrated out. Therefore, the exact likelihood function based on the joint density of rating histories involves an integral of high dimension, increasing with the number of observations over time. Due to the presence of the multiple integrals, the exact maximum likelihood needs to be replaced by an approximation in practice. This paper introduces the conditional Maximum Composite Likelihood (MCL) as an alternative estimation method for the stochastic factor ordered Probit model. The MCL estimators have been widely used in the statistical literature to handle complex likelihood functions [see, Lindsay1988, Varian2008, VRF2011, GM2018]. The conditional composite likelihood functions are obtained by multiplying a collection of conditional component likelihoods, each depending on some integrals. In the one factor framework, these integrals are of dimension 1.
There exist alternative approximation methods, most of which involve a set of arbitrary control parameters, having a significant impact on the associated required capital. These parameters are, for example, the discretization steps [Farmer2021], tuning parameters, penalties, etc. The effect of the statistical approximation and optimization method can go as far as to partly eliminate the need for keeping an internal capital reserve, which is called a “statistical regulatory arbitrage”. Therefore these approximations are often not validated by the supervisory authorities who are regularly auditing the internal databases and estimation techniques.\footnote{“Any estimation technique should be duly justified and documented” [EBA2012].}
So far, the banking supervisory authority has validated selected standardized approximation methods, such as the granularity adjusted approach, that is valid and efficient if both the cross-sectional and temporal dimensions are large [see GG2014,GG2015, for general discussion] and the Simulated Maximum Likelihood (SML) method with a large number of simulations [Feng2008]. Both these methods circumvent the high-dimensional integration. Under the SML estimation employed in Feng2008, the integral is approximated by simulations, allowing for the latent factor values to be filtered out ex-post. The quality of the simulation-based approximation depends on the number of simulations, which can become high, depending on the number of time units considered and the complexity of factor dynamics. This makes this method computationally intense. The granularity-based approach [GG2015] is a two-step estimation method that eliminates the burden of simulations and provides the estimates of the unknown parameters and unobserved factor values. However, the granularity-based estimator depends in the first step on a set of “nuisance” parameters of size $T-1$, which increases the computational complexity of this method. We show that the proposed MCL estimators are computationally less intense than the granularity-based approach, and are reliable in finite sample.
In the panel analysis of credit ratings, the number of firms determines the cross-sectional dimension $n$, and the number of observed time units determines the dimension $T$. Two different asymptotics are considered, when $n,T$ both tend to infinity, or $n \rightarrow \infty$ and $T$ fixed, i.e. finite sample in $T$. In practice, $n$ is often large, while $T$ can be rather small. Therefore, these two types of asymptotics are compared. When both $n$ and $T$ tend to infinity, the new conditional composite maximum likelihood estimators are shown to be consistent, but not fully efficient, while the granularity-based estimator is consistent and asymptotically efficient. However, when $n$ is large and $T$ is fixed, all estimators converge to stochastic limits that depend on the latent factor values and differ from the true values of the parameters.
This paper is organized as follows. Section 2 compares the credit rating models that exist in the literature and describes the ordered probit model of credit rating transitions. Section 3 introduces the conditional composite maximum likelihood estimators and the granularity approach. The order and rank conditions for identification are also provided. Section 4 derives the asymptotic properties, i.e. the consistency, rates of convergence and asymptotic normality when both $n$ and $T$ tend to infinity, and also when $n$ tends to infinity and $T$ is fixed. In Section 5, the performance of MCL and granularity-based estimators in finite sample is examined in a simulation study. Section 6 includes the empirical application. The observed transition probabilities are computed from the Compustat Standard and Poor's (S&P) rating database from 1985Q4 to 2016Q4, available through Wharton Research Data Services. We analyze the estimated parameters, probabilities of defaults, and the downgrade probabilities at different horizons. Section 7 concludes the paper. Proofs are given in Appendices A-C and the simulation details and additional simulation results are presented in the online Appendix D. Throughout the paper, variables indexed by $t$, $i$ in parenthesis denote random sequences.
In this section, we focus on the stochastic factor ordered probit model and its state space representation. The expression of the complete likelihood function is derived, highlighting the presence of multiple integrals of large dimension. Next, we discuss the models of joint evolution of individual ratings that already exist in the literature.
Let $y_{i,t}$ denote the rating of firms $i, \; i=1,...,n$ at time $t, \; t=1,...,T$. The ratings are qualitative variables that take $K$ values associated with different rating categories. The sequences of variables $y_{i,t} \; t=1,...,T $ for $i=1,...,n$ represent the panel of qualitative individual histories of credit ratings. The migration model defines the joint distribution of the qualitative variables $y_{i,t}, \; i=1,...,n, \; t=1,...,T$ and provides information on the transitions (migrations) of individuals (firms) between the ratings.
Let us consider the following extention of the standard Vasicek model of default risk [Vasicek1991, GL2013, GG2016] to a migration model with a latent factor. The unobserved (latent) stochastic common factor is denoted by $f_t$.
We assume that the conditional transition probability of the factor given the whole past information depends on $f_{t-1}$ only , i.e.,
$$l(f_t|f_{t-1}, f_{t-2},... ;\;y_{i, t-1}, y_{i, t-2},...; i=1, ...,n) = l(f_t |f_{t-1}).$$
In addition, we assume conditional on the path of the common factor, the individual rating histories are independent, heterogeneous Markov chains with
$$P[y_{i,t} = k| y_{i,t-1} = j, f_t] = p_{jk}(f_t).$$
Under the above assumptions, the joint process $(y_{i,t}, i=1,...,n, f_t)$ is a Markov process with an exogenous evolution of factor process $(f_t)$.
Since the factor is unobserved, its evolution has to be integrated out to get the joint distribution of individual histories. This creates migration (and default) correlation because the factor is common to all individuals (firms). The interdependence of risks has to be included as an incremental risk, accounted for by additional required capital [see bas3, EBA2012]. It also implies non-Markovian features after integrating the conditional transition probability with respect to factor $f = (f_t)$. In particular, the conditional transition probability
$P[y_{i,t+1} = k, y_{i,t+2} = k| y_{i,t} = j] $, for example, is not equal to $P[y_{i,t+2} = k| y_{i,t+1} = k] P [y_{i,t+1} = k| y_{i,t} = j] $.
Indeed, we have $$P[y_{i,t+1} = k, y_{i,t+2} = k| y_{i,t} = j] = \frac{ E P[y_{i,t} = k, y_{i,t+1} = k, y_{i,t+2} = j|f] } {E P[y_{i,t} = j| f]},$$
where the expectation is taken with respect to the stochastic evolution of $f$ over the period $(0,t+2)$, which has a different impact on the probability of staying in state $k$ depending on the last transition being an up- or down-grade, and the date of that transition.
Therefore, the model account for the rating momentum effect, i.e. the fact that the intensity of transitions out of a given state is influenced by previous transitions into that state, and more generally for non-Markovian features [GK2009].
This unobserved factor model can be viewed as an infinite mixture model at time $t$, with stochastic weights. An example of this type of model is the stochastic factor ordered probit model [GG2005a, GG2014, Feng2008, HZ2015, CLP2021] examined in this paper.\footnote{or its continuous time counterparts, i.e. the dynamic marked point processes with common systemic factors (see Creal2013, Section 4.3, Koopman2008 for a continuous time approach without systemic factor).} Its state-space representation is given below.
Let $y_{i,t}^*$ and $y_{i,t}$ denote the (credit) score and rating of firm $i$, $ i = 1, . . . , N$ at time $t$, $t = 1, \ldots, T$. The latent continuous quantitative score $(y^{*}_{i,t})$ determines the individual qualitative rating $y_{i,t}$. More precisely, the quantitative score is discretized in order to obtain the individual qualitative ratings. Therefore, an observed rating is determined as follows,
where $c_1< \cdots <c_{K+1}$ are the thresholds. Relation (2.1) shows how the observable endogenous credit rating $(y_{i,t})$ is linked to the latent score function $(y_{i,t}^{*})$. By convention, we have $c_{1}=-\infty$ and $c_{K+1}= +\infty$. Relation (2.1) defines the measurement equation of the state space representation of the model.
The conditional distribution of the quantitative scores given the factor path and the previous scores depends on the common latent factor $f_t$ and on the past individual ratings $y_{i,t-1}$, such that
and $y_{i,1}$ is the first observed rating for firm $i$. The multivariate, continuous, latent processes $y^{*}_{i,t}$, are generated by individual level effects $(\delta_{j})$, volatility effects $(\sigma_{j}), \sigma_j>0$, factor effects where the components of $\beta_{j}$ define the factor sensitivities. When coefficient $\beta$ is large (small, resp.), the effect of systemic risk carried through the factor is strong (weak, resp.). All the parameters $\delta_j$, $\beta_j$, $\sigma_j$ depend on the previous rating $j$. While the idiosyncratic risks $(u_{i,t})$ can be diversified, the systemic risk $(f_{t})$ cannot be diversified. Thus the presence of systemic risk generates risk interdependence in the model. Because parameters $\beta$ are different in each rating category, the risk interdependence varies across rating transitions resulting in risk momentum. Among these parameters, $\delta_j$ and $\sigma_j$ summarize the effect of idiosyncratic risk, and $\beta_j$ is the sensitivity to systemic risk.
The following autoregressive model of order 1 (AR(1)) represents the common factor dynamics,
where $\eta_t$ defines the shock to the common factor and $f_1$ is drawn in the stationary distribution. Alternatively, a multidimensional factor can be considered to distinguish between the dynamic migration patterns of firms with good and poor credit quality, respectively. See, for example, GG2014, for more details.
The system of equations (2.2)-(2.3) defines the state equations of the state-space model. Let us introduce the following assumptions to obtain a migration model with migration correlation and rating momentum.
{\bf Assumption A.1:} The errors $u_{i,t}, \; \eta_t, \; i=1,...,n, t=1,...,T$, are independent, standard normal variables.
The independence assumption allows for performing impulse response analysis by shocking separately the idiosyncratic and systematic innovations, to perform a stress-test in particular. The assumption of identical distribution and the fact that coefficients in (2.2) are independent of the firm implies that we consider a homogeneous set of firms, obtained by crossing the country, industrial sector and firm size, in compliance with the current regulation.
{\bf Assumption A.2:} The factor process $(f_t)$ is the strongly stationary solution of autoregressive equation (2.3).
As the processes ($f_t$), ($u_{i,t}$), $i=1,\ldots,n$, are independent and strictly stationary, it follows that the joint $n-$dimensional process $y_t^\ast=\left(y_{1,t}^\ast,\ldots,y_{n,t}^\ast\right)^\prime$ is also strictly stationary, and so is its state discretized version $y_t=\left(y_{1,t},\ldots,y_{n,t}\right)^\prime$. However, the individual components ($y^\ast_{i,t}$), $i=1,\ldots,n$ are not independent due to the effect of the common factor $f_t$.\footnote{In this respect this model differs from Tuzcuoglu2019, where the state equations (2.2)-(2.3) are replaced by $ y_{i,t}^* = \rho y_{i,t-1}^* + \beta' x_{i,t} + \alpha_i + \epsilon_{i,t},$ with independent $(\alpha_i, (\epsilon_{i,t})),\; i=1,...,n,$ given $x$. This specification does not contain systemic risk and does not allow for risk interdependence.}
It is important to notice that the error variance in equation (2.3) has been set equal to $1-\rho^2$. This implies that factor $f_t$ is marginally normally distributed with mean 0 and variance 1: $E(f_t)=0, Var(f_t) =1$. These moment restrictions are introduced to solve the factor identification issue, since the factor is defined up to a linear affine transformation.
In practice, the underlying quantitative scores are computed by a credit institution and each individual (firm) can request the records of its own score history. However, the complete score database is, in general, proprietary and the information on the quantitative scores is not available to an outsider econometrician/data scientist. The factor $f_t$ is assumed unobserved for the following two reasons: First it creates the cross-sectional correlation between individual risks. Second, it provides a dynamic model that can be used to predict the future defaults. A bias could result from directly replacing factor $f_t$ by an observed proxy $\hat{f}_t$, such as the VIX market volatility index, a consumer sentiment index, consumption growth, a business cycle indicator [see e.g. BDDF2018, AGS2018], or the slope of the yield curve. Moreover, if factors are observed, their predictions cannot be computed without specifying an additional model of the dynamics for all the observed factors in $\hat{f}_t$, and checking that these observed factors are exogenous.
In order to derive the joint density of observations $y_{i,t}, \; i=1,...,n, \; t=1,...,T$, the unobserved factor path has to be integrated out. As a consequence, observations $y_{i,t}$ are cross-sectionally dependent and serially dependent with a non-Markovian serial dependence. More precisely, the stochastic migration probabilities between dates $t-1$ and $t$, conditional on $f_t$, are given by
where $\Phi$ denotes the cumulative distribution function (c.d.f.) of the standard normal. Thus each row of the transition matrix conditional on $(f_t)$ contains an ordered polytomous probit model with a common explanatory factor $f_{t}$. When factor $f_{t}$ is unobserved stochastic and serially correlated as in (2.3), the transition matrices are stochastic and serially dependent.
Let us now define the log-likelihood function of the stochastic migration model. The vector $\theta$ includes the parameters of the state space model, which are parameters $\delta_j,\beta_j, \sigma_j,j=1,\ldots,K$ in the quantitative score, and parameters $c_k, k=2,\ldots,K$ defining the states. As the conditional migration matrices are functions of parameter vector $\theta$ as well as of the common factor values $f=(f_{t})$, the likelihood function conditional on $f$ and the initial rating $y_1$ is
where $n_{jk,t}$ denotes the number of firms which migrate from $j$ to $k$ between $t-1$ and $t$, $Y=(y_{i,t})$ for $i=1, ..., n$ and $t=2, ..., T$, $f=(f_t,\; t=2,\ldots,T)$, and $y_1=\left(y_{1,1},\ldots,y_{n,1}\right)^\prime$.
Since the factor history is not observed, after integrating out the factor values $(f_{2}, ..., f_{T})$, the log-likelihood function, given the initial value $y_1$ only, is
where $\psi$ refers to the joint probability distribution function of factor values. The above log-likelihood function contains a multivariate integral. The dimension of this integral is of order $T-1$, as there is a common factor value for each transition at time $t$. Therefore the exact computation of this likelihood is infeasible and its approximation is often not sufficiently robust.\footnote{See Feng2008 for the discussion of robustness when simulations are used.} The MCL estimators are convenient alternatives for complicated nonlinear dynamic state-space models allowing for circumventing the high-dimensional integral.
In the literature, there exist alternative models of rating transitions that rely on simplified assumptions. For comparison, these models can be classified as the Markov chain models and models with a rating momentum adjustment.
The basic migration model assumes that the histories of individual transitions (migrations) of firms between the ratings (states) are independent, identically distributed Markov chains. The assumption of independent Markov chains implies that the individual histories ($y_{i,t}, \; t=1,...,T$) are independent across individuals (firms) $i$, $i = 1,...,n$. These Markov chains are homogeneous with transition probabilities $P[y_{i,t} = k| y_{i,t-1} =j] = p_{jk}$ independent of the firm.\footnote{We choose the transition matrix $P=[p_{jk}]$ such that the elements of each row sum up to one.} This model is examined in LS2002, Section 3, and reviewed in DPS2020, for example. It does not account for the rating momentum and it does not assume migration correlations among the individuals, implying no default correlation among the firms.
The rating momentum is accounted for in a non-Markovian model, such as the hazard model introduced in LS2002 and DPS2020, Section 4, based on the assumption of Semi-Markov chains. This assumption implies that the individual histories are independent. Then, the transition probabilities are given by $$P[y_{i,t} = k| y_{i,t-1}, y_{i,t-2},,...] = p_{jk} exp( c Z_{i, t-1}),$$ where $y_{i,t-1}=j$ and $Z_{i,t-1} = 1$, if firm $i$ was downgraded to its current state, $Z_{i,t-1} = 0$, if firm $i$ was upgraded to its current state.
Alternatively, the non-Markovian feature can be introduced by considering a mixture of Markov chains [FS2008]. Under the assumption of Markov Chain mixture, the individual histories are independent and the transition probabilities are obtained from a mixture of Markov chains.
Other non-Markovian assumptions can be introduced by considering time varying exogenous variables, such as the observed regime of business cycle [BDKSS2002, GS2014], or time itself. These extensions are characterized by the independence of migrations, implying no migration correlation among the individuals (firms). Hence, drawback b) is the common issue of all these models.
Other models without a latent factor can be found in RC2017, MT2012, and HVH2022. Alternatively to the frequentist approach in this paper, Bayesian procedures like the one in STT2009 could be used if the researcher is willing to impose priors on the different parameters.
The process of transition matrices $\{P_{t}, t=1, ..., T\}$ has component matrices $P_{t}=\left(p_{jk,t}\right)$, which provide the probabilities of transitions from state $j$ to state $k$ between times $t-1$ and $t$ given $f_t$. From (2.4), it follows that the elements of matrix $P_{t}$ are $$p_{jk,t}=p_{jk}(f_{t};\theta)= \mathbb{P}[y_{i,t}=k|y_{i,t-1}=j,f_{t}] =\Phi\left(\frac{c_{k+1}-\beta_{j}f_{t}-\delta_{j}}{\sigma_{j}}\right)-\Phi\left(\frac{c_{k}-\beta_{j}f_{t}-\delta_{j}}{\sigma_{j}}\right),$$ $\; k, j =1,...,K .$
Let us now compute the product of two successive transition matrices $P_{t}^{(2)}=P_{t}P_{t-1}$ to obtain the probabilities of transition at horizon 2 from state $j$ to $k$ between times $t-2$ and $t$ given the factor history $\left(f_t\right)$. The elements of matrix $P_{t}^{(2)}$ depend on $f_t, f_{t-1}$ and are given by
They can be computed from the elements of matrices $P_{t}$ and $P_{t-1}$. Let us denote by $P$ and $P^{(2)}$ the expectations of matrices $P_{t}$ and $P_{t}^{(2)}$ with respect to the common factor history, i.e,
The elements of matrix $$ P=[p_{jk}]=[p_{jk}(\theta)]=E_{f_{t}}[p_{jk}(f_{t},\theta)],$$ are obtained by integrating out the unobserved factor value $f_t$.
Proof. See Appendix A.1.
It is easy to see that the expected matrix $P$ is also a transition matrix. Indeed, its elements are non-negative and additionally $Pe = E(P_t)e = E (P_t e) = e$, where $e$ is a vector of ones. Matrix $P$ is a quasi-transition matrix for $y_{i,t}$, since $p_{jk}$ is not equal to the conditional probability of $y_{i,t }=k$ given $y_{i,t-1}=j$, except when the $f_t$'s are i.i.d., i.e., when $\rho=0$. In fact, it is computed as if $f_t$ was independent of $y_{i,t-1}$. Each row of this quasi-transition matrix corresponds to another ordered probit model.
The elements of matrix $P(2)$ are obtained by integrating jointly with respect to $f_t, f_{t-1}$. We have
Proof. See Appendix A.2.
$P^{(2)}(\theta, \rho)$ is a quasi-transition matrix at horizon 2 computed as if $(f_t, f_{t-1})$ were independent of $y_{i, t-1}$. The quasi-transitions at horizon 2 involve one-dimensional integrals only, which are easy to compute numerically.
This section presents the conditional composite likelihood functions for the migration model with an unobserved AR(1) factor. The composite likelihoods are often based on misspecified likelihoods, which are easier to calculate [see CR2004, VRF2011]. In our framework, the conditional composite likelihoods are constructed from the quasi-migration probabilities at horizons 1 and 2 to reduce the dimension of the integrals. We also present the conditional likelihood used in the first step of the granularity approach.
As mentioned earlier, the parameter vector $\theta=(\beta_k, \delta_k, \sigma_k, c_k)$ includes the parameters characterizing the latent quantitative score, representing the systemic and idiosyncratic risks, and the thresholds that define the qualitative rating category associated to the latent quantitative score. The additional parameter $\rho$ allows for predicting the future systemic risk. Some among the estimation methods given in this section are focused on the rating parameters $\theta$, while others concern both $\theta$ and serial dependence parameter $\rho$.
i) The Conditional Composite Log-Likelihood at Lag 1
The conditional composite log-likelihood function at lag 1, called CL(1), is focused on parameter $\theta$. The associated log-likelihood $L_{cc}(\theta)$ is defined as
where $\hat{p}_{jk,t}= n_{jk,t}/n_{j, t-1}$ is the observed transition frequency from $j$ to $k$ in one step over the period $(t-1,t)$, $n_{j, t-1}$ is the count of firms with rating $j$ at the beginning of period $t$, and $\pi_j, \; j=1,...,K$ is a given set of weights. The log-likelihood $L_{cc}$ is calculated as if the observed ratings $(y_{i,t})$, $i=1, ..., n$, were independent across the individuals, while in reality they are linked by the common factor. Moreover, $L_{cc}$ considers the rating processes $(y_{i,t})$, $i=1, ..., n$, as if these were components of a Markov chain with quasi-transition matrix $P$, although $(y_{i,t})$, $i=1, ..., n$, are not Markov because integrating the factor out increases the memory of the process. It also assumes a time independent rating structure $(\pi_j, \; j=1,...,K)$. Therefore, the CL(1) is a quasi (pseudo) log-likelihood. The conditional composite log-likelihood CL(1) depends on parameter vector $\theta$ only, and cannot be used to estimate the factor dynamics, i.e, the autoregressive coefficient $\rho$. For that purpose, it is necessary to increase the lag.
ii) The Conditional Composite Log-Likelihood at Lag (2)
The conditional composite log-likelihood at lag (2), called CL(2), depends on both parameters $\theta$ and $\rho$. The log-likelihood $L_{cc,2}(\theta,\rho)$ is given by
where $\hat{p}_{jk,t}^{(2)}$ is the observed transition frequency from state $j$ to $k$ in two steps over the period $(t-2,t)$ and $\pi=(\pi_j,\; j=1,...,K)$ is a fixed structure of ratings.
The composite log-likelihood function $L_{cc,2}(\theta,\rho)$ is computed from the density of $(y_{i,t})$ conditional on $(y_{i,t-2})$ as if the rating histories $(y_{i,t})$ were cross-sectionally independent from one another, $(y_{i,t-2})$ were containing all information about the past and were based on quasi-transitions over 2 steps. Therefore, the CL(2) is a quasi (pseudo) log-likelihood too.
An important difference between $L_{cc}$ and $L_{cc,2}$ is the set of identifiable parameters. As mentioned above, we can expect to identify $\theta$ from $L_{cc}$, but we cannot identify parameter $\rho$ characterizing the cross-sectional dependence. $L_{cc,2}$ provides additional information that is sufficient to identify $\rho$.
iii) The Conditional Composite Likelihood up to Lag 2
The conditional composite log-likelihood up to lag 2, CL(1,2), is defined by summing up the previous composite log-likelihoods at lags 1 and 2,
where $a$ is a constant to be selected.\footnote{$a$ could be optimally selected to increase the efficiency of the estimator in a two step approach [Cox, CR2004, p.730].} It concerns both parameters $\theta$ and $\rho$. This objective function cannot be interpreted as a quasi-likelihood.
iv) The Granularity-Based Conditional Log-likelihood
Let us now introduce another type of log-likelihood for the estimation of parameter $\theta$. As shown in Section 2.2, the complete log-likelihood has a complicated expression including a high-dimensional integral of a dimension increasing with $T$. The granularity approach [GG2005a,GG2015] replaces the complete log-likelihood by an appropriate expansion for large $T$. This leads to a two step estimation method where, in the first step, the factor values are considered as fixed time effects. This log-likelihood conditional on $(f_2,...,f_T)$ is
where $n_{jk,t}$ (resp. $n_{j,t-1}$) counts all transitions from $j$ to $k$ (resp. is the structure of ratings at $t-1$). It is maximized with respect to both parameter $\theta$ and factor path $f_2,...,f_T$ subject to the identification restrictions
These restrictions on time fixed effects $(f_t)$ approximate the identification restrictions on the latent stochastic factor, i.e. $E (f_t) =0, Var(f_t) =1$. This conditional constrained log-likelihood resembles the composite log-likelihood $L_{cc}$ except that in the composite log-likelihood, $p_{lk}(\theta)$ was made independent of $f_t$ by marginalizing and the observations are introduced with a fixed rating structure $\pi=(\pi_j, j=1,....,K)$. Since this objective function is maximized with respect to $\theta, f_2, \ldots ,f_T$, it provides not only an estimator of $\theta$, but also an approximation $\hat{f}_t$ of the factor values. In the second step, an estimator of $\rho$ is obtained by regressing $\hat{f}_t$ on $\hat{f}_{t-1},t=2,\ldots,T$.
Let us focus on parameter $\theta$ and briefly discuss the expected properties of the above estimation methods. In the panel framework involving both $n$ and $T$, various notions of asymptotics can be considered. When $n$ and $T$ both tend to infinity, the granularity approach provides consistent and asymptotically efficient estimators [GG2014]. In Section 4, we prove that both CL(1) and CL(2) methods also provide consistent estimators of $\theta$. In practice, the cross-sectional dimension $n$ is large, but $T$ is much smaller. Therefore, we expect finite sample effects in $T$ affecting all the estimators. When $n$ tends to infinity, $T$ is fixed, all estimators converge to pseudo-true values $\theta_{\infty}(f_{0,2},...,f_{0,T})$ depending on the latent factor values $f_{0,2},...,f_{0,T}$, including the first step of the granularity approach due to replacing the factor identification restrictions $E (f_t) =0, Var(f_t) =1$ by their fixed effect counterparts. Hence, for $T$ fixed, all the estimators considered are asymptotically (in $n$) biased both conditionally on factor values and after re-integrating the factor. Let us now discuss their variances conditional on $f_{0,2},...,f_{0,T}$. In this respect, it is important to consider the number of “nuisance” parameters in each of the estimation method: no nuisance parameter in CL(1), one nuisance parameter $\rho$ in CL(2), $T-1$ nuisance parameters $f_{0,2},...,f_{0,T}$ in the first step of the granularity approach.\footnote{The number of nuisance parameters quickly increases if more systemic risk factors are introduced.} The variances can increase with the number of nuisance parameters. Moreover, the bias and variance trade-off\footnote{In statistics, the variance-bias trade-off is through the quadratic loss = variance + squared bias. In credit portfolios, it is through a Value-at-Risk of the type VaR= bias + 1.96 $\sqrt{\mbox{variance}}$.} would depend on the dynamic pattern of the true factor values $f_{0,2},...,f_{0,T}$, in particular if they approximately satisfy the restrictions of zero sample mean and unit sample variance for the granularity approach, or are more or less erratic for the CL(1) and CL(2) methods.
In this section, the order and rank conditions for identification of each of the conditional composite log-likelihoods are discussed. The identification of $\theta, \rho$ in the conditional granularity approach has already been examined in GG2005a, GG2015.
The parameters to be identified and their respective numbers are as follows,
The total number of independent parameters to identify is $4K-2$. The negative two is due to the score $y_{i,t}^{*}$ being defined up to an increasing function. As we have supposed that it was a linear function of factor $f_{t}$, the score $y_{i,t}^{*}$ is defined up to a linear affine increasing function. The intercept and slope of that linear function are not identifiable.
In this subsection, the order conditions for each conditional composite log-likelihood are discussed. These conditions are derived from the probabilities $p_{jk}(\theta), p_{jk}^{(2)} (\theta, \rho)$ that appear in the composite log-likelihoods. These probabilities can be consistently estimated if $n$ and $T$ tend to infinity\footnote{They cannot be consistently estimated otherwise, in particular when $n\rightarrow \infty$, $T$ fixed. Indeed, in such a panel framework, the identification does not necessarily imply the existence of a convergent estimator.} (see Section 4).
i) Identification of $\theta$ under CL(1)
The identifying functions are the reduced form parameters in the CL(1) objective function, i.e. the elements $p_{jk}(\theta)$ of the quasi-transition matrix $P$. There are $K(K-1)$ of these elements that are linearly independent because of the unit mass restriction on each column. Hence, the order condition is
by taking into account the absence of parameter $\rho$ in the objective function. This order condition is satisfied for $K\geqslant 5$.
ii) Identification of $\theta$ under CL(2)
The identifying functions are determined by observing that the factor $f$ varies within the integral expression of $p_{jk}^{(2)} (\theta, \rho)$ (see Lemma 2). These identifying functions and their respective numbers are as follows:
where $\epsilon=\pm 1$ is an unknown sign, since the distribution of $f$ is symmetric. This implies that the integral expression in Lemma 2 is also valid with $f$ replaced by $-f$. There is only one such invariance property and therefore the sign $\epsilon$ is equal for all $j$. The total number of identifying functions of parameters is $2K(K-1)+2K=2K^2$. Hence, the order condition is
The order condition holds for any $K$.
iii) Identification of $\theta, \rho$ under CL(1,2)
The total number of functions available is equal to the sum of functions available for each component of the total composite log-likelihood. Therefore, the order condition is
The order condition is satisfied for any $K$.
The rank conditions are important for the local identifiability. They are derived for the CL(1) and CL(2) approaches and are similar to the rank condition derived for the granularity approach in GG2005a and GG2015, p.84.
Proof. See Appendix B.1.
Proof. See Appendix B.2.
In order to identify the unknown sign $\epsilon$, an additional constraint needs to be introduced such as $$\beta_{1}>0.$$ The unknown sign $\epsilon$ is a problem of global identification and not of local identification. Hence, when the asymptotic properties of the estimators are derived (see Section 4), this positivity constraint has to be taken into account to obtain the consistency of the estimator. It has no effect on the asymptotic normality. The asymptotic properties of the composite log-likelihood estimators are discussed in the next section.
In a panel data framework, the asymptotic analysis can be performed with respect to the cross-sectional dimension $n$ and time dimension $T$ that can tend to infinity as follows.\footnote{The last case (iii) $ n$ fixed, $T \rightarrow{\infty}$ is less relevant for applications to credit rating.}
The double asymptotics in case (i) has been developed for applications to big data [GG2015, GG2014, BJR2017]. It corresponds to a long panel of high dimensional time series.
When $n \rightarrow \infty, T \rightarrow \infty$, the asymptotic properties of the conditional maximum composite likelihood estimators are much easier to derive than the asymptotic properties of the complete ML estimator. Indeed, the conditional composite log-likelihood functions are finite sums of products of summary statistics and functions of parameters. This simplifies the proof of uniform convergence with respect to the parameters. The next section examines the asymptotics (i) -(ii) and describes the properties of the conditional composite maximum likelihood estimators.
This section examines the consistency of the maximum conditional composite likelihood estimators of the identifiable parameters when $n \rightarrow \infty, T \rightarrow \infty$. To prove the consistency, we need the following additional assumption.
{\bf Assumption A3}
a) The parameter set of $(\theta, \rho)$ is compact, and strictly included in the set $\sigma_j>0, \forall j, |{\rho}|<1$.
b) The model is well-specified and the true value $(\theta_0, \rho_0)$ is in the interior of the parameter set.
The condition $\sigma_j>0, \; \forall j$, ensures that the transition probabilities $p_{jk}(f_t; \theta)$ [resp. $p_{jk} (\theta), p_{jk}^{(2)} (\theta, \rho)$] are infinitely continuously differentiable with respect to $f_t$ and $\theta$ (resp. with respect to $\theta, \rho$).
(i) Double asymptotics: $n\rightarrow{\infty}, T\rightarrow{\infty} $
Let us consider the double asymptotics with CL(1) approach. We have:
$$L_{cc}(\theta)=\sum_{k=1}^{K}\sum_{j=1}^{K}\sum_{t=2}^{T} \left[\pi_j \hat{p}_{jk,t} \ log \ p_{jk}(\theta)\right].$$
The conditional composite likelihood $L_{cc}$ depends on $n$ and $T$, although it is not indexed by $n$ and $T$ to simplify the notation. Since $T$ is varying, we need uniform a.s. convergence of the ratios $\hat{p}_{jk,t}$ to $p_{jk}(f_t; \theta_0)$ with respect to $t$, not only their pointwise a.s. convergence.
{\bf Assumption A.4}
i) $P [Max_{m \geq n} | \hat{p}_{jk,t}(m) - p_{jk}(f_t; \theta_0)| > \epsilon \,|f_t) < \frac{g_{jk} (f_t; \theta_0)}{ n \epsilon^2}$,
$\forall j,k,\epsilon, n, f_t$, where the notation $\hat{p}_{jk,t}(n)$ is introduced to indicate the dependence of the transition probabilities on the number of observations, and it exists a function $g_{jk}$ of $f_t$ and $\theta_0$, which is integrable with respect to the marginal distribution of $f_t$.
ii) $n,T \rightarrow \infty$ with $T/n \rightarrow 0$.
Assumption A.4 i) is a domination condition. Assumption A.4 ii) means that we have a panel with the cross-sectional dimension much larger than the time dimension. This allows for disregarding the uncertainty in $n$ with respect to the uncertainty in $T$ in the double asymptotic [see Appendix C].
Let us consider the objective function normalized by $T$, $$\frac{1}{T} L_{cc}(\theta) = \sum_{j=1}^K \pi_j \sum_{k=1}^K \left[ \left( \frac{1}{T} \sum_{t=2}^T \hat{p}_{jk,t}\right) \log \, p_{jk} (\theta) \right].$$ When $n$ tends to infinity, this quantity tends to \\ $\sum_{j=1}^K \pi_j \sum_{k=1}^K \left[ ( \frac{1}{T} \sum_{t=2}^T p_{jk}(f_t, \theta_0)) \log \, p_{jk} (\theta) \right].$ If moreover $T$ tends to infinity, the limiting objective function is
by using the ergodicity of the factor process and the Strong Law of Large Numbers in time dimension.
By the property of the Kullback-Leibler divergence measure applied to each row of the transition matrix, we know that the associated limiting conditional composite log-likelihood is maximized at $\theta_0^*$, with $$p_{jk}(\theta_0^*)=p_{jk}(\theta_{0}), \forall j,k.$$ Then, by the identifiability of $\theta=(c, \delta, \gamma)$ (see Proposition 1), we get $\theta^{*}_{0}=\theta_{0}$, and the consistency follows.
(ii) Short panel asymptotics: $n\rightarrow\infty,T$ fixed.
Like for the granularity approach, we cannot expect the conditional composite ML estimators to be consistent for $n \rightarrow \infty$, $T$ fixed. This is a consequence of the cross-sectional dependence due to the common systemic factor $f_{t}$. To clarify this point, let us assume $T=2$ and consider the maximum conditional composite likelihood CL(1) estimator. For $T=2$, the conditional composite log-likelihood is $$L_{cc}(c,\delta,\gamma)=\sum_{k=1}^{K}\sum_{j=1}^{K} \left[\pi_j \hat{p}_{jk,2} \ log \ p_{jk}(\theta)\right],$$ where $\theta = (c, \delta, \gamma)$ is the identifiable parameter satisfying the identification restriction in Proposition 1, that are $c_2=0,\gamma_1=1$. By Assumptions A.1, A.2 and the fact that the rating indicators are nonnegative and bounded, we can apply the Strong Law of Large Numbers to individuals. The conditional composite log-likelihood tends a.s. to $$\lim_{n\rightarrow{\infty}} \mbox{a.s.} L_{cc}(c,\delta,\gamma)=\lim_{n\rightarrow{\infty}} \mbox{a.s.} \sum_{j=1}^{K} \pi_j \Bigg( \sum_{k=1}^{K} \Bigg[ p_{jk} (\theta_0, f_{02}) \ log \ p_{jk}(\theta) \Bigg] \Bigg),$$
Then, this limiting objective function admits at least a maximum on the parameter set by Assumption A4 ii). Let $ \theta_0^*$ denote the pseudo-true value, i.e. a solution of the asymptotic optimization problem, we have $$ \theta_0^* = \underset{\theta}{\operatorname{argmax}}\sum_{j=1}^{K} \Bigg[ \pi_j \bigg[ \sum_{k=1}^{K} p_{jk}(\theta_0,f_{0,2}) \ log \ p_{jk}(\theta) \bigg] \Bigg]. $$ This pseudo-true value is a function of $\theta_{0}$ and $f_{0,2}$. Therefore, it cannot be equal to the true value, that does not depend on $f_{0,2}$. In other words, the MCL estimator $\hat{\theta}_n$ converges to a stochastic limit whose distribution depends on the distribution of $f_2$.
For expository purpose, we continue the discussion of the CL(1) approach for $n \rightarrow \infty, T \rightarrow \infty$. As mentioned above, the conditional composite log-likelihood is continuously differentiable. Since the estimator $\hat{\theta}_{n,T} = (\hat{c}_{n,T}, \hat{\delta}_{n,T}, \hat{\gamma}_{n,T})$ tends to the true value $\theta_0= (c_0, \delta_0, \gamma_0)$, which is in the interior of the parameter set, the estimator will also be asymptotically in the interior of the parameter set and will satisfy the necessary first-order conditions for large $T$. Therefore, we have $$ \frac{\partial L_{cc} (\hat{\theta}_{n,T})}{\partial \theta}=0 \iff \sum_{k=1}^{K}\sum_{j=1}^{K}\sum_{t=2}^{T} \left[ \pi_j \hat{p}_{jk,t} \ \frac{\partial \ log p_{jk}(\hat{\theta}_{n,T})}{\partial \theta}\right]=0. $$ We can perform a Taylor-McLaurin expansion with respect to $\hat{\theta}_{n,T}$ in the neighborhood of $\theta_{0}$. Let us assume:
{\bf Assumption A.5} The parameter set $\Theta$ for $\theta$ is convex.
We get
where $\tilde{\theta}_{n,T}$ is an intermediate value between $\hat{\theta}_{n,T}$ and $\theta_0$.
By applying the same argument as for the uniform a.s. convergence of the composite log-likelihood function, we deduce that
$\frac{1}{T} \sum_{k=1}^K \sum_{j=1}^K \sum_{t=2}^T \left[\pi_j \hat{p}_{jk,t} \frac{\partial^2 log \ p_{jk} (\tilde{\theta}_{n,T}) }{ \partial \theta \partial \theta'} \right]$ will converge a.s. to $\sum_{k=1}^K \sum_{j=1}^K \left[\pi_j p_{jk} (\theta_0) \frac{\partial^2 log \ p_{jk} (\theta_0) }{ \partial \theta \partial \theta'}\right]$,
$\frac{1}{T} \sum_{k=1}^K \sum_{j=1}^K \sum_{t=2}^T \left[\pi_j \hat{p}_{jk,t} \frac{\partial log p_{jk} (\theta_0) }{ \partial \theta}\right]$ will converge a.s. to
$ \sum_{k=1}^K \sum_{j=1}^K \left[\pi_j p_{jk} (\theta_0) \frac{\partial log \ p_{jk} (\theta_0) }{ \partial \theta}\right]=0$,
since $\theta_0$ is the maximizer of the limiting objective function (4.1), and
where $o_p(1)$ is a negligible term in probability by Assumption A.4. Let us assume:
{\bf Assumption A.6} The matrix $J_0 = \sum_{k=1}^K \sum_{j=1}^K \left[ \pi_j p_{jk} (\theta_0) \frac{\partial^2 log \ p_{jk} (\theta_0) }{ \partial \theta \partial \theta'} \right]$ is positive definite.
Then, by normalizing the expansion (4.2) by $1/( \sqrt{T})$, we get
where $vec$ denotes the vectorization that stacks the columns of the transition matrix. Note that $$ vec [\pi_j p_{jk} (f_t; \theta_0)] = vec[P(f_t; \theta_0) diag \, \pi] = [diag \, \pi \otimes Id] \, vec P(f_t; \theta_0),$$
where $diag \, \pi$ is the diagonal matrix with terms $\pi_j$ on the main diagonal and $\otimes$ denotes the Kronecker product.
The common factor $f_{t}$ is strictly stationary and geometrically mixing. Thus, the same property holds for the $K^{2}$ dimensional process $ vec\left[\pi_j p_{jk} (f_t; \theta_0)\right]$. We deduce the asymptotic normality of
As expected, we obtained the following results.
(a) \ The speed of convergence of $\hat{\theta}_{n,T}$ is $1/\sqrt{T}$ instead of $1/\sqrt{nT}$, which characterizes the granularity approach. This is a consequence of the crude cross-sectional aggregation of the data in the composite approach as if the observations $y_{i,t}$ were cross-sectionally independent. However, a drawback of the granularity approach is the need to estimate the $T-1$ latent factors, which increases its computational complexity compared to the conditional MCL methods.
(b) \ The asymptotic variance is obtained from the “sandwich” formula, as it is common in a mis-specified (pseudo) maximum likelihood approach [see, Godambe1960, Hubert1967, White1982].
(c) \ The terms $p_{jk} (f_t, \theta_0)$ and $p_{jk}(f_{t-h}, \theta_0)$ depend on $f_{t}$, and $f_{t-h}$, respectively. They are correlated because of the factor dynamics (except when $\rho=0$, that is the case of an i.i.d. factor). Therefore, the covariances have to be taken into account even if we consider only a small number of values of lag $h$. It is important to notice that the sum $\sum_{h=-\infty}^{\infty} I_{0h}$ always exists due to the geometric ergodicity of the factor process.
(d) \ The asymptotic variance-covariance matrix of the MCL estimator depends on the selected set of weights $\pi$. It is out of the scope of this paper to discuss the optimal choice of weights that likely reduce the robustness of this estimator. Instead, to facilitate the comparison with the granularity-based approach, this set of weights has to be close to the structure of ratings at the different dates. As it is assumed to be time independent, a solution is to choose the set of weights close to the true unconditional structure of ratings. More specifically, we choose the fixed rating structure to be the average of the observed rating structures across time in the simulation and the empirical application. For the different states, these observed rating structures are given by the number of firms in each state at the end of each time period divided by the total number of firms transiting from non-default states. The simulation results in Appendix D suggest that this choice yields rating structures that can approximate the unknown unconditional rating structure (see Sections 5 and 6).
The above asymptotic analysis is different from the main literature on composite likelihood that usually considers either i.i.d. individuals, or finite dimensional time series [see e.g. CR2004, VRF2011].
The asymptotic variance-covariance matrix of the conditional composite maximum likelihood estimator is consistently estimated by considering appropriate sample counterparts of components $J_0, I_{0h}$.
In this section, we perform a Monte-Carlo experiment to assess the finite sample properties of estimators based on the conditional composite likelihood function and step one of the granularity approach.
The designs include $K=8$ ratings, with a higher $k$ indicating a lower capacity to repay debt, and $k=8$ denoting the absorbing state of default. These rating categories can be interpreted as AAA, AA, A, BBB, BB, B, CCC/CC and D, respectively (according to the Standard and Poor's (S & P) terminology).
Given the rating at time $t-1,$ i.e. $y_{i,t-1}=j\in\left\{ 1,\ldots ,7\right\} $, suppose that the underlying latent continuous quantitative score $y_{i,t}^{\ast}$ can be written as \[ y_{i,t}^{\ast}=\delta_{j}+\beta_{j}f_{t}+\sigma_{j}u_{i,t},\text{ }u_{i,t}\sim i.i.d.N\left( 0,1\right) , \] where the rating is determined by \[ y_{i,t}=k,k=1,\ldots,8\Longleftrightarrow c_{k}\leq y_{i,t}^{\ast} <c_{k+1},k=1,\ldots,8, \] with the thresholds $\left( c_{k}\right) $ described in Table 1 and the intercepts $\left( \delta_{j}\right) $ described in Table 2.
The thresholds and intercepts are ranked in an increasing order, and their values are chosen to get higher transition probabilities on the main diagonal and decreasing probabilities when a firm transits to other states. The treatment of the “absorbing state D” corresponding to $j=8$ is discussed later in Section 5.1.3.
The uncertainty on migrations is driven by rating-specific shocks $u_{i,t}$ and the common systematic shocks $f_{t}$. To see the effects of risk on the systematic and rating-specific components, we consider two designs for $\sigma_{j},\beta_{j}, j=1,\ldots,7$.
Design 1: $\rho=0$; The idiosyncratic and systemic components have, for each $j$, the same impact, that is: $\sigma _{j}=\beta_{j}=\frac{(1+r)^{j-1}}{\sqrt{2}}$, with $r=0.05$. Thus, when the rating is lower, the risk of downgrading is higher.
Design 2: $\rho = 0.4, 0.7$ and $0.95$; The autocorrelation parameter is taken into account and the impact of the systemic component relative to the idiosyncratic one decreases with $l$. This means that the idiosyncratic errors largely explain the junk bonds in non crisis environment. To capture this feature, we consider the ratios $\frac{\beta_{j}}{\sigma_{j}}=\frac{1}{(1+r)^{j-1}}$, with $r=0.05$, where $\beta _{j}=\frac{1}{\sqrt{2-\rho^{2}}}$, $\forall j$.
There is also a persistence of the systematic factor $f_{t}$ satisfying \[ f_{t}=\rho f_{t-1}+\sqrt{1-\rho^{2}}\eta_{t}, \eta_{t}\sim i.i.d.N(0,1), \] where the autocorrelation parameter $\rho$ measures the persistence and $f_{1}$ is drawn in th estationary distribution $N(0,1)$.
We consider the following four values for the autocorrelation parameter $\rho$: $\rho=0$, that corresponds to independent migration matrices. This is the basic assumption of the Value of the Firm model introduced in Vasicek1987; $\rho=0.4$ is used to reflect a moderate amount of autocorrelation at lag 1 of the systematic factor; $\rho=0.7$ corresponds to a high amount of autocorrelation at lag 1 of the systematic factor, while $\rho=0.95$ allows for some persistence in the systematic factor.
The state of default is an absorbing state. Therefore, if we follow a given population of corporates, all of the corporates will default at some date, and the number of still alive corporates (the so-called Population-at-Risk (PaR)) will diminish. Theoretically, the process of observed ratings is asymptotically stationary with a stationary distribution equal to a point mass on default. This difficulty is solved by assuming that newly created corporates offset the corporates entering into default, thus ensuring a PaR of constant size. This corresponds to the model with equal birth and death rates used in epidemiological studies (see e.g. HLM2014). As at the time of new firms arrival their rating are high, we replace the last row of the migration matrix at the individual level,
\[ 0, 0, 0,0,0,0,0,1, \] corresponding to a standard absorbing state, by the row of assignment of new entries at the population level, \[ 0.5,0.3,0.2,0,0,0,0,0. \]
Thus we have to distinguish individual migration matrices $P_{t}$, from the population migration that could be adjusted by taking into account the newly created firms. When the newly created firms are taken into account, the migration matrix is indexed as $P^{a}_t$.
Let us consider the design with $\rho=0.4$. For each individual $i$, we compute and compare the time series of underlying scores, ratings as well as the series of expected stability measures in the current rating. These series are denoted by $y_{i,t}^{\ast}$, $y_{i,t}$ and $s_{i,t}$, where \[ s_{i,t}=\Phi\left( \frac{c_{k+1}-\beta_{j}f_{t}+\delta}{\sigma_{j}}\right) -\Phi\left( \frac{c_{k}-\beta_{j}f_{t}+\delta}{\sigma_{j}}\right) ,\text{ with }y_{i,t-1}=j. \] These series are displayed in Figure 1 for an initial factor value of $f_{1}=0$ and initial rating of $y_{i,0}=2$, which is equivalent to AA.
The displayed trajectories correspond to three different corporate bonds. At time 0, a bond with rating 2 (AA) is issued. It is subject to downgrading after time 10 down to default in time 21. At that time a new bond with rating 1 (AAA) is issued to balance the defaulted bond. It is gradually downgraded to default at time 44. Then, a new bond is issued at time 45 and so on. In such an environment of births and deaths occurring with equal rates, each trajectory corresponds to a stochastic number of firms, rather than a single firm. This stochastic number is equal to the number of observed defaults plus one. This approach ensures the stationarity of the process and provides the rating histories of equal length $T$.
In this section we present the quasi-transition (migration) matrices with $\rho=0.4$. The time unit can be viewed as one month or one quarter, and the horizons of one and two correspond to one and two time units, respectively. Matrix $P^{a}$ is evaluated at the true parameter value from the formula in Lemma 1 and given in Table 3.
{
}
We observe a commonly reported feature of a migration matrix, i.e. the largest rates located on the main diagonal and the two adjacent diagonals with larger rates of downgrades than of upgrades. Moreover, there are significant rates of default from ratings 6 and 7, corresponding to the “junk bonds”. The last row corresponds to the new firms introduced to compensate for the defaulted corporates. Next, we determine the nondegenerate stationary distribution $\mu^{a}$, solution of
Because of the absorbing state, without the equal birth-death rates each corporate bond would default and the asymptotic stationary distribution of individual ratings would be a point mass at 8 (D). The interpretation of the stationary distribution $\mu^a$ is different and concerns the population ratings. It provides the long run rating structure of the population of corporate bonds under rebalancing. This long run structure, that does not depend on the initial rating structure, is given in Table 4. In practice, the stationary distribution provides the information on how the ratings agencies determine the thresholds of scores to define the ratings. In our experimental design, the unobserved quantitative scores are discretized to obtain close proportions of bonds across ratings. {
}
Let us now consider the quasi-migration matrix at horizon 2. Table 5 shows the matrices $P^{(2)a}$ and $(P^a)^{2}$. Matrix $P^{(2)a}$ is computed by Monte-Carlo integration with $S=50,000$ replications of the latent factor values $f_t$ at the true parameter value (see Lemma 2). The changes in $P^{(2)a}$ and $(P^{a})^{2}$ compared to $P^a$ are due to the aggregate effect of both the rating-specific and systemic shocks. Both matrices have non zero elements on the diagonals up or down by 2 from the main diagonal because of time aggregation. The matrices $P^{(2)a}$ and $(P^{a})^{2}$ are not equal. This difference is caused by the systemic risk. As expected, we observe larger diagonal elements in the matrix $P^{(2)a}$ for $j=2,...,7$ because of the persistence of the factor, which leads to more stability in the ratings.
{
}
To give some insights into the accuracy of the MCL estimation, in terms of the number of months $T$ and the factor autocorrelation parameter $\rho$, we conduct Monte-Carlo experiments. The estimation is performed with $n=1,000$ firms, including the adjustment for the newly created firms and the designs described above. The numbers of observation periods are $T=60$ ($5$ years for monthly data and $15$ years for quarterly data), $T=120$ ($10$ years for monthly data and $30$ years for quarterly data), $T=240$ ($20$ years for monthly data and $60$ years for quarterly data). In each experiment, we perform $S=500$ simulations of individual trajectories, with initial ratings $y_{i,0}$ drawn from the adjusted stationary distribution $\mu^{a}$, conditional on non-default ratings. The stationary probability computed conditional on non-default ratings provides the initial rating structure. The condition for this initial draw is satisfied by starting the simulations from another initial structure, $20$ dates before $t=1$ in our case. This structure is chosen equal to the average observed structure of ratings.\footnote{Figures (ref) - (ref) in online Appendix D.4 presents the rating structures at each time period $t$, and the chosen fixed rating structure, represented by their means over the $S$ simulated data. The figures show that the rating structures at each date $t$ are generally close to the fixed structure for most ratings. We find an improvement in the ability of the average rating structures to capture their dynamics over time as the autocorrelation in the systematic factor increases. However, as mentioned above, the choice of optimal rating structure is beyond the scope of this paper.}
The stochastic migration model depends on a large number of parameters that are $19$ identifiable parameters ($c_{3}\ldots,c_{8}$, $\delta_{1}\ldots,\delta_{7}$, $\gamma_{2}\ldots\gamma_{7}$) for the CL(1) method, where $\gamma_j=\sqrt{\sigma_j^2+\beta_j^2}$, and $27$ identifiable parameters ($c_{3}\ldots,c_{8}$, $\delta_{1}\ldots,\delta_{7}$, $\beta_{2}\ldots,\beta_{7}$, $\sigma_{1}\ldots,\sigma_{7}$ and $\rho$) for the CL(2) method, counted by taking into account the two identification conditions $c_{2}=0,\gamma_{1}=1$. For the granularity approach, further identification conditions are needed. The next sections discuss the CL(1) estimation results based on design 1, the CL(2) estimation results based based on design 2, and compare the findings to the granularity estimation results. In particular, we discuss the finite sample distributions of the estimated $c_{k+1},k=2,\ldots,7$ and $\delta_{j},j=2,\ldots,7,$ which are common for the CL(1), the CL(2) and the granularity methods. In addition, we analyze the distribution of the estimated $\beta_{j},j=2,\ldots,7$ and $\sigma_{j},j=2,\ldots,7,$ obtained from the CL(2) and the granularity approaches.
The CL(1) approach depends on the selected set of weights $(\pi_j,\; j=1,...,7)$ with weight zero on the default rating. This set of weights is set equal to the average rating structure from the simulated data conditional on non-default ratings. We illustrate the finite sample behavior of the parameter estimators by plotting the empirical pdf of the estimates (Figures (ref) - (ref) in online Appendix D.1), when $\rho=0, 0.4, 0.7$ and $0.95$ for $T=60, 120$ and $240$. Figures (ref) - (ref) present the empirical probability distributions depicted by the histograms for the estimated threshold parameters $\left(c_{3},\ldots,c_{8}\right)$. Figures (ref) - (ref) show the empirical probability distributions of the estimated intercepts $\left( \delta_{1},\ldots,\delta_{7}\right) $, while Figures (ref) - (ref) present these distributions for the unconditional variances $\left(\gamma_{2},\ldots,\gamma_{7}\right) $. In each figure, the $x$-axis shows the values of estimators, while the $y$-axis presents their frequencies. The red vertical line shows the true value of the estimated parameters and allows us to analyze to which extent the estimators are biased and discuss their distributions. When $\rho=0, 0.4$ and $0.7,$ a common feature observed in these figures is that, when $T$ varies, the distribution of the parameter estimators remains centered around their true values. The mode of estimates from the simulated data takes values close to the true value. As the sample sizes increase, the range of values taken by the estimates tends to decrease. This indicates a smaller dispersion, and, therefore, an improvement in the precision of the CL(1) estimation. In the extreme case, where $\rho=0.95$ is close to the non-stationarity in the latent factor, with $\rho=1$, the estimations are less accurate. However, the accuracy also improves with an increase in the time dimension. The results are in line with the asymptotic results on the $\sqrt{T}$-consistency of the CL(1) estimates given in Proposition 3.
We conduct further analysis based on the CL(2) method with the same set of weights. Figures (ref) - (ref) in online Appendix D.2 show the distributions of the CL(2) estimators when $\rho=0, 0.4, 0.7$ and $0.95$ for $T=60, 120$ and $240$. Like for the CL(1), Figures (ref) - (ref) provide the histograms of the estimated thresholds, and Figures (ref) - (ref) show the histograms of the estimated intercepts. In addition to those estimated parameters, Figures (ref) - (ref) display the histograms of the estimated slopes, while Figures (ref) - (ref) present these distributions for the estimated volatilities. Under the CL(2), the thresholds and intercepts are generally well estimated when $\rho=0, 0.4$ and $0.7$. The results show that the estimates are centered around the true parameters. The accuracy of the estimation improves with increasing sample sizes. We observe similar results for the factor sensitivities and volatilities.
The granularity estimation consists of two steps. First, the log-likelihood from the micro-density given by
is maximized with respect to parameter $\theta$ and factor values $f_2,...,f_T$. In the second step, the values $\hat{f}_t$ are regressed on their lagged values to get an estimator of $\rho$. We examine the first step of the granularity approach providing the estimates of $\theta$.
The estimators are presented in online Appendix D.3. Figures (ref) - (ref) provide the histograms of the estimated thresholds, and Figures (ref) - (ref) show the histograms of the estimated intercepts. Figures (ref) - (ref) show the histograms of the estimated slopes, while Figures (ref) - (ref) present these distributions for the estimated volatilities. By comparing the finite sample distributions, we observe that the CL(1) method provides slightly more accurate results than the CL(2) and the first step of the granularity approach for the parameters identifiable under the CL(1). The results obtained from the CL(2) and granularity for parameters $\beta_j, \sigma_j$ are close, although slight asymmetries arise in the histograms of some of these parameters estimated by the granularity.
In addition, the granularity approach is computationally more intensive, and the computation burden increases with the number of factor values, given that the factor is estimated for each time period at the first step of the granularity estimation.
{{{0.2cm}
} To give an idea of the computation time, the average estimation times are presented in Table (ref). The estimation is carried out on a virtual machine running in a VMware (Virtual Machine) cluster.\footnote{Intel Xeon Gold 6140, multi threaded CPU with 20 real cores, 192GB ECC RAM, and 1TB enterprise grade SAS drive disk space with RAID-6, and dual power supply.} The code is written in Matlab and the optimizations are performed by using the fminsearch. The results show that the maximum CL(1) likelihood estimation is at least $5$ up to $7$ times faster than the CL(2) estimation and granularity approach when $T=120$, respectively. Furthermore, it is at least $6$ and $18$ times faster than the CL(2) estimation and granularity approach when $T=240$. Overall, the granularity procedure is more intensive (and requires maximizations with respect to higher numbers of parameters) than the CL(1) and the CL(2) methods as the number of time periods increase, while the CL(1) is always the fastest in terms of computation.
In this section, we present the empirical results. Subsection 6.1 introduces the data. In Section 6.2, we estimate the model from the conditional composite log-likelihood at lag one proposed in Section 3.2, given that CL(1) performs well in the finite sample experiments and is computationally less intensive. We analyze the estimated parameters, transition probabilities, probabilities of defaults, and the downgrade probabilities at different horizons.
The observed transition probabilities are computed from the Compustat S&P rating database over the period 1985Q4 to 2016Q4. The average number of firms over all the $125$ quarters is $13,599$. In this section, we describe the data set and explain the necessary transformation due to missing data on non-rated companies.
We use the domestic long-term issuer quarterly credit ratings classified in eight categories: AAA, AA, A, BBB, BB, B, CCC/CC, and D, ranked from the lowest up to the highest risk, as in Section 5: $k=1,\ldots,8$, respectively. Each transition matrix summarizes all rating movements across these categories over one quarter.
As explained by Feng2008, non-rated firms arise when the relevant debt is extinguished and there is a lack of balance sheet information to determine the firm rating due to a merger or an acquisition. We follow the approach of Feng2008 to correct for non-rated firms. More precisely, we use for our analysis the transition probability conditional on being rated at the end of the quarter. We divide the frequency of migrating from any rating $j=1,\ldots,7$ to $k=1,\ldots,8$ by one minus the frequency of migrating from $j=1,\ldots,7$ to the non-rated state $j=9$. The resulting non-rated adjusted transition matrix is presented in Table (ref).
{
}
The frequencies of firms remaining in the same ratings over 1987Q2 is close to $100\%$ because the changes of ratings do not occur very often. The next highest transition probabilities are on the two diagonals below and above the main diagonal, as firms generally move away by one category from their initial rating if they are up-graded or down-graded by the rating agencies. The last column shows the probability that a rated firm defaults. The migration matrices in the sample have a similar pattern to the one documented in Table (ref).
The S$\&$P database provides the rating structure of issuers at each date (see the first column of Table (ref)). This structure changes over time due to rating migration, and also because of defaults of some issuers and the arrivals of new issuers. Therefore, there is another type of rebalancing of the population of firms. We use the time averages of these structures to define the weights in the estimation.
We report in Table (ref), the estimated parameters and their bootstrap confidence intervals.\footnote{The firms are randomly drawn with replacements. Consequently, the firms' histories are kept unchanged once a firm is drawn. Therefore, we estimate parameters from the bootstrap data and use them to find the bootstrap intervals of the estimated parameters based on $B=399$ replications.} Table (ref) contains the estimated transition matrix, and Table (ref) presents some resulting downgrade probabilities and probabilities of default.
{{{0.3cm}
}}
The estimated thresholds in Table (ref) increase with the ratings, so that firms with higher latent scores receive higher ratings. Furthermore, firms with higher ratings have higher estimated intercepts and therefore increased scores. The largest estimated value of unconditional volatility is obtained for $j=6$ and $j=7$, which corresponds to firms facing major uncertainties, currently vulnerable, or which have filed for bankruptcy protection.
In practice, the parameters of interest can be nonlinear functions of $c_{j}$, $\delta_j$ and $\gamma_j$. For instance, one might be interested in the quasi-transition probabilities, the quasi-downgrade probabilities and the quasi-probabilities of default. We used the estimates to compute the migration probabilities and illustrate the prediction of downgrade probabilities and probabilities of default as follows.
{
}
Table (ref) shows the estimates of the expected migration probabilities from any state $j=1,\ldots,7$ to $k=1,\ldots,8$, after plugging in the parameter estimates. The estimated quasi-migration matrix reproduced the aforementioned features of an observed transition matrix. Table (ref) presents the downgrade probabilities and the probabilities of default computed from the estimated transition matrix. We observe that the downgrade probability and the probability of default often increase as the horizon increases. At horizon 1, the downgrade probability is $2.72\%$ for a firm initially rated $j=1$, while it increases to $10.26\%$ for a firm initially rated $j=7$. At horizon 2, the estimated probability increases to $18.35\%$ for $j=7$ from $5.37$ for $j=1$. At horizon $1$, the probabilities of default are zero for ratings $j=1,\ldots,5$. This is consistent with the fact that the probability that a firm with a better capacity to repay its debt defaults during one quarter is negligible. The estimated probabilities increase as the horizon increases.
{
}
The stochastic factor ordered Probit model has been introduced for dynamic analysis of credit risk, as it is sufficiently flexible to account for the rating dynamics, the presence of systemic risk and the stylized fact of rating momentum. This paper proposes three maximum composite likelihood estimation methods of different complexity for this model: the conditional composite log-likelihood function at lag 1, the conditional composite log-likelihood at lag 2, and the conditional composite likelihood up to lag 2. The paper discusses the identifiability of the model parameters and establishes the asymptotic properties of these estimators when both the cross-sectional dimension $n$ and the time dimension $T$ tend to infinity. In this asymptotic setup, the MCL methods are less efficient than the two-step granularity-based approach existing in the literature. However, in practice $n$ is large, but $T$ is much smaller. The considerable advantage of the MCL approaches is that they reduce the computational burden of the granularity approach as they avoid estimating a large number of nuisance parameters. We illustrate the finite sample properties of the conditional composite log-likelihood at lag 1 and at lag 2 by conducting Monte-Carlo experiments and compare with the granularity approach. Our results indicate that the MCL methods are reliable at T=60 and are computationally less demanding than the granularity-based estimator at finite $T$. An empirical study illustrates the application of the proposed method to credit rating data.