Simulating Historical Ratings Transition Matrices for Credit Risk Analysis in Mathematica

Mark S. Coleman
Director of Research, Financial Services Consulting Practice

Quantitative methods for evaluating credit risk have gained importance the past several years. Many of the most common analytical credit risk procedures use historical ratings transition matrices published by the three major U.S. rating agencies. In many cases, however, the specific approaches used are too limiting in that the observed historical variation in the Markov transition probabilities is not fully utilized to evaluate the risk of credit instrument portfolios. We argue that this variation, which in part is caused by changes in the underlying economic and financial factors that determine credit market behavior, is an important element of true underlying risk of these financial instruments. This paper uses Mathematica to present an approach for simulating ratings transitions that accounts for the underlying  variation in credit transition behavior. We use a copula-based Monte Carlo approach to help account for the historical cross-correlation of ratings probability changes over a 20-year sample of data spanning 1981 - 2000.

. Introduction

Over the past decade, quantitative methods for evaluating, predicting, and managing the risk of portfolios of credit-based securities have become increasingly important for financial institutions.  This need for improved risk assessment is a direct consequence of several factors. First, financial institutions are subject to increasingly stringent formal regulatory requirements, such as The New Basel Capital Accord  []. Second, there has been an enormous increase in the use of credit derivatives and other credit-related instruments for financing and hedging  []. Third, investors in credit-sensitive instruments are demanding greater liquidity, lower trading costs, and more "transparency" for their portfolios — requirements that can only be met by the use of sophisticated valuation models. The end result has been a proliferation of credit risk analytis methods  [] [].

This paper uses Mathematica to develop a generalization of a widely used credit risk estimation procedure. Our method addresses an important limitation to the existing procedure in that we more fully account for the implicit variation seen in published historical ratings transitions. We argue such variations flow from fundamental economic factors and are an important source of underlying credit volatility that are not adequately captured by existing methods. Using published annual transitions data spanning the period 1981 - 2000, we employ copula functions to build a multivariate distribution of ratings probabilities based upon the assumption that each element of the mean transtion matrix has a unique (marginal) beta distribution. This allows for correlations across the transition probabilities — a feature consistent with the historical data. We also present a basic Monte Carlo procedure for drawing correlated random vectors from this multivariate distribution.

. Modeling Credit Risk Using Ratings Transition Matrices

Some of the most widely used credit risk methodologies employ ratings transition matrices  [], published regularly by the major U.S. credit rating agencies  []. In essence, these matrices depict the average probability of a rated credit security moving from one rating category to another over a predetermined time interval. The probabilities are tabulated using an historical database of ratings issuance conducted by each agency over a specific period of time. For each period of interest, e.g., one year, a cohort of existing credits is first defined and then tracked through the end of the period. Any changes in ratings issued by the agency over the period is noted and then used to update the proportion of credits in each cell of the matrix at the end of period. This method of tabulation leads to a rather obvious Markov interpretation; hence formal Markov methods are typically used to model transition behavior and credit risk.

For instance, consider the following simple hypothetical example  []. This example calculates the expected value for a single bond one year ahead. For the purposes of illustrating our method, we adopt the conventional Standard & Poor's ratings categories, AAA, AA, A, BBB, BB, B,CCC,D, where D denotes default and is an absorbing state. We begin with a $100 par value bond, having a 6% annual coupon and a five year maturity with a beginning-of-year rating of BBB. Nominally, An investor holding this bond would receive $6 at the end of each of the first four years. When the bond matures at the end of year 5, the investor recieves $106 (the $100 face value plus the final $6 coupon). The value of such a bond today is simply the present discounted value of this five year earnings stream, assuming some appropriate market discount rate. This simple calculation, however, ignores the credit risk of the bond. To determine the expected value of a risky bond, we need to calculate a discounted value for each possible coupon stream that the bond would pay should it migrate to another rating category (for simplicity we assume that only one transition can occur). This requires a vector of one-year transition probabilities and a vector of discount rates corresponding to each possible rating category for each year of the bond's life (such rates are referred to as forward zero curves). The following table presents these inputs

For instance, if the bond is upgraded from BBB to A status, then the value of the bond at the end of one year is simply

