Numerical estimates of risk factors contingent on credit ratings

Assuming a favorable or an adverse outcome for every combination of a credit class and an industry sector, a binary string, termed as a macroeconomic scenario, is considered. Given historical transition counts and a model for dependence among credit-rating migrations, a probability is assigned to each of the scenarios by maximizing a likelihood function. Applications of this distribution in financial risk analysis are suggested. Two classifications are considered: 7 non-default credit classes with 6 industry sectors and 2 non-default credit classes with 12 industry sectors. We propose a heuristic algorithm for solving the corresponding maximization problems of combinatorial complexity. Probabilities and correlations characterizing riskiness of random events involving several industry sectors and credit classes are reported.


Motivation
Intuitively, migrations towards more secure credit classes should be more frequent during economic upturns and less frequent during downturns, while migrations towards riskier states should exhibit the opposite pattern. With this observation in mind, we analyze statistically the variety of economic conditions affecting different industries and credit classes. Given the actually observed over a period of time credit-rating migrations, we identify the most likely distribution D of upturns and downturns for this period. We illustrate how the distribution can be used in financial risk analysis.
The two phases of a business cycle-a downturn or, more specifically, a contraction and an upturn or an expansion-, affecting all debtors in the economy, render their credit-rating migrations dependent. This mechanism of dependence is implemented in several models. See among others Bangia et al. (2002), McNeil and Wendin (2007), Frydman and Schuermann (2008), Stefanescu et al. (2009), Xing et al. (2012). In the same vein, Fei et al. (2012) consider three phases: an expansion, a mild recession and a severe recession. The GDP dynamics is necessary for identifying the corresponding time periods. Knowing them, the transition matrices, if discrete time is considered, or the respective infinitesimal generator matrices, in a continuous time setting, are estimated for each of the phases. In sum, the dependent credit-rating migrations are modelled relying on the conventional business cycle theory or its extension to more than two phases of a business cycle. As a consequence, the migration probabilities are not industry specific. Therefore, neither are intensities of credit-rating events considered, nor their interdependence. In this paper we propose a novel model, where the intensities as well as the interdependencies can be industry specific.
Consider a pool of debtors. Assume that there are M non-default credit classes and S industry sectors. Let every debtor be completely characterized by a couple m ∈ {1, 2, . . . , M} and s ∈ {1, 2, . . . , S}. We postulate, first, that the economic conditions affecting a debtor can be favorable or adverse and, second, that the conditions are specific for every combination of a credit class and an industry sector. Then, keeping 1 for the favorable outcome and 0 for the adverse one, a state of the economy, or a macroeconomic scenario, is encoded with a binary string having M S positions. There are 2 M S such strings. For a practically interesting choice of M and S, this can be a huge number. The conventional setting, assuming just two phases of a business cycle for the whole economy, has neither the sectoral classification nor the riskiness one. Therefore, M and S takes value 1 and there are just two strings in this case, 1 and 0.
Given a model for the probabilistic dependence among individual credit-rating migrations, a discrete time dynamic of the credit rating can be specified and the corresponding likelihood function can be written down. We formulate the model through the conditional on a macroeconomic scenario transition probabilities. Then the likelihood function depends upon a distribution D defined on the 2 M S binary strings. Given the actually observed migrations, the maximum likelihood principle allows to identify the most likely D. Should the support of this distribution contain just two sample points-the binary string formed by 1 repeated M S times and the binary string with 0 at all M S positions-this setting would reduce to the conventional model with just two phases of a business cycle. In fact, all industries and credit classes would be affected by the same economic conditions in this case. Otherwise, if there are more than two sample points in the support, analyzing the distribution D and its moments can, on the one hand, shed light on hidden, because the corresponding market outcomes are at most partly observed, dependencies among credit events in different credit classes and industry sectors and, on the other hand, imply more specific and precise, because the sectoral classification is taken into account, estimates of the risk associated with a particular portfolio. In sum, the classical two-phases regime can emerge for some input data. However, the estimates given below show that the support of D contains typically more than two sample points. Consequently, a fine-grained differentiation of the phases of a business cycle is possible. We present several new numerical estimates of financial risk factors based on this more general view of a business cycle.
Within our approach, the sectoral classification is not the only one possible. In particular, the model can be reformulated classifying debtors according to their capitalization level. Such a specification is typical for some macroeconomic risks factor models. See Aretz et al. (2010), Boons (2016), Fama and French (2017), Cooper and Maio (2019), Goncalves et al. (2020) among others. On the one hand, they employ more particular than the phases of a business cycle macroeconomic indicators-exchange rate, employment, inflation rate, GDP, etc.-, on the other hand, their predictions address the whole economy rather than a credit class. In other words, conceptually such models are richer, but, ceteris paribus, our analysis is more fine-grained.
For testing the approach, we use the same S&P's dataset as Boreiko et al. (2017). See https://doi.org/10.1371/journal.pone.0175911.s001. Two couples of M and S are considered: M = 7 with S = 6 and M = 2 with S = 12. Since 2 42 and, respectively, 2 24 macroeconomic scenarios are involved, the number of unknown parameters greatly exceeds the number of available observations. A heuristic algorithm is suggested for solving the corresponding likelihood maximization problems of combinatorial complexity.
A numerical technique for evaluating the risk factors contingent on the dependence among credit-rating migrations is the main contribution of the paper. The key element is an estimate of the hidden distribution that shapes the migrations.
In the next session we specify the probabilistic dependence among individual migrations, the corresponding likelihood function and the non-linear programming problem for estimating the unknown parameters. A heuristic algorithm for solving this problem with a combinatorial number of unknowns is described in Sect. 3. Several new risk estimates are suggested in Sect. 4. A detailed characterization of the input data and the heuristic solutions is presented in Sect. 5. Section 6 contains the main results of the paper-several risk estimates based on the distribution D. Section 7 concludes. Auxiliary technicalities are given in Appendix 1. Appendix 2 contains some statistical characteristics of the input data.

