Portfolio optimization of credit risky bonds: a semi-Markov process approach

This article presents a semi-Markov process based approach to optimally select a portfolio consisting of credit risky bonds. The criteria to optimize the credit portfolio is based on l∞-norm risk measure and the proposed optimization model is formulated as a linear programming problem. The input parameters to the optimization model are rate of returns of bonds which are obtained using credit ratings assuming that credit ratings of bonds follow a semi-Markov process. Modeling credit ratings by semi-Markov processes has several advantages over Markov chain models, i.e., it addresses the ageing effect present in the credit rating dynamics. The transition probability matrices generated by semi-Markov process and initial credit ratings are used to generate rate of returns of bonds. The empirical performance of the proposed model is analyzed using the real data. Further, comparison of the proposed approach with the Markov chain approach is performed by obtaining the efficient frontiers for the two models.


Introduction
The problem of credit risk, which consists of finding the likelihood of default of an obligor going into debt, is one of the most important problems in financial world. The credit risk, also known as default risk, is the risk of lender that borrower may not be able to meet its debt obligations. Credit risk analysis consists of finding the likelihood of default of an obligor going into debt. Credit risk models are basically divided into reduced-form models and firm value models (also known as structural models). Firm value models consider the model in Merton (1974) as the base model, which gives a mechanism of default in terms of the relation between liabilities and assets at maturity time T. On the other hand, reduced-form models do not specify the actual mechanism of default but model it as a non-negative random variable with distribution depending on the economic co-variables, interested readers can refer Duffie and Singleton (2012). There are many parameters associated with bond issuer or bond itself, which quantify the credit risk associated with it. Credit rating is one of the important parameters. Credit rating of a credit risky bond, assigned by a company, is an evaluation of its likelihood of default and ability to pay back the loan. Better the credit rating of a bond, safer it is. Banks and the firms that issue bonds are most concerned to quantify the default risk. International organizations like Standard & Poor's assign ratings to the companies that issue bonds in order to evaluate the credit risk. A credit rating is given to each company that issues a bond, specifying its capacity to repay the debt. Clearly, a higher interest rate is expected from a firm whose rating is lower. Credit ratings have become important since they are used as an input in many models for calculating economic capital for banks. Therefore, estimating the transition matrices of credit ratings is at the core of risk management and has applications in pricing derivatives and credit portfolio optimization.
Models for credit quality based on ratings are also frequently used in the pricing and risk management literature, see for example Jarrow et al. (1997); Kijima and Komoribayashi (1998); Akutsu et al. (2003). In 1997, Jarrow et al. (1997 applied for the first time Markov processes to capture the time evolution of credit ratings. These models are called "migration models". One of the drawbacks of this model is that it gives, in small time interval, zero probability of default to bonds with high credit ratings. Other papers (Hu et al. 2002;Nickell et al. 2000;Baillo and Fernandez 2007;Grimshaw and Alexander 2011)) followed the same approach to generate the transition matrix. More recent contributions include the quantification of the impact of business cycles on the dynamics of credit ratings and corresponding computation of conditional migration matrices (Boreiko et al. 2018), the proposal of parsimonious higher order Markov chains models able to reproduce downgrade rating momentum effects (Baena-Mirabete and Puig 2018) and the modeling of credit default data with intensity-based model stemmed from hidden Markov chains (Yu et al. 2017;Yu et al. 2019). In 2017, Dharmaraja et al. (2017) proposed a Markov chain model with catastrophe to determine the mean time of default of a risky asset. Centanni et al. (2017) modelled the dynamics of defaults for a dynamic set of firms in a given time period under the Markovian assumption. They assumed that there are some observable variables such as the total number of defaults, the total number of firms operating in the market at time t and some unobservable variables such as the number of firms alive or defaulted in each class at time t. Then, they using particle filtering techniques they obtained an approximation of the distribution of unobserved variables. Later, Tardelli (2018) extended this idea to find the probabilistic prediction of the actual partition of the population, and of the conditional distribution of the distance to defaults. In other direction on empirical analysis, Kou et al. (2014) studied the effectiveness of multiple criteria decision making methods in evaluating clustering algorithms by performing experiments on real-life credit risk and bankruptcy risk data sets.
In the papers by Nickell et al. (2000); Kavvathas (2001); Lando and Skødeberg (2002), appropriateness of Markov process for the description of accurate rating dynamics was addressed. Carty and Fons (1994); Nickell et al. (2000) proved that the current rating of a company depends not only on its last rating but on all the previous ones, the effect called rating momentum. Thus, the evolution of credit rating dynamics is non-Markov. Carty and Fons (1994); Duffie and Singleton (2012), showed that a complete information of the time spent inside the states is of major interest in the credit risk problem. The credit migration probability depends on the time spent by a company in a particular rating. In a continuous-time homogeneous Markov chain, durations in states follow an exponential distribution. An exponential distribution has a constant hazard function. But for credit rating dynamics, the hazard function is not constant. Hence, Markov model for credit rating is not appropriate (Frydman and Schuermann 2008). The other issue is time depen-dence. It means that in general, transition probabilities tends to change with the state of the economy, being low during periods of economic expansion and high during recession. Rating evaluation at two different time points is different (Nickell et al. 2000) and hence the process describing the evolution of rating dynamics is time non-homogeneous.
The literature proposing models that addresses the non-Markov behavior of credit rating dynamics is of recent origin. Korolkiewicz and Elliott (2008) proposed a hidden Markov model assuming that the Markov chain governing the true credit quality evolution is hidden in noisy or incomplete observations about credit ratings. In the paper by D' Amico et al. (2005), they have considered the credit risk problem as a reliability problem. They applied time-homogeneous semi-Markov processes (SMP) to solve the first issue. The second issue has been solved by extending the state space by the same authors (D' Amico et al. 2016) and many further developments have been made in their model (D'Amico et al. 2011;D' Amico et al. 2010;D'Amico et al. 2012). In 2017, Pasricha et al. (2017) proposed a credit rating model based on a Markov regenerative process which is a generalization of semi-Markov process model.
By optimization of portfolio consisting of credit risky bonds, we mean to find the optimal allocation of the wealth to the risky bonds in the portfolio such that the overall credit risk is minimum. In order to formulate a optimization problem, we need to define the risk measure to be used. For example, in equity portfolio optimization, the pioneering work by Markowvitz (1952) considers the variance of portfolio return as the risk measure. Konno and Yamazaki (1991) introduced mean absolute deviation as measure of risk and Young (1998) considered l ∞ norm as a risk measure. In credit portfolio optimization, value at risk was considered to quantify the portfolio credit risk in Morgan (1997). Later, Andersson et al. (2001) proposed a credit portfolio optimization model based on conditional value at risk (CVaR). Further Kalkbrener et al. (2007) considered risk measures VaR and CVaR (referred therein as Expected shortfall) for credit risk optimization formulation with importance sampling technique for calculation of risk measures. Their formulation overcomes the numerical problems associated with calculation of CVaR allocation in credit portfolio models. The other models proposed in the literature considered l 1 norm i.e., mean absolute deviation and l 2 norm i.e., variance as a measure of risk. In (2017), Singh and Dharmaraja (2017) proposed a l ∞ based credit portfolio optimization using credit rating dynamics modeled by a Markov chain. In 2015, Ma et al. (2015) proposed a model based on extreme value theory to evaluate the default risk of bond portfolios.
On similar lines to Singh and Dharmaraja (2017), this article considers a l ∞ norm i.e., min-max absolute deviation credit portfolio optimization model. However, due to the limitations of the Markov chains to model the credit rating dynamics, we assume that the credit rating dynamics of the risky bonds are assumed to follow a semi-Markov process and propose new risk premium adjustments for the pricing purposes. The transition probabilities of the rating migration is obtained using the historical data. Based on the transition probabilities, future credit ratings are obtained and hence following Jarrow et al. (1997), the price paths of credit risky bonds are obtained. These generated paths of returns of bonds are input parameters to the optimization model. Extensive numerical illustrations are given to show the applicability of the proposed model and is compared with Markov chain model proposed by Singh and Dharmaraja (2017). This article uses semi-Markov process to obtain the transition probabilities of the rating migration, how-ever, one can consider the models by Centanni et al. (2017); Tardelli (2018) to obtain the dynamics of credit quality of firms and generate the path of bond returns.
The rest of the paper is organized as follows. "The proposed framework" section introduces the min-max absolute deviation model followed by the semi-Markov credit rating model. "Proposed methodology" section proposed the algorithm to generate the future credit rating scenarios credit rating followed by the methodology to solve the portfolio optimization model based on the generated ratings. Numerical illustrations are given in "Empirical analysis" section, and "Conclusions and future work" section concludes the paper with future work.

The proposed framework
In the first subsection, we give a brief review of the l ∞ norm risk measure based portfolio optimization model. In the next subsection, a semi-Markov credit rating model proposed in D' Amico et al. (2005) is presented. In the reminder of the paper, we shall be using the following notations: N Number of bonds in the portfolio x n dollar amount spent on nth bond in the portfolio T length of time horizon t each period over the time horizon t = 1, 2, . . . T r j the expected return E(R j ) of asset j r jt the observed return of asset j during the period t C the total portfolio expenditure d user defined value minimum rate of return required by investor.
Note that here we assume that the bonds are zero coupon bonds. l ∞ or Min-Max absolute deviation model Young (1998) proposed a l ∞ -norm risk measure based portfolio optimization model of equities. Following this approach, Cai et al. (2000) considered maximum of absolute deviation of individual asset's return as measure of risk, that is l ∞ -norm risk measure. In 2017, Singh and Dharmaraja (2017) proposed a l ∞ based credit portfolio optimization. In this section, we present the mathematical formulation of minimizing the maximum absolute deviation as follows: The objective here is to minimize the maximum of the absolute deviation of portfolio return. In other words, instead of minimizing mean of the absolute deviations, we propose to minimize the worst possible performance of the portfolio. This is more conservative portfolio selection rule and is suitable for risk-averse investors. Further, the first constraint restrict the mean portfolio to be greater than a threshold percentage d of the total budget C, given by the second constraint. Finally, the last constraint indicates the restriction on short selling. This original problem can be reformulated into a linear programming problem using auxiliary variable y t as follows: Let Let us assume that Y = max t=1,2,...,T y t . This implies that for all t = 1, 2, . . . , T, we have Using these transformations we get the following linear programming problem:

Semi-Markov credit rating model
The level of credit rating changes from time to time because of random credit risks and thus need to be modeled by an appropriate stochastic process. Jarrow et al. (1997) The kernel Q(t) =[ Q i,j (t)] associated with the process is defined by and it follows that where P =[ p ij ] i,j∈ is the one-step transition probability matrix of the embedded discrete time Markov chain with state space . The probability that the process will leave state i in time t is given as, It can be observed that Next, consider the distribution function of the waiting time in each state i, given that the next state is known These probabilities can be obtained as follows The main difference between a continuous-time Markov chain and an SMP is the distribution functions G ij (t). In a Markov environment this distribution function has to be cumulative distribution function a negative exponential distribution. On the other hand, in the semi-Markov case the distribution functions G ij (t) can be cumulative distribution of any general distribution. Thus accounting for the effect of duration inside a rating class. Now, a time-homogeneous semi-Markov process {Z(t), t ≥ 0}, which represents, for each waiting time, the state occupied by the process is defined as The transition probabilities for {Z(t), t ≥ 0} are defined by They can be obtained by solving the Markov renewal equation (Kulkarni 2016) where δ ij represents Kronecker delta. This equation can be solved numerically using discretization to numerically evaluate the integrals. Let h > 0 be the step size of discretization, then we have the countable linear system given by In the credit risk environment, the first part of above equation can be interpreted as the probability of the firm to remain in rating i from time 0 to t given that the rating organization gave a new rating evaluation at time 0. In the second part of above equation, φ γ j (t − τ )q iγ (τ ) represents the probability that firm will get a rating γ in time τ given that starting at time 0 from state i and then firm will migrate to rating j in remaining time (t − τ ) following one of the possible paths. φ i,j (t) are the actual transition probabilities of the observed discrete time-homogeneous semi-Markov process.

Proposed risk premium adjustments
For the pricing purposes, we need to define the probabilities under risk neutral measure. In this section we consider a change of measure that is compatible with the general theory developed by Vasileiou and Vassiliou (2006): whereQ ij (t) is the kernel in the risk neutral world and i (t) are the risk premium adjustments. Since Q(t) is only asymptotically stochastic, we have the following constraint lim t→∞ j∈ This implies that where lim t→∞ i (t) = * i . Hence, we proposed the following change of measure for t ≥ 0, Since we need to obtainQ(·) as a semi-Markov kernel, then some restrictions should be considered for the function (t) other than lim t→∞ i (t) = * i . More precisely, 1Q ij (0) = 0 ∀i, j ∈ . This means that is an increasing function of t, i.e., This means that Similarly, for i = j, we haveQ This means that Summarizing, we have that the function i (t) should satisfy the following conditions These risk premium adjustments can be obtained from the available market prices of risk free bonds and credit risky bonds of all the rating classes and of all the maturities..

Valuation of bond portfolio
Let μ 0 (s, T) and μ j (s, T) denotes the price of risk-free discount bond and the price of risky discount bond with rating j at time s and with maturities T respectively. The price of risky bond is given by Jarrow et al. (1997) whereẼ j,s denotes the conditional expectation (under the risk neutral probability measurẽ P), given the information that at time s rating is j, τ j represent the time of default of credit rating processZ whenZ s = j, I A denotes the indicator function. It follows that is the probability that no default occurs till date T given the information that at time s rating is j, δ is the recovery rate of the risky discount bond. The claim holders get δ amount at the maturity of the contract if default occurs and receive face value if no default occurs. Here, it is assumed that recovery rate is a constant. If s ≥ T, and if the credit rating of bond at time s i.e., Z s = j is known, we have

Proposed methodology
Throughout this section, it is assumed that the present time is s and the risk horizon, i.e., the time at which portfolio is evaluated is t = s + 1. Let Z n (s) denote the credit rating of asset n, n = 1, 2, . . . , N at time s and assume that they are independent to each other. Although the assumption of independence among the credit ratings of the bonds is very simplistic, several articles in the literature have justified this assumption. For instance, Kalkbrener et al. (2007) argued that the credit ratings of the bonds can be assumed uncorrelated if they are from different sectors of the financial markets. Incorporating the correlation among the credit ratings of the bonds could be a possible extension of the proposed framework.

Algorithm to generate the credit rating scenarios
In this paper, using Monte Carlo method, we generate future credit rating dynamics of a portfolio of credit risky bonds with the assumption that credit rating dynamics of bonds in the portfolio are independent of each other. Let Z = (Z 1 (t), Z 2 (t), . . . , Z N (t)) be the N-tuple vector where Z k (t) is the credit rating of the kth bond at time t. The algorithm for random sample generation of N-tuple credit rating vector, for time periods t = s + 1, s + 2, . . . , T, is as follows: Step 1: Set t = s + 1. Given initial credit rating Z n (s) = z n s , find the pmf and hence CDF of nth bond as follows:

Hence, the CDF of nth bond is given by
where F z n s ,i (t) represent the probability of nth bond being in state less than equal to i at time t. Hence, pmf and CDF of nth bond at time 1 can be calculated once we know the initial rating and transition probability function.
Step 4: Set t = s + 2 and following the Step 1, obtain the pmf and CDF of each bond as in Step 1.
Step 6: Increment t to t + 1. Repeat Step 4 to get the next period ratings.
Step 7 If t + 1 = T, then stop. Otherwise repeat Steps 4, 5 and 6. Repeat the mentioned algorithm until enough number of scenarios are obtained. Let {Z n k (r); n = 1, 2, . . . , N, k = 1, 2, . . . , S} be the credit rating of nth bond generated at time r in the kth scenario. Each generated scenario gives N × 1 vector of credit ratings of N bonds.

Methodology
This subsection explains the methodology to obtain the optimal portfolio of credit risky assets. Let S be the number of simulated scenarios of joint credit rating process.
1 Estimate the one-step transition probability matrix φ(1) =[ φ i,j (1)] i,j∈ of continuous time-homogeneous semi-Markov model using the risk neutral measure. 2 Given initial rating z j (0) of j th asset, j = 1, 2, . . . , N, find the one-step condition CDF as follows 3 Simulate S number of scenarios of the credit ratings following the methodology in "Algorithm to generate the credit rating scenarios" subsection. We have N × S matrix with (j, k) entry representing rating of j th asset at time 1 in k th scenario i.e., z j k (1). 4 Given z j k (1) and z j (0), we obtain the price of bond at time 0 and 1 i.e. μ z j (0) (0, T) and μ z j k (1) (1, T) for each j = 1, 2, . . . , N, k = 1, 2, . . . , S using Eq. (4). 5 Find return of j th asset in k th scenario as follows where μ z j k (1) (1, T) are obtained in Step 4 above. 6 Repeat the above process by moving the window 1 step ahead (from time 1 to time 2 and so on) until T time is reached. 7 Solve the optimization problem (P) o and obtain the weights of each asset.

Empirical analysis
In this section, we illustrate the applicability of the proposed model using real data. We consider two different sectors namely industry sector and service sector. Further, we consider 10 credit risky bonds with same maturity of 9 years and these bonds are grouped into two sets of five bonds each and these two sets comes from two different sectors. We assume that the five bonds belonging to each sector have same transition probability matrix that is obtained assuming they follow the semi-Markov process. However, they differ in their initial credit rating (Tables 3 and 4). Similarly, we assume that the two bonds belonging to different sectors have different credit rating dynamics both obtained using semi-Markov process. We assume that the recovery rate for each of the bond is 0.
Considering the real data of historical ratings for two sectors, we estimate the parameters of the two semi-Markov models for two sectors. After estimating the parameters, the nth step transition probability matrices are obtained solving the Markov renewal equation. For instance, 2 step transition probability matrix for two sectors are given in Tables 1 and 2. Using these transition probabilities matrices, we obtain CDFs for each rating category as discussed in the previous sections. We obtain 1000 scenarios of the future credit ratings for 10 bonds and for each period t = 1, 2, 3, . . . , 9 using CDFs obtained for the two sectors. Using the simulated scenarios, we obtain the price of the credit risky bonds and hence we find the returns. In the method of obtaining price of a bond for a given credit rating state, we need the transition probabilities in a risk neutral world. We apply the proposed change of measure to the historical kernel estimated from the data and obtain risk neutral transition probabilities. For the illustration purpose we have considered current yield for all the credit risky bonds and the risk free government bond from Akutsu et al. (2003) and risk premium adjustments from Kijima and Komoribayashi (1998). Average returns for each period is obtained taking the average of the returns simulated for the two sectors over all the scenarios. The average returns obtained are shown in Tables 3 and 4.  We obtain the optimal portfolio using the R software solving the min-max absolute deviation optimization problem. For the comparison purpose, we also obtain the future credit ratings scenarios and hence returns of credit risky bonds using the Markov chain model. The optimal weights of the assets are obtained by solving the optimization problem. The efficient frontier of min-max absolute deviation optimization model is shown in Fig. 1. It is obtained by solving the problem for different values of minimum return d. It is evident from Fig. 1 that the returns are better in case of semi-Markov model as compared to the Markov chain model. Further, for the comparison purpose, the average return, average risk and the average Sharpe ratio (the returns per unit risk taken) are obtained for both the portfolios and are shown in Table 5. It is evident that semi-Markov model performs better than Markov model with respect to all three measures.
The model possesses relevant practical managerial implications. First, we notice that recent investigations proved that credit bond ratings can be jointly considered with market measures to optimally select financial portfolios, see Choi et al. (2020)). This strategy may overperform classical portfolios choices based only on market variables. We conjecture that the benefits may increase with the adoption of adequate rating models and the semi-Markovian proposal appears to be justified by credit rating data. Indeed, the results of the evaluation of the bond portfolio acknowledge that if the semi-Markov processes are incorporated to model the credit rating of the firms, the portfolios obtained performs  better in comparison to the model where Markov chains are used (Singh and Dharmaraja 2017). Second, the evidence denotes that credit rating dynamics exhibit durations effects and accordingly managers could track important benefits by incorporating credit bond duration information in the portfolio selection procedure by avoiding for example too early or even unnecessary bond reallocations. Further, the other contribution of this article is the new risk premium adjustments, proposed for the pricing purposes, based on the general theory developed by Vasileiou and Vassiliou (2006). Since the credit risk management is one of the most important problems for investors in practical risk management, this article may be of great interest to the investors as it explicitly incorporates a firm's credit rating, modelled through semi-Markov processes, for valuing debt securities.
Finally, the proposed model can prove useful to banks and other financial intermediaries most concerned about the portfolio credit risk evaluation.

Conclusions and future work
In this article, we considered a portfolio optimization problem when portfolio consists of credit risky bonds. We considered l ∞ norm as the risk measure. To obtain the future returns of the bonds that are input to optimization model, we followed the approach of Jarrow et al. (1997) using credit ratings. In order to obtain the future credit ratings, we applied a semi-Markov credit rating model by D' Amico et al. (2005). The risk premium adjustments to obtain the risk neutral transition probabilities are proposed. A detailed algorithm is given to generate the credit ratings followed by the methodology to obtain the optimal portfolio using credit ratings. Extensive numerical analysis is given based on the real data in order to show the applicability of the model. Further, the returns obtained from the two models are compared based on the average return, average risk and average Sharpe ratio. One direction of the further extension of the present work is the incorporation of correlated credit rating dynamics of bonds.
Abbreviations SMP: Semi-Markov processes; CVaR: Conditional value at risk; VaR: Value at risk; CDF: Cumulative distribution function; pmf: Probability mass function

Authors' contributions
We have no conflicts of interest to disclose. All the authors contributed equally to this work. All authors read and approved the final manuscript.

Funding
This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

Availability of data and materials
The datasets used and analysed during the current study are available from the corresponding author on reasonable request.

Competing interests
The authors declare that they have no competing interests