\(TraditionalForm\`V(
        BBB → 
          A)\  = \ \($6\  + \ $6\/\((1\  + \ 0.0372)\) + $6\/\((1\  + \ 0.0432)\)\^2 + $6\/\((1\  + \ 0.0493)\)\^3 + $6\/\((1 + 0.0532)\)\^4 = \ $108  .66\)\)

Table 1 repeats this calculation for each possible rating migration, yielding a set of future values for the bond (note: we use an industry-standard assumption about the recovery value of the bond should it migrate to the default state). Each of these values is then weighted by the corresponding migration probabilitiy,  giving an expected value of the bond at the end of one year. The volatility of the one year ahead value is then calculated as a standard deviation in the usual manner. In this case

E(V)  =  $107  .09,  σ(
        V)  =  $2  .99

. Economic Variation and Correlations in Ratings Transitions

Though the example shown above is admittedly simplified, it is not far from procedures used in practice. In particular, most transition-based credit risk models rely upon a single point estimate of transition probabilities. This has the obvious defect of ignoring the time series behavior of transition probabilities, as well as any cross-category correlations. This is all the more curious in that, within the context of borader macorecnomic business cycles, credit cycles are quite pronounced and exhibit strong empirical regularity and predictability over time  []. Credit quality is strongly determined by borrower earnings, which in turn are tied to broader macroeconomic factors such as sales, wages, productivity, intereset rates, exports, etc. The historical correlation between key credit market factors — including the values and volatilties of the assets underlying the credit obligation — and general economic indicators is well documented  [] []. In addition, credit market spreads —  differences between yields on comparable maturity bonds of different credit quality — and spread volatilities also have strong statistical ties to broader capital market factors  []. Another more subtle reason to expect economic correlation is that ratings analysts are also monitoring economic conditions and will periodically adjust rating crieteria to reflect their perceptions of future economic risk. As a result, a reasonable estimate of expected credit risk should reflect prospective changes in underlying economic factors.

The historical volatility of selected transition probabilities can seen in Figure 1. The graph depicts the time series of the "no-migration" (the main diagonal) probabilities take from the series of one-year transition matrices published by Standard & Poor's over the period 1981 - 2000. There is obvious cyclical variation over this period, which includes a full business cycle. The statistical evidence indicates significant first order autocorrelation for several ratings levels, as well as  significant correlation across several categories. This is the case for any combination each of the rating classes (a summary of the sample statistics and correlations are available from the author upon request)

Figure :  Diagonal transition probabilities by rating category

. A Note on Calculating Markov Generators

As noted above, transition matrices are often modeled using Markov processes. Published transition matrices, however, are not always well suited for valuing many types of credit securities. For example, transition matrices are typically available only for annual frequencies, with the shortest period being one year. Many credit instruments have maturities under one year, making accurate valuation problematic. In addition, since many published matrices are at best only approximations to true Markov processes  [] , complicating the valuation of longer-term securities as well since multi-year transition probabilities cannot always be accurately calculated by the product of the one-year transitions. In a recent article, Israel, Rosenthal, and Wei [] (IRW) identify a set of conditions under which a valid Markov generator exists and present an algorithm for estimating an approximate generator when a true generator does not exist. If P is a time-homogenous N x N Markov transition matrix, then a generator Q has the property that P(t) = exp(tQ). Using a Markov generator implies that proper transtion probabilities can be estimated for any t ≥ 0, enabling the evaluation of credit securites of any maturity. We have implemented the IRW procedure in Mathematica (available upon request) and applied it to each of the one-year transition matrices in our sample. The results that follow reflect this adjustment.

. Simulating Transitions

. Generating  Uncorrelated Transition Probabilties

Markov models of credit risk treat default rates as a random variable. Default rates are often modeled as gamma distributed random variates, which follows from the idea that the underlying number of defaults in any transition cell can be modeled as a Poisson arrival process. We generalize this approach by employing a Monte Carlo procedure , where we draw from a unique beta distribtution for  each transition probability. The beta distribution is well suited to modeling set of proportions, and can be shown to generalize the gamma distribution is certain circumstances  []. In the discussion that follows, we begin by assuming no correlation between probabilities; we then relax this assumption using a proecedure shown in Clemen and Reilly [] that employs a multivariate normal copula function to define a joint distribution of the probabilities using only marginal distributions and pairwise correlations.

The beta distribution is defined as

\(TraditionalForm\`F(
        x)\  = \ \(\(Γ(α + β)\)\/\(\(Γ(α)\)\ \(Γ(β)\)\)\) \(\(x\^\(α - 1\)\)(1 - x)\)\^\(\(β - 1\)\(
\)
          \)\)

for α > 0, β > 0, and 0<x<1, where Γ(▪) is the gamma function. The mean and variance of (3) can be written in terms of α and β as

\(TraditionalForm\`μ = \(α\/\(α + β\)\ and\ \ \ \ \ σ\^2 = \(α\ β\)\/\(\(\(\((α + β)\)\^2\) \((α + β + 1)\)\)\(
\)
          \)\)\)

By solving (4), we can parameterize the specific beta distributions that we need for the Monte Carlo simulation. Solving this system yields

betaParams = 
      Solve[{μ == α/(α + β),  σ2 == (α β)/((((α + β))^2) (α + β + 1))}, {α, β}]

As an example, consider the mean of the mean and standard deviation of the AAA⟶AAA transition probability. The corresponding beta parameters are

bp  =  
      Reverse[Flatten[
          betaParams   /.  {μ -> 0.92, σ2 -> 
                0.00222}]]

So to generate a random sample,

RandomArray[
          BetaDistribution[bp[[1]][[2]], bp[[2]][[2]]], {1000, 7}]; // Timing

To implement this on our data sample, we first calculated the mean one-year transition probabilities by averaging across each of the cells in our historical sample. We then calculated the variance of each cell in a similar manner.  The resulting matrices are

NumberForm[MatrixForm[meanfinalTrans], 2]

NumberForm[MatrixForm[varfinalTrans], 3]

Now calculate the beta parameters for each element of the mean transition matrices (note: warnings supressed)

betaParamsMat = 
      betaParams  /. {μ -> meanfinalTrans, σ2 -> 
            varfinalTrans}

Now extract the α and β parameters into separate lists:

αMat = betaParamsMat[[1]][[2]][[2]];
βMat = betaParamsMat[[1]][[1]][[2]];

At this point, we need to handle the indeterminate cases. W assign the indeterminate entries values α -> 0.0000001 and β->0.15. This will force the beta probabilities to be effectively zero.

αMat =  αMat  /.  
          Indeterminate -> 0.00000001;
βMat  =  βMat  /.  
          Indeterminate -> 0.15;

Recall from the example above that to value a specific bond we first consider its initial credit rating. This corresponds to a selecting a specific row of αMat and βMat. For this row, we calculate the the two beta parameters of each element (corresponding to one of the eight possible ratings to which the bond could migrate).  Thus we need to generate N trials for each of the elements of

\(\(ratingsList = {"AAA", "AA", "A", "BBB", "BB", "B", "CCC", "D"};\)\)

The following module generates the Monte Carlo draws from the set of beta distributions

transMonteCarlo[rating_, n_]  :=  
      Module[
                    {r, αVec, βVec}, 
                     r = 
          Flatten[Position[ratingsList, 
              rating]]; 
                     αVec = Flatten[αMat[[r]]];   βVec = 
          Flatten[βMat[[r]]]; 
                     Map[RandomArray[#, n] &, 
          MapThread[
            BetaDistribution, {αVec,  βVec}]]
                     ]

For instance, to generate 10,000 trials of the vector of migration probabilities for a bond initially rated "A",

y = transMonteCarlo["A", 10000]; // Timing

. Generating Correlated Transitions

As noted above, transition probabilities show evidence of  correlation across rating categories. In general, to simulate a vector of correlated variates requires specification of the underlying joint density. This can be an unweildy procedure if the univatiate marginal distributions are complex, or if the number of variates is large. An alternative procedure is suggested by Clemen and Reilly  [], in which an arbitrarily complex multivariate density function can be constructed using only marginal distributions and a standard measure of pairwise dependence among the random variables. Recently, copulas have emerged as an important tool to model correlations and dependence in quantitative finance  [] [] [].

The core of the Clemin and Reilly approach makes use of copula functions. Let  F((Subscript[x, 1]), …, (Subscript[x, n])) be a multivariate distribution function for random variables (Subscript[X, 1]), …, (Subscript[X, n]), with corersponding marginal distributions (Subscript[F, 1])((Subscript[x, 1])),…,(Subscript[F, n])((Subscript[x, n])). Then

F((Subscript[x, 1]), …, (Subscript[x, n])) = C[(Subscript[F, 1])((Subscript[x, 1])),…,(Subscript[F, n])((Subscript[x, n]))]

where C((Subscript[u, 1]),…,(Subscript[u, n])) is the joint distribution with uniform marginals. Under the assumption that each (Subscript[F, i]) and C are differentiable, then the corresponding joint probability density function f((Subscript[x, 1]), …, (Subscript[x, n])) can we written as

f((Subscript[x, 1]), …, (Subscript[x, n])) = (Subscript[f, 1])((Subscript[x, 1]))×⋯×(Subscript[f, n])((Subscript[x, n])) × c[(Subscript[F, 1])((Subscript[x, 1])),…,(Subscript[F, n])((Subscript[x, n]))]

where (Subscript[f, i])((Subscript[x, i])) is the density corresponding to the distribution (Subscript[F, i])((Subscript[x, i])), and c = \(TraditionalForm\`\(\(\ \)\(∂\^n C\)\)\/\((∂F\_1  … ∂F\_n)\)\)is called the copula function. The copula function c essentially "couples" the marginals to a joint distribution and incorporates information about the dependence among the underlying random variates. The copula approach has several important advantages over the conventional approach of expressing a joint density function as the product of conditional densities and marginal densities. First, arbirarily complex marginal distributions can be accomodated. Second, dependencies across the variates can often be expressed as pairwise correaltions, which are often simple for analysts to assess, even for large problems. Third, the copula-based approach is often computationally more efficent to sample from than drawing directly from the complete joint density, making it useful for Monte Carlo procedures.

To make this procedure operational, we need to define a copula that embeds the dependence among the transition probabilities. One the simplest and most common measures of statistical dependence is the Pearson product-moment correlation. Following Clemen and Reilly, we make use of the multivariate normal copula, (Subscript[c, N]), which makes use of the correlation matrix in the same was as the familiar multivariate normal distribution. Let R be an n x n positive definite real matrix giving the pairwise correlations of the (Subscript[X, i])'s, and let Φ denote the univariate standard normal distribution. If y = ((Subscript[y, 1]),…,(Subscript[y, n])) is an n-dimensional vector, the multivariate normal copula can be written

(Subscript[c, N])[Φ((Subscript[y, 1])),…,Φ((Subscript[y, n]))|R] = exp[-(y^T)((R^(-1))-I)y/2]/|R\( | \^\(1/2\)\)

where I is the n x n identity matrix. To develop the multivariate density, we use the normal inverse transformation (Φ^(-1))and define (Subscript[Y, i])=(Φ^(-1))[(Subscript[F, i])((Subscript[X, i]))] for i =1,…,n. Subsituting into (6) and (7) gives

f((Subscript[x, 1]),…,(Subscript[x, n])|R)  =(Subscript[f, 1])(
              (Subscript[x, 1]))×⋯×(Subscript[f, n])(
              (Subscript[x, n])) × exp{-\((\(Φ\^\(-1\)\)[\(F\_1\)(
                          x\_1)], …, \(Φ\^\(-1\)\)[\(F\_n\)(x\_n)], )\) × ((R^(-1))-I) × \((\(Φ\^\(-1\)\)[\(F\_1\)(
                          x\_1)], …, \(Φ\^\(-1\)\)[\(F\_n\)(x\_n))\^T/2}\)/|R((|)^(1/2))

To implement this method in Mathematica we first define a the marginal density function and CDF for each of the cells. Since each is a beta distribution, we can write

betaDensity[alpha_, beta_, x_]  :=  
        PDF[BetaDistribution[alpha, beta], x];
betaCDF[alpha_, beta_, x_]  :=  
        CDF[BetaDistribution[alpha, beta], x];

The first term in (8) is a product of the underlying beta marginals. Using the convention discussed above for a specific initial rating category, the following module calculates the marginal density product

densityProduct[rating_, x_] :=  
        Module[{r, αVec, βVec}, 
r = 
            Flatten[Position[ratingsList, 
                rating]]; 
αVec = 
            Flatten[αMat[[r]]];   βVec = 
            Flatten[βMat[[r]]]; 
Apply[Times,  
            MapThread[
              betaDensity, {αVec,  βVec, 
                x}]]
];

Next we define the input vector for the normal copula, using the inverse normal CDF

copulaInput[rating_, x_] :=  
        Module[{r, αVec, βVec}, 
r = 
            Flatten[Position[ratingsList, 
                rating]]; 
αVec = 
            Flatten[αMat[[r]]];   βVec = 
            Flatten[βMat[[r]]]; 
Print[
            MapThread[
              betaCDF, {αVec,  βVec, 
                x}]]; 
Map[
            Quantile[NormalDistribution[0, 1], #] &, 
            MapThread[
              betaCDF, {αVec,  βVec, 
                x}]]
];

Now set up the multivariate normal copula function, incorporating the correlation matrix R,

mvnCopula[y_, R_] := 
        Exp[(-y . (Inverse[R]  -  
                        IdentityMatrix[Dimensions[R][[1]]]) . 
                    y)/2. ]/Sqrt[Det[R]];

Finally, we define the complete multivaiate joint density for the initial rating category,

mvBeta[rating_, x_, R_]  :=  
        densityProduct[rating, x]  .  
          mvnCopula[copulaInput[rating, x], R];

Since our convention applies to a prespecified initial rating, we also need a way of calculating R corresponding to each initial rating state. Assuming the historical data is kept in finalTransList, then

corrCopula[rating_]  := 
        Module[{r}, 
r = 
            Flatten[Position[ratingsList, 
                rating]]; 
CorrelationMatrix[
            Map[Extract[#, r] &, 
              finalTransList]]
];

For instance, for the AAA rating, the historical correlation matrix is

corrCopula["AAA"]

Using  mvBeta[ ] it is straightforward to evaluate any function of the multivariate distribution, such as conditional densities, expected values, etc. For purposes of generating correlated beta variates for a Monte Carlo simulation of transition probabilities, we generate a vector of multivariate standard normal variates with correlation matrix R, apply the standard normal distribution function Φ(▪) to each element of the vector, and then transform these elements to the appropriate beta distrubution using (Subscript[F, i]^(-1))[Φ(▪)] (see the module copulaInput[ ]).