Dependence model and optimization problem
Modeling credit-rating migrations as trajectories of time-homogeneous Markov chains has a long tradition in financial risk analysis. It is implemented in CreditMetrics, a toolkit for understanding and managing credit risk. See Gupton et al. (1997). In the model considered next, a time-homogeneous discrete-time Markov chain is as well one 123 of the two components governing credit-rating migrations. This component is termed as idiosyncratic because, should it be the only driver of the credit rating migrations, they would not depend on each other and the resulting dynamic would be Markovian.
Let us formulate precisely the assumptions regarding the credit-rating dynamics. Consider a portfolio involving debtors classified into S ≥ 1 industry sectors. There are M + 1 levels of creditworthiness. Non-default credit classes are numbered in descending order of creditworthiness: the most secure assets are labelled by 1, while the next to default credit class is indexed by M. Defaulted debtors receive the index M + 1. They will never return to business.
Denote by {0, 1} M S the set of binary strings (or vectors) V with M S positions (or coordinates). They encode all possible macroeconomic scenarios in the economy. We need a rule for assigning coordinates of a binary vector V to industries and credit classes. To this end, let the coordinate V M(s−1)+m characterize the economic conditions affecting the credit class m of the industry s. That is, the industries occupy blocks of M coordinates each. The blocks are numbered in ascending order of s. Within a block, the credit classes are listed in ascending order of m.
Let X (t) be the credit rating of a debtor at time t = 1, 2, . . .. The rating randomly changes in time, becoming X (t + 1) at t + 1, while the assignment to an industry remains the same. We postulate that X (t + 1) is a weighted sum of an idiosyncratic ξ t and a common η t component: First, we characterize distributions of the random variables δ t , ξ t and η t . Let the debtor belong to an industry s and let X (t) = m for some m ≤ M. Then: • δ t stands for a Bernoulli random variable whose probability of success equals Q m,s . The M × S matrix Q with entries Q m,s is not known. It has to be estimated. • Migrations in industry s that cannot be attributed, at least directly, to market mechanisms are characterized by an M × (M + 1) migration matrix P (s) . It is known. Its entry P (s) i, j equals the probability that a debtor moves from a credit class i to a credit class j in one time instant. The distribution of ξ t is: • The conditional on V M(s−1)+m distribution of η t reads: Here, P Then the conditional on a macroeconomic scenario V distribution of X (t + 1) can be written down in the following way: Note that, defining the conditional migration probabilities, our model does not require identification of the time periods corresponding to the phases of a business cycle. Consequently, no GDP dynamics is necessary to estimate the probabilities. There are several possibilities to relate the trend of credit rating migrations to the phases of a business cycle. In the simplest case, it is required that the probabilities of migrating towards more secure credit classes are higher for economic upturns than for downturns. Incorporating such linear inequality constrains into a likelihood maximization problem, two migration matrices, one for upturns and the other for downturns, are estimated as well as the probability of an upturn. For survey analysis, this approach was used in Hölzl et al. (2019). Second, let us characterize the dependence among the random variables involved in (1). For a fixed time instant t: • for every debtor, the random variables δ t , ξ t and η t are independent; • the random variables δ t , ξ t are independent across debtors, while the random variables η t are conditionally, given a macroeconomic scenario, independent across debtors; • δ t , ξ t and η t do not depend upon their realizations in the past.
Considering a period of observation from t = 1 to t = T , these assumptions imply the following likelihood function: In the formula for L we ignore the multiplier 123 that does not contain the unknown parameters. Transition counts I t (s, m 1 , m 2 ) for the period of observation is the only input required. By I t (s, m 1 , m 2 ) we denote the number of debtors in the industry s that migrated from the credit class m 1 to the credit class m 2 in the period t. (The industryspecific Markovian migration matrices P (s) are estimated from these counts as time averages.) All the unknowns have to satisfy the box constraints: D(V ) ∈ [0, 1], Q m,s ∈ [0, 1] and m,s ∈ [0, 1]. Since D is a probability distribution, the following equality has to be satisfied: The value Q m,s determines the impact of the idiosyncratic component in (1). In particular, Q m,s = 1 implies that migrations of debtors belonging to the credit class m and the industry s are independent on migrations of other debtors. The existing models assume that m,s = 1 for all possible m and s. See Kaniovski and Pflug (2007), Wozabal and Hochreiter (2012) or Boreiko et al. (2017). Suggesting more general formulas for the conditional probabilities and estimating m,s , we attempt to test empirically this assumption. All other things being equal, a larger m,s implies smaller probabilities P (s) m, j (1), j < m, and larger probabilities P (s) m, j (0), j > m. Consequently, the known models may overestimate the default rates P (s) m,M+1 (0). The conditional distribution of X n (t + 1) defined by (2) typically deviates from the m-th row of the corresponding P (s) . To guarantee that the unconditional distribution of X n (t + 1) in (1) equals P Here,P Having a distribution D, probabilities of fine-grained macroeconomic events associated with the phases of a business cycle can be estimated. For example, probabilities of adverse for a combination of credit-classes and industries events. Since a time period in the past is considered, frequencies is a more precise denomination for such values. Numerical characteristics other than probabilities can be evaluated as well.
In particular, correlations between the economic factors affecting different industries and credit classes.
In what concerns credit risk, our approach, accounting for the sectoral dimension, implies a more specific and, consequently, precise analysis of credit events. In fact, the migration matrices are industry-specific and, in the course of a business cycle, they are adjusted according to the industry-specific economic conditions.
Assessing how likely is a contagion in the financial sector, probabilities of adverse events affecting several credit classes can be used. Unlike the majority of models addressing systemic risk, our estimates do not require any information regarding interconnections in the banking sector. See Glasserman and Young (2016) for a comprehensive review of systemic risk models.
Let us show that, under natural assumptions regarding the migration matrices P (s) , the feasible set defined by the linear relations (3) and (4) is not empty. Since the probabilities Q m,s are not involved in these relations, it is enough to indicate a distribution D and probabilities m,s satisfying (3) and (4).
Let p = max s,mP  m,m > 1 2 for all combinations of m and s in the migration matrices quoted in Appendix 2, that the assumptions of Proposition hold true for our input data.
Since quarterly P (s) m,m typically do not fall below, all other things being equal, their annual counterparts, verification of assumptions of Proposition is easier dealing with quarterly transition counts.
Even if concavity of ln L cannot be established, its derivatives are known and, apart from the box constraints, there are just linear equality constraints. If 2 M S is not too large, a standard solver can be used for finding a solution. Restarting the algorithm for different initial points, the solution can be improved.
This course of action was implemented by Boreiko et al. (2017) for M = 2 and S = 6. The Interior Point (IP) method was used. (Postulating that all m,s are equal to 1, the model in Boreiko et al. (2017) is slightly simpler.) On the one hand, this approach looks like a classical optimization technique. In fact, the classification in Gilli and Schumann (2012) requires for such a technique "at least well-behaved objective functions" so that a kind of gradient descent can be implemented. See p. 130. On the 123 other hand, a heuristic-restarting the gradient descent from different initial points-is used because concavity of the likelihood function cannot be guaranteed. (As a matter of fact, all of the restarts implied the same solution.) In sum, a combination of a classical technique and a heuristic is applied even in this simplest case.
A more detailed analysis of credit-rating migrations considers two non-default credit classes as well, but the number of industries is doubled. Therefore, the total number of unknowns becomes 24 + 2 12 . This level of computational complexity requires special algorithms. We suggest and test one of such methods.
Maximizing ln L, there are two challenges. First, the number 2 M S of unknown probabilities D(V ) is typically combinatorial for a practically interesting choice of M and S. Second, the estimated distribution D can be nested in too many sample points. In such a case, it is difficult to analyze the corresponding macroeconomic scenarios.
Confronting with a large M S, the sample space {0, 1} M S can be split into smaller parts so that the likelihood maximization problem could be solved for each of them with a standard algorithm. Since there can be too many such subsets, instead of analyzing all of them, a random sample can be used. Passing from one of the subsets to the next one, we need a rule for identifying the sample points that are retained. In sum, there are two optimization processes: a continuous space search for inputs Q, and D maximizing the likelihood function defined on a subset of {0, 1} M S and a discrete space search for a subset of {0, 1} M S with a greater maximum likelihood value. According to the classification given in Gilli and Winker (2009), this is a combination of a local search method and a constructive method. See p. 87. A heuristic algorithm based on these principles is described next.

Heuristics
The heuristic search for a better solution relies on several concepts. Let us explain and motivate them.
We begin with characterizing the subsets that can be used in a partition of the sample space. We call V ⊆ {0, 1} M S a suitable set, if there exists a probability distribution D V such that the support of D V belongs to V and the relations hold true for some i,s ∈ [0, 1]. Note that every extension V of a suitable set V, that is a subset of {0, 1} M S containing V, is suitable as well. In fact, let Keeping the values i,s unchanged, relations (6) for V will hold true with this distribution D V .
Next, let us describe a continuous space optimization problem for identifying optimal inputs Q, and D corresponding to a given suitable set. Consider a suitable set V, the likelihood function That is, our analysis is restricted now to the distributions D that are nested in V.
Ceteris paribus, the values of L and L V coincide for such a D. Since V is a suitable set, the feasible set defined by the linear equations (7) and (8) is not empty. In fact, the distribution D V involved in the definition of V, see relations (6), satisfies these constraints.
If the cardinality of V is low enough, maximizing ln L V subject to the linear constraints (7) and (8), the corresponding optimal distribution D * V can be estimated. A standard solver can be used, say the IP method. Since all unknowns are probabilities, the box constraints for Q i,s , i,s and D(V ) are necessary as well. Restarting the algorithm from different initial points, the solution can be improved.
If D * V (V ) is sufficiently small, the macroeconomic scenario encoded by V is not (statistically) significant and, therefore, it can be ignored passing from V to its extension V . A heuristic for retaining significant elementary outcomes is described next.
For a practical numerical search algorithm in {0, 1} M S , consider a threshold ∈ (0, 1) and the -support V of D * V : Note that V ∈V D * V (V ) ↑ 1 as decreases and D * V satisfies constraints (7) and (8). Therefore, V is a suitable set for all sufficiently small . In fact, the linear constraints (6) and (8), where V and D are substituted with V and D * V , can be satisfied with any precision.
If decreases, the number of binary strings in V can increase. In fact, if the strings were equally probable, the total number of them would grow as −1 . A non-negligible D * V (V\V ) can be another effect of a sparse D * V . In particular, D * V (V \ V ) > . Keeping the threshold constant, the following heuristic can be used to obtain a more concentrated than D * V optimal distribution. As a measure of concentration, we require that the total probability of the significant macroeconomic scenarios exceeds 1 − .
To this end, consider We interpret R ≥ 0 as a penalty parameter. It determines the relative importance of the two components inL V . Let us maximizeL V subject to the linear constraints (7) and (8). For a sufficiently large R, we expect that the optimal distributionD V obtained in this way should be more concentrated than D * V . (The latter corresponds to R = 0.) As a heuristic argument supporting this intuition, observe that the function of N variables N k=1 x 2 k attains its maximum value 1 under the constraints if one of the addends equals 1. (Its minimum value obtains if all of them are equal.) Adding such a term to ln L, we expect, given the box constraints and (8), to get a more concentrated than D * V optimal distributionD V . Since the weight of the second term inL V increases in R, we expect the same forD V (V ). Therefore,D V (V ) > 1 − for all sufficiently large R. ByV we denote the -support ofD V . Now the rule for identifying the significant macroeconomic scenarios can be formulated more precisely. Fix a maximum mismatch in the linear constraints (7) and (8). Denote this value by . The setV such thatD V (V ) > 1 − is a suitable set. In fact, substitutingV andD V instead of V and D V into the linear constraints (7) and (8), the maximum error will not exceed . Therefore, for this numerical precision,V is a suitable set.
The structure of our discrete space search is as follows: 1. At the beginning, a suitable set V with a low cardinality has to be identified and a threshold value has to be chosen. 2. Solve the optimization problem for specifying V . Decrease if V is not suitable.
Keep as large as possible. If V contains too many elements or/and D V (V \ V ) ≥ , consider the correspondingL V and identifyV . Increasing R, achieve thatD V (V ) > 1 − . Keep R as small as possible. 3. Better solutions (in terms of the likelihood value) can be found by extendingV to a larger set. Since any set containing a suitable set is also suitable, the extension preserves suitability. Let V =V ∪Ṽ, whereṼ is a subset of {0, 1} M×S \V. This is a suitable set, an extension ofV .Ṽ should not be too large so that the optimization problem with V instead of V could be solved. 4. We repeat steps 1.-3. with V replaced by V and keep on doing extensions.
There are two clarifications. First, we can always argue about an extension ofV . In fact, V =V for R = 0. Second, since V ⊃V , the extension cannot lead to a smaller maximum value of L.
If the assumptions of Proposition hold true, every set V ⊆ {0, 1} M S containing the string with 1 at all positions and the string with 0 at all positions is suitable. Consequently, the discrete space search can always be initiated. We suggest a particular extension process exploiting the interpretation of a binary string as a macroeconomic scenario.
Let us describe this particular scheme for generating extensions. Initially, we assume that every industry is affected by the same economic conditions. The corresponding binary M S-string is termed as a block-structure. In a block-structure, all binary Msubstrings allocated to different industries coincide. There are 2 M block-structures. Containing the macroeconomic scenario favorable for all industries and credit classes as well as the macroeconomic scenario adverse for them all, the set V of all blockstructures is suitable.
Solve the optimization problem for identifying D * V . This is not a hard task. In fact, just 2M S + 2 M unknowns are involved. Let V k , k = 1, 2, . . . , K , be the blockstructures whose probabilities D * V (V k ) exceed . It is convenient to number them in descending order of the probabilities. In order to reduce K , the modified objective functionL V can be used instead of L V . Denote by v k the binary M-substring such that V k is formed by S copies of v k . We call v k a block.
A scheme for generating new macroeconomic scenarios as mutations of the blockstructures V k , k = 1, 2, . . . , K , is motivated by elementary facts from genetics. In particular, substituting one block v k in V k by a block v i , i = k, we get a mutant with a single mutation of the block-structure V k . Conceptually, this is a macroeconomic scenario, where all industries, except for one, are affected by the economic conditions encoded with v k , while the remaining industry is affected by the economic conditions summarized in v i . Trying all possible V k , k = 1, . . . , K , industries s = 1, . . . , S or, equivalently, positions of the mutated block, and its type v i , i = 1, . . . , K , i = k, all mutants V 1 with a single mutation will be obtained. Similarly V n , can be defined for n > 1. That is, V n contains all possible mutants with n mutations. The mutations can coincide: defining V n for n > 1, a block v i can appear more than one time in a mutant stemming from the block-structure V k . Then V n ∩ V n = ∅, if n, n < S/2, n = n . Defining the sets V n for n ≥ S/2, some of the mutants, originating from different block-structures, can be listed several times. Such repetitions have to be avoided. For example, if S is an even number, a mutant containing S/2 blocks v k and S/2 blocks v i can be regarded as originating from the block-structure V k as well as originating from the block-structure V i . Such mutants have to be listed in V S/2 just one time.
There are K S possibilities for allocating K building blocks among S positions. If K < 2 M , this number is smaller than 2 M S . Consequently, restricting the analysis to all possible combinations of blocks, on the one hand, reduces the scope of macroeconomic scenarios, but, on the other hand, the corresponding solution should not be necessarily optimal even if all of the combinations were considered. A smaller K implies a stronger reduction of the search space.
Instead of dealing with all binary vectors from V n , a random sample can be used. Generating a binary vector from V n , for the first mutation there are S equally probable positions, for the second-the remaining S − 1, etc. The simplest way for choosing mutations, is to sample them with equal probabilities 1 K −1 . The resulting distribution of mutants is referred to as uniform.
The S&P's dataset used for estimating unknown parameters contains 103723 transition counts. See https://doi.org/10.1371/journal.pone.0175911.s001. We consider two classifications: M = 7 with S = 6 and M = 2 with S = 12. The corresponding numbers of unknowns are 84+2 42 and 48+2 24 . They are huge and, more importantly, these numbers greatly exceed 103723-the number of available observations. Common sense suggests that no classical approach can be used. Looking for a way out, we observe that there are 2 42 (2 24 ) probabilities D(V ) among the unknowns. Dealing with a real life problem, unlikely outcomes can be ignored. A lower bound for probability of a significant outcome implies that at most −1 such outcomes have to be considered. In particular, if = 10 −3 then at most 1000 macroeconomic scenarios are statistically significant. This number is approximately 100 times smaller than 103723. Therefore, taking into account that the upper bound of 1000 corresponds to a practically impossible situation when all of the significant outcomes are equally probable, the task of estimating D does not look hopeless. In sum, a classical solution is not possible given the input, instead a conceptually plausible approximation for D and, consequently, for the hidden dependence structure among industries and credit classes can be suggested. Motivating this approach, we rely on the principles formulated in Gilli and Winker (2009): "Often the term 'heuristic' is linked to algorithms mimicking some behavior found in nature, ... a heuristic should be able to provide high-quality (stochastic) approximation to the global optimum at least when the amount of computational resources spent on a single run of the algorithm ... is increased." See p. 83.
Searching for a plausible approximation, we first hypothesize that all industries are affected by the same economic conditions. This assumption corresponds to the standard business cycle theory. A new element is a fine-grained classification of creditworthiness. Under such a restriction on the class of admissible macroeconomic scenarios, we identify the significant outcomes, referred to as block-structures. Since no penalty is necessary and the restarts of the IP method from different initial points result in the same likelihood value, all criteria of Gilli and Schumann (2012) for a classical optimization technique are satisfied. See p. 130. Therefore, this is a classical solution. Next we consider a richer set of admissible macroeconomic scenariosblock structures and their mutants having just one mutation. It is a further extension of the standard business cycle theory because industries can be affected by different economic conditions. This extension is the simplest possible since at most one variation per macroeconomic scenario is allowed. For all choices of M and S considered in this paper, the number of unknowns allows, given our computational resources, to use the IP method. Since a penalty was necessary, an approximation to the solution was obtained for this extension. Next we turn to macroeconomic scenarios where at most two of the industries can be affected by different economic conditions. Conceptually, this is a further generalization of the standard business cycle theory. We list all mutants with exactly two mutations. Given our computational resources, the corresponding number of unknowns for the couple M = 7 and S = 6 is too high. Thus, we are forced to split this set into equally large subsets. First, we consider one of the subsets and the significant binary strings with at most one mutation. The corresponding continuous space maximization problem is solved with the IP method. The significant binary strings identified at this step are considered together with the next subset, etc. Having tested all of the subsets, we obtain a heuristic approximation to the most likely, given the input, distribution on the binary strings with at most two mutations of the block-structures. This is just an approximation. In fact, non-zero penalties were used in several extensions and, more importantly, we cannot guarantee that the greatest likelihood value achieved by the sequential analysis of the subsets coincides with the maximum value attained on the union of them, should the corresponding maximization be done. If M = 2 and S = 12, all binary strings with two mutations are used for a single extension. Further extensions are done in the same way with randomly generated mutants having 3, 4, . . . , S − 1 mutations. No binary string may be considered more than one time. All necessary details are given in Appendix 1.
As a standard measure of performance for the heuristic, we present in Appendix 1 the percentage of increase of the likelihood value. Also, since for M = 2 and S = 6 a classical solution exists, we compare it in Appendix 1 with the heuristic solution. This is an instructive example illustrating the particularities of both the likelihood maximization problem and the heuristic algorithm.
The discrete space search described above mimics biological evolution: first, only macroeconomic scenarios more probable (or more frequent, because we refer to a period in the past) than a threshold are retained to the next round of selection (where they compete with the new scenarios of the extension) and, second, more complex macroeconomic scenarios obtain as mutations of the original forms-the block-structures. The minimum viable population argument is a particular justification for the selection rule. Business cycles is another natural point of reference: we begin with the standard setting and attempt to arrive at a more fine-grained view. Some of the reported solutions are exact, but typically an approximation is found. The quality of a solution increases as the computation efforts do. In particular, widening the scope of admissible macroeconomic scenarios, a more realistic view of a business cycle can be obtained. However, it requires solving more complicated optimization problems. The random search in the space of mutants with n ≥ 3 mutations exhibits, according to our experience, path dependence. Moreover, the shape of L depends upon the transition counts available. For this reason, any estimate of the empirical distribution of the maximum likelihood value would be input-specific as well. Therefore, such an empirical convergence analysis, even if it is recommended in the literature on heuristic methods in econometrics, see among others Gilli and Winker (2009) or Gilli and Schumann (2012), cannot be used. A bootstrap procedure can be suggested instead. Then, for a set of estimated parameters, new credit-rating migrations are generated according to formula (1). Applying the algorithm, this input is transformed by into a new set of estimates. If they match reasonably well the original values, we conclude that the algorithm works correctly. For a simple example of such a convergence analysis see Boreiko et al. (2016).
Given D, Q and , some risk estimates are suggested in the next section.

4 Methodology
The following quantitative characteristics of risk can be evaluated: • Probability of a favorable, if I m i s i = 1, (an adverse, if I m i s i = 0,) for the industry s i and the credit class m i macroeconomic scenario, i = 1, 2, . . . , L. Both indexes can repeat. That is, this probability can concern more than one credit class from the same industry or/and several industry sectors having the same creditworthiness.
• Correlation between the indicators of the following macroeconomic outcomes: if I m 1 s 1 = 1 (I m 1 s 1 = 0), one of them is favorable (adverse) for the credit class i 1 of the industry s 1 , the other one is favorable (adverse) for the credit class i 2 of the industry s 2 if I m 2 s 2 = 1 (I m 2 s 2 = 0).

Inputs and heuristic solutions
If M = 7, the S&P's credit classes A A A, A A, A, B B B, B B, B and C are numbered by 1, 2, . . . , 7. The following 6 industry sectors are considered: 1-agriculture, mining and construction; 2-manufacturing; 3-transportation, technology and utility; 4trade; 5-finance; 6-services. Wozabal and Hochreiter (2012) as well as Boreiko et al. (2017) dealt with the same choice of industries.
In the case of two non-default credit classes, index 1(2) is assigned to the investmentgrade (non-investment-grade) debtors. The investment-grade debtors occupy the S&P's ratings from A A A to B B B. The non-investment-grade debtors are those whose creditworthiness is lower: B B, B or C. There are 12 industry sectors: 1industry sectors are consideredaero, auto, capital goods, metal; 2-consumer, service; 3-energy, natural resources; 4-financial institutions; 5-forest and building products, homebuilders; 6-health care, chemicals; 7-high technology, computers, office equipment; 8-insurance, real estate; 9-leisure time, media; 10-telecommunications; 11-transportation; 12-utilities. Nagpal and Bahar (2001) estimated default correlations for these 12 industries and two non-default credit classes.

Input
We use annual transition counts I t (s, m 1 , m 2 ) covering the period 1991-2015: t = 1 corresponds to 1991 and T = 25 corresponds to 2015. These transition counts are available at https://doi.org/10.1371/journal.pone.0175911.s001. Distribution of the counts among credit classes and industries is characterized in Appendix 2.

Solution
All continuous space optimization problems were solved with the IP method. The significance threshold equals 10 −3 everywhere.
Consider first the couple, where M = 7 and S = 6. There are K = 14 blockstructures whose probabilities exceed 10 −3 . Total probability assigned to the remaining block-structures does not exceed 0.0506%. Numbering the significant V k in descending order of their probabilities, the corresponding blocks v k are listed in Table 1. Searching for the heuristic solution, mutants with up to 5 mutations were considered. (The most important technicalities regarding this search are given in Appendix 1.) Table 2 contains the 22 significant macroeconomic scenarios and their probabilities. Interestingly enough, the block v 14 does not appear in these scenarios. Total probability of the remaining sample points falls below 0.0223%. In particular, probability of the significant macroeconomic scenario numbered by n = 1 is 0.2661. It can be thought of as a mutant with four mutations of the block-structure V 1 . In fact, the block v 1 is assigned to industries 1 and 5, while the remaining industries received all different blocks: v 13 , v 8 , v 11 and v 4 . Since this macroeconomic scenario allocates the block v 8 to the industry 3, economic conditions are favorable for the transportation, technology and utility firms rated at A A A, A A, B B B and B B, while the debtors in this industry sector rated at A, B and C are affected by adverse economic conditions. In fact, the row of Table 1 allocated to v 8 contains 1 at the positions 1, 2, 4 and 5, while 0 is assigned to the remaining cells. The significant scenarios listed in Table 2 contain up to five mutations. (In this case all blocks in the corresponding row are different.) The minimum number of mutations in the macroeconomic scenarios quoted in Table 2 123 is three. Therefore, the simpler outcomes containing zero, one or two mutations were wiped out in the course or calculations. The content of Tables 1 and 2 suggests that the significant macroeconomic scenarios do not exhibit the pattern assumed within the classical business cycle theory.
A referee, whose insight was crucial for streamlining this paper, observed that some of the blocks in Table 1 encode scenarios that look counterintuitive: being favorable for lower rated debtors they are adverse for higher rated ones. For example, the scenario encoded by v 1 is adverse for debtors rated at A, but it is favorable for those rated at B B B. He suggested to consider only "monotone" blocks. In particular, dealing with 7 non-default credit classes, there are 8 such blocks: 1111111, 1111110, 1111100, . . ., 1000000, 0000000. This proposal is promising because, on the one hand, the counterintuitive outcomes are excluded and, on the other hand, the total number of unknowns can be substantially reduced. In fact, dealing with M non-default credit classes and S industries, the total number of sample points will be (M + 1) S instead of 2 M S . We are going to test the corresponding model numerically.
The following matrices Q and were estimated: If M = 2 and S = 12, all four block-structures are significant macroeconomic scenarios. In Table 3 blocks are listed in descending order of probabilities assigned to the corresponding block-structures. The main facts regarding the heuristic search in this case are summarized in Appendix 1. Table 4 contains the significant macroeconomic scenarios and their probabilities. Total probability assigned to the remaining sample points does not exceed 0.0506%. The matrices Q and read: In the considered so far sectoral models of dependent credit rating migrations, all entries of are assumed to be equal to 1. See Boreiko et al. (2017), for example. Our Having Q, and D, transition paths mimicking the actually observed historical migrations can be simulated. Using the Monte-Carlo method, losses generated by a portfolio can be estimated. Rather than moving in this direction, we concentrate in the next section on the risk characteristics given by formulas involving the estimates of Q, and D.

Estimates
If Q m,s = 1, macroeconomic factors do not affect migrations in the credit class m of the industry s, so such combinations are not considered further. As a curious fact regarding the combination of M = 7 with S = 6, note that Q 1,5 = 1 can be a quantitative argument supporting the "too big to fail" theory in the financial sector. In particular, an evidence of implicit too-big-to-fail bail-out guarantee policies of the regulatory authorities. If M = 2 and S = 12, Q 1,4 = 1 can be interpreted in the same way. In sum, the financial sector debtors belonging to the most secure credit classes seem to enjoy a special treatment.
Let us analyze first the results for the combination of M = 7 with S = 6. Given in Table 5 percentages regarding the financial sector imply that the credit classes AA, BBB, B and C exhibit under adverse conditions much stronger increase of default rate than A and BB. (Since Q 1,5 = 1, the variation equals zero for the credit class AAA. See the formula for v   resort for investment grade debtors. Analyzing the respective entries of Q, provides with a quantitative argument supporting such an explanation. In fact, Q 2,5 = 0.5809 and Q 3,5 = 0.6219 imply that credit rating migrations of the debtors rated at A and at AA are strongly affected by market forces, while Q 4,5 = 0.9482, the second largest entry of Q estimated for the financial sector, causes almost idiosyncratic migrations of the debtors rated at BBB. Of course, a list of the downgraded investment grade debtors, rather than anonymized data, would be necessary for a more comprehensive explanation. Arguing about how likely is a contagion, as a factor of systemic risk, in the European financial networks, Glasserman and Young (2015) "consider the possibility that the failure of a bank causes the next two largest banks to default". See p. 396. Depending upon the country, the largest bank can be an investment grade debtor, like Deutsche Bank in Germany, a well as a non-investment grade debtor, like Alpha Bank in Greece. The frequencies quoted in Table 8 can be an input for such an analysis. We consider three out of the four credit classes with the highest increase of default probabilities 123 Table 8 Frequencies π m 1 ,m 2 ,m 3 5,5,5 (0, 0, 0) m 1 = 2, m 2 = 6, m 3 = 6 m 1 = 2, m 2 = 6, m 3 = 7 m 1 = 6, m 2 = 6, m 3 = 7 m 1 = 6, m 2 = 7, m 3 = 7 0.2839 0.2839 0.3743 0.3743   Tables 7 and 8, we conclude that the frequency of an adverse outcome for a triple of credit classes coincides with the frequency of an adverse outcome for the most creditworthy entity involved in the triple. Given that the creditworthiness is typically positively related to the size, the pattern exhibited by the frequencies in Table 8 supports indirectly the argument of Glasserman and Young. (For a less approximative analysis, a dataset, where the counts regarding banks are separated from the counts characterizing the remaining financial institutions, is necessary.) Table 9 characterizes the increase of default rate under adverse conditions in the (non-investment grade) credit class B for six industries. For a comparison, Table 10 quotes analogous values for the non-investment grade debtors classified into twelve industry sectors.
Tables 11 and 12 contain correlations between indicators of favorable outcomes, a measure of dependence between the corresponding couples of a credit class and an industry sector. These dependencies are not directly observable, but they affect occurrence of credit events involving these couples. The values quoted in Table 11 concern all seven non-default credit classes of two industry sectors: correlations regarding agriculture, mining and construction are quoted below the main diagonal, while those for transportation, technology and utility age given above the main diagonal. Every diagonal element equals one as the correlation coefficient of an indicator with itself. Table 12 contains correlations characterizing couples of an investment grade credit debtor and a non-investment grade debtor classified into twelve industry sectors. Only combinations of industry sectors corresponding to the entries of Q that fall below 1 are included. (Remember, if Q m,s = 1, macroeconomic factors do not affect migrations in the credit class m and industry s.) In a portfolio, strongly positively correlated assets   can provoke cascades of defaults. The simulations reported in Kaniovski and Pflug (2007) illustrate this possibility. To the contrary, combining strongly negatively correlated assets can imply smaller losses. To this end, the cells with correlations given in italic mark the combinations that should be avoided, while the cells with correlations given in bold correspond to the couples of credit classes that can mitigate losses in the corresponding industries. The cutoff levels are ±0.5.

Conclusion
A numerical technique is suggested for estimating hidden and observable risk characteristics. It exploits a probability distribution on macroeconomic scenarios. A scenario characterizes the conditions affecting every combination of a credit class and an industry sector. The conditions can be favorable or adverse. Given historical migrations, the maximum likelihood principle is used to estimate the distribution. As a result, a fine-grained extension of the standard business cycle theory emerges. Since throughthe-cycle ratings were used, this fine-grained structure is rather unexpected. See Kiff et al. (2013) for an analysis of the through-the-cycle rating approach. We guess that, applying our technique to point-in-time ratings, even a more interesting pattern can emerge. Dealing with such an input, constraints (4) can be too restrictive. Transition counts is the only required input. For practically interesting combinations of the number of non-default credit classes M and the number of industry sectors S, the number of unknown parameters is combinatorial and it greatly exceeds the number of observations. Therefore, the corresponding non-linear estimation problem cannot be solved with a classical method. Instead, a heuristic algorithm was suggested. It entails a local search method and a constructive method. Performing the local search, a set of macroeconomic scenarios is fixed. This is a continuous space optimization problem. All results presented in this paper were obtained with the IP method as a continuous space solver. The constructive search in the space of macroeconomic scenarios resembles a genetic algorithm (GA). Unlike in a classical GA, recombination is not used and an allele entails more than one digit in our case. Instead of dealing with all mutants, a random sample can be used in some instances.
Two combinations of M and S are considered: M = 7 with S = 6 and M = 2 with S = 12. The corresponding numbers of unknowns are 84 + 2 42 and 48 + 2 24 .
In Hölzl et al. (2019) a similar numerical technique was applied for analyzing survey data.
Funding Open access funding provided by Libera Università di Bolzano within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix 1
The threshold for identifying the significant macroeconomic scenarios is 0.001. Solving a continuous space optimization problem, up to five different initial points were tried in several cases. In all instances the restarts did not affect the solution.
First, the set V 0 0.001 of significant block-structures is extended with V 1 . If M = 7 and S = 6, there are 14 binary strings in V 0 0.001 and 6552 binary strings in V 1 . The respective numbers equal 4 and 1728 for the combination of M = 2 with S = 12. The continuous space optimization problems for this extension, containing 6650 and, respectively, 1780 unknowns, were solved with the IP method. Since R = 10 was used for both combinations of M and S, these are approximations rather than exact solutions. Denote byV 1 0.001 the corresponding set of all statistically significant macroeconomic scenarios.
Second,V 1 0.001 is extended using V 2 . For M = 2 and S = 12 there are 1584 strings in V 2 . The corresponding continuous space optimization problem was solved with the IP method. Since R = 20 was used, an approximation was obtained. If M = 7 and S = 6, there are 35490 strings in V 2 . For our computational resources, the corresponding continuous space optimization problem is too hard. Therefore, we had to split V 2 into 10 subsets, each containing 3549 elements. They were used one after another for extensions. Two values of R were applied: 20 and 50. At the end, an approximation to the solution for the sample space limited to binary strings with at most two mutations was obtained. Third, the set of significant macroeconomic scenariosV 2 0.001 with at most two mutations is extended using a random sample from V 3 . In the sample, binary strings may not repeat. If M = 7 and S = 6, the sample contains 35000 strings. This set comprise 10 subsets, each formed by uniformly, that is with the probabilities specified in Sect. 3, sampled 250 mutants of a block-structure. The subsets are used one after another for extensions. The penalties term with R = 20 and R = 50 was applied in some of the extensions. At the end, an approximation to the solution for the sample space limited to binary strings with at most three mutations was obtained. If M = 2 and S = 12, sampling uniformly 350 mutants for each block-structure, 1400 binary strings were generated. They were used for a single extension. Since R = 20 was applied, an approximation to the solution for the sample space limited to binary strings with at most three mutations was obtained.
Thereafter extensions employed random samples from V n , n = 4, 5, . . . , S − 1. Always 250 mutants were sampled uniformly for each block-structure and these subsets were used sequentially for extensions if M = 7 and S = 6, while, sampling uniformly 350 mutants per block-structure, a single bulk of 1400 binary strings was used if M = 2 and S = 12. With such numbers of mutants, the total number of binary strings used for extensions remains approximately equal to the cardinality of V 2 . In several cases, the value of R was 10, 20 or 50 dealing with S = 6 and M = 7, while in the case of S = 12 and M = 2 the three values, 10, 20 and 30, of R were applied. Table 13 characterizes the improvement of the maximal likelihood value L H achieved by the heuristic algorithm. As the benchmark the maximal value L 0 attained on the block-structures is used. A smaller relative increase of the likelihood value for a larger M S is consistent with a smaller ratio of the number of binary strings tested by the heuristic algorithm to the total number 2 M S of mutants. We observed that the values Q i,s stabilize relatively early in the course of calculations. Testing mutants with at most three mutations is typically sufficient to stabilize the third digit after the decimal point. To the contrary, the entries of keep adjusting to the variations of D.
Since for S = 6 and M = 2 there are only 2 12 = 4096 binary strings involved, the corresponding likelihood maximization problem can be solved exactly with the IP method. We consider this case in order to illustrate efficiency of the heuristics and some of the complications discussed above. The threshold 0.001 remains. This threshold implies 64 significant macroeconomic scenarios according in the exact solution. Total probability of insignificant scenarios amounts to 72.59%. On the one hand, there are too many, for a conceptual analysis, significant macroeconomic scenarios. On the other hand, the total probability of all insignificant scenarios is too large.
The heuristic solutions was evaluated according to the scheme explained above. The four building blocks and their probabilities are given in Table 14.   Since there are 360 binary strings in V 2 in this case, sampling 90 mutants for every block-structure, random samples from V n , n = 3, 4, 5 were generated, each containing 360 binary strings. Using penalty values 5 and 10, the approximate distribution D given in Table 15 was estimated. Probability assigned to all insignificant macroeconomic scenarios was 0.05%. Table 16 characterizes the relative standing of L H against L 0 as well as against the maximum value L * found with the exact algorithm. The heuristic solution appears to be just slightly worse than the exact one, being more convenient for a conceptual analysis of the corresponding macroeconomic scenarios.  The following migration matrices were estimated for S = 6 and M = 7: