Robust and Consistent Estimation of Generators in Credit Risk

TL;DR

This paper improves the estimation of generator matrices in credit risk models by providing closed-form formulas for the EM algorithm, ensuring convergence, and demonstrating superior performance over existing methods in benchmark tests.

Contribution

It introduces directly computable formulas for the EM algorithm in CTMCs, proves convergence under weak conditions, and benchmarks the method showing improved accuracy and efficiency.

Findings

Faster computation of generator matrices with closed-form formulas.

Proven convergence of the EM algorithm under weak conditions.

Outperforms existing algorithms in credit risk benchmarks.

Abstract

Bond rating Transition Probability Matrices (TPMs) are built over a one-year time-frame and for many practical purposes, like the assessment of risk in portfolios or the computation of banking Capital Requirements (e.g. the new IFRS 9 regulation), one needs to compute the TPM and probabilities of default over a smaller time interval. In the context of continuous time Markov chains (CTMC) several deterministic and statistical algorithms have been proposed to estimate the generator matrix. We focus on the Expectation-Maximization (EM) algorithm by Bladt and Sorensen (2005) for a CTMC with an absorbing state for such estimation. This work's contribution is threefold. Firstly, we provide directly computable closed-form expressions for quantities appearing in the EM algorithm and associated information matrix, allowing to easily approximate confidence intervals. Previously, these quantities had to be estimated numerically and considerable computational speedups have been gained. Secondly, we prove convergence to a single set of parameters under very weak conditions (for the TPM problem). Finally, we provide a numerical benchmark of our results against other known algorithms, in particular, on several problems related to credit risk. The EM algorithm we propose, padded with the new formulas (and error criteria), outperforms other known algorithms in several metrics, in particular, with much less overestimation of probabilities of default in higher ratings than other statistical algorithms.

Tables

Table 1. Table 4.1: True unstable generator

	AAA	AA	A	BBB	BB	B	C	D
AAA	-0.146371	0.085881	0.04549	0.015	0	0	0	0
AA	0.018506	-0.166337	0.114831	0.033	0	0	0	0
A	0.0276	0.047012	-0.198043	0.09043	0.023001	0.01	0	0
BBB	0.011469	0.010734	0.088133	-0.243046	0.077569	0.044407	0.010734	0
BB	0	0	0.019159	0.184699	-0.323077	0.106166	0.013053	0
B	0	0	0.012280	0.034822	0.093489	-0.296265	0.134273	0.022401
C	0	0	0	0	0.02	0.140209	-0.600939	0.440730
D	0	0	0	0	0	0	0	0

Table 2. Table 4.2: True stable generator

	AAA	AA	A	BBB	BB	B	C	D
AAA	-0.061371	0.055881	0.005490	0	0	0	0	0
AA	0.013506	-0.096337	0.074831	0.008	0	0	0	0
A	0	0.037012	-0.097442	0.06043	0	0	0	0
BBB	0	0.000734	0.058133	-0.120843	0.057569	0.004407	0	0
BB	0	0	0.009159	0.104699	-0.190024	0.076166	0	0
B	0	0	0	0.024822	0.083489	-0.174985	0.064273	0.002401
C	0	0	0	0	0	0.080209	-0.300939	0.220730
D	0	0	0	0	0	0	0	0

Table 3. Table 4.3: Order of time taken to execute the various algorithms. Note that MCMC also depends on the level of information i.e. obligors in each rating. We also note that BS 09 algorithm is faster than the other MCMC algorithms but still takes 10 4 superscript 10 4 10^{4} seconds in the case of 1000 obligors per rating.

Algorithms	Deterministic	EM	MCMC
Time (seconds)		$\sim 10$	$\sim {10}^3$ to $\sim {10}^4$

Table 4. Table 4.4: Observed TPM used to estimate the generators in probability of default plots.

	AAA	AA	A	BBB	BB	B	C	D
AAA	0.8824	0.1176	0	0	0	0	0	0
AA	0.0064	0.9111	0.0813	0.0008	0.0001	0	0.0003	0
A	0.0003	0.0559	0.8836	0.0499	0.0079	0.0015	0.0002	0.0007
BBB	0	0.0116	0.1585	0.7640	0.0528	0.0070	0	0.0061
BB	0	0	0.0213	0.1193	0.7746	0.0623	0.0099	0.0127
B	0	0	0.0062	0.0199	0.1669	0.7017	0.0730	0.0322
C	0	0	0	0	0.0417	0.2083	0.4544	0.2956
D	0	0	0	0	0	0	0	1

Table 5. Table 4.5: Mixed portfolio

AAA	100, 500, 1500, 750
AA	200, 750, 2000, 650
A	150, 400, 400
BBB	300, 500, 150, 1500
BB	500, 250, 700
B	200, 500
C	100, 150, 200

Table 6. Table 4.8: Risk charge results for the true generators.

		Stable			Unstable
	Mixed	Investment	Speculative	Mixed	Investment	Speculative
IRC	702	0.32	3395	1251	0.41	5057
IRC ES	508	0.20	2409	842	3.78	3826
IDR	750	0	3400	1750	200	4600

Table 7. Table 4.9: Risk charge results for each algorithm as a %.

			Stable			Unstable
		Mixed	Investment	Speculative	Mixed	Investment	Speculative
	EM	7.3	7.5	1.5	22.5	29 195	2.6
	DA	11.9	8.1	2.4	36.9	66 829	4.3
	WA	11.8	8.1	2.3	37.3	69 293	4.1
IRC	QOG	11.6	7.8	2.3	26.7	38 976	4.1
	MCMCBS05	154	306 000	2	49.6	478 000	4.1
	MCMCBS09	24.9	18.4	14.4	68.3	264 000	14
	MCMCMode	12.5	8.1	3.6	34.9	39 000	3.9
	EM	5.3	115	3.4	8.6	375	2.7
	DA	8.2	235	5.1	16.6	1130	3.9
	WA	7.8	210	5	16.4	1109	3.8
IRC ES	QOG	7.3	123	4.9	12.6	622	3.8
	MCMCBS05	35.4	135 000	4.7	19.7	5315	4.1
	MCMCBS09	21	610	15.5	67.7	6693	13.1
	MCMCMode	9.2	235	6.1	19.1	1063	3.5
	EM	6	0	0.3	4.3	113	3.5
	DA	10	0	1.2	8.6	295	5.7
	WA	9.3	0	0.6	8.6	295	5.2
IDR	QOG	7.3	0	0.6	5.4	185	5.3
	MCMCBS05	139	1580	0.3	12.6	530	4.7
	MCMCBS09	20	40	9.3	33.7	775	13.2
	MCMCMode	10	10	0.9	8	278	4.7

Equations

$$g(y;\Psi )={\int }_{\mathcal{X}(y)}f(x;\Psi )dx.$$

$$\displaystyle R(\Psi^{\prime};\Psi):=\mathbb{E}_{\Psi}\Big{[}\ln\big{(}f(x;\Psi^{\prime})\big{)}\big{|}y\Big{]}\,\qquad\textrm{for }\leavevmode\nobreak\ \Psi^{\prime},\Psi\in\Lambda,$$

$$L_t(\mathbf{Q})=\exp \left\{\sum _{i=1}^h\left[\sum _{j\ne i}\log (q_{ij})K_{ij}(t)-S_i(t)\sum _{j\ne i}{q}_{ij}\right]\right\},$$

$$R({\mathbf{Q}}^{\prime };\mathbf{Q})=\sum _{i=1}^h\left[\sum _{j\ne i}\log (q_{ij}^{\prime }){\mathbb{E}}_{\mathbf{Q}}[K_{ij}(t)|y]-{\mathbb{E}}_{\mathbf{Q}}[S_i(t)|y]\sum _{j\ne i}{q}_{ij}^{\prime }\right],$$

$$q_{ij}^{\prime }=\frac{{\mathbb{E}}_{\mathbf{Q}}[K_{ij}(t)|y]}{{\mathbb{E}}_{\mathbf{Q}}[S_i(t)|y]}.$$

$$\displaystyle V^{*}_{ij}(\mathbf{c},\mathbf{Z};t)=\mathbb{E}_{\mathbf{Q}}\left[\exp\left\{-\sum_{\mu=1}^{h}c_{\mu}S_{\mu}(t)\right\}\prod_{\mu,\nu=1}^{h}Z_{\mu\nu}^{K_{\mu\nu}(t)}\mathbbm{1}_{\{X(t)=j\}}\biggl{|}X(0)=i\right]\,,$$

$${\mathbf{V}}^{\ast }(\mathbf{c},\mathbf{Z};t)=\exp \{[\mathbf{Q}\bullet \mathbf{Z}-\Delta (\mathbf{c})]t\},$$

$$\mathbf{C}^{(\alpha\beta)}_{\gamma}=\left[{\begin{array}[]{cc}\mathbf{Q}&q_{\alpha\beta}\mathbf{e}_{\alpha}\mathbf{e}_{\beta}^{\intercal}\\ 0&\mathbf{Q}\end{array}}\right]\qquad\text{ and }\qquad\mathbf{C}^{(\alpha)}_{\phi}=\left[{\begin{array}[]{cc}\mathbf{Q}&\mathbf{e}_{\alpha}\mathbf{e}_{\alpha}^{\intercal}\\ 0&\mathbf{Q}\end{array}}\right]\,\quad\alpha,\beta\in\{1,\cdots,h\}.$$

$${\mathbb{E}}_{\mathbf{Q}}[K_{ij}(t)|y]=$$

$${\mathbb{E}}_{\mathbf{Q}}[S_i(t)|y]=$$

$${\mathbb{E}}_{\mathbf{Q}}[K_{ij}(T)|\mathbf{P}]=$$

$${\mathbb{E}}_{\mathbf{Q}}[S_i(T)|\mathbf{P}]=$$

$$R({\mathbf{Q}}^{(k+1)};{\mathbf{Q}}^{(k)})-R({\mathbf{Q}}^{(k)};{\mathbf{Q}}^{(k)})\ge \sigma (||{\mathbf{Q}}^{(k+1)}-{\mathbf{Q}}^{(k)}||).$$

$$\displaystyle P^{u}_{ij}\frac{\epsilon q_{ij}}{h}\leq\mathbb{E}_{\mathbf{Q}}[K_{ij}(T)|\mathbf{P}]\leq h^{2}N\frac{ht}{\epsilon\min\big{\{}\epsilon^{h}t^{h}\exp\{-th^{2}/\epsilon\}\ ,\ \exp\{ht/\epsilon\}\big{\}}}\,.$$

$$\displaystyle\mathbb{E}_{\mathbf{Q}}[S_{i}(T)|\mathbf{P}]\geq P^{u}_{ii}t\exp\big{\{}-\frac{ht}{\epsilon}\big{\}}\,.$$

$$R\big{(}\mathbf{Q}^{(k+1)};\mathbf{Q}^{(k)}\big{)}-R\big{(}\mathbf{Q}^{(k)};\mathbf{Q}^{(k)}\big{)}\geq\lambda\|\mathbf{Q}^{(k+1)}-\mathbf{Q}^{(k)}\|^{2}\,,$$

$$\displaystyle\sum_{i=1}^{h}\sum_{j\neq i}\left[\mathbb{E}_{\mathbf{Q}^{(k)}}[K_{ij}(t)|\mathbf{P}]\big{(}\log(q^{(k+1)}_{ij})-\log(q^{(k)}_{ij})\big{)}-\mathbb{E}_{\mathbf{Q}^{(k)}}[S_{i}(T)|\mathbf{P}]\big{(}q^{(k+1)}_{ij}-q^{(k)}_{ij}\big{)}\right]\,.$$

$$\displaystyle\mathbb{E}_{\mathbf{Q}^{(k)}}[K_{ij}(T)|\mathbf{P}]\Big{(}$$

$$-{\mathbb{E}}_{{\mathbf{Q}}^{(k)}}[S_i(T)|\mathbf{P}]\left(q_{ij}^{(k+1)}-q_{ij}^{(k)}\right)\ge \lambda {\left(q_{ij}^{(k+1)}-q_{ij}^{(k)}\right)}^2,$$

$$\log (q_{ij}^{(k+1)})-\log (q_{ij}^{(k)})=\frac{-1}{q_{ij}^{(k+1)}}(q_{ij}^{(k)}-q_{ij}^{(k+1)})+\frac{1}{2Z^2}(q_{ij}^{(k)}-q_{ij}^{(k+1)}{)}^2,$$

$$\frac{{\mathbb{E}}_{{\mathbf{Q}}^{(k)}}[K_{ij}(T)|\mathbf{P}]}{2Z^2}(q_{ij}^{(k)}-q_{ij}^{(k+1)}{)}^2\ge \lambda (q_{ij}^{(k+1)}-q_{ij}^{(k)}{)}^2.$$

$$\frac{{\mathbb{E}}_{{\mathbf{Q}}^{(k)}}[K_{ij}(T)|\mathbf{P}]}{2Z^2}(q_{ij}^{(k)}-q_{ij}^{(k+1)}{)}^2\ge \frac{{\mathbb{E}}_{{\mathbf{Q}}^{(k)}}[S_i(T)|\mathbf{P}{]}^2}{2{\mathbb{E}}_{{\mathbf{Q}}^{(k)}}[K_{ij}(T)|\mathbf{P}]}(q_{ij}^{(k)}-q_{ij}^{(k+1)}{)}^2.$$

$$\frac{{\mathbb{E}}_{{\mathbf{Q}}^{(k)}}[K_{ij}(T)|\mathbf{P}]}{2Z^2}(q_{ij}^{(k)}-q_{ij}^{(k+1)}{)}^2\ge \frac{{\mathbb{E}}_{{\mathbf{Q}}^{(k)}}[K_{ij}(T)|\mathbf{P}]}{2{\left(q_{ij}^{(k)}\right)}^2}(q_{ij}^{(k)}-q_{ij}^{(k+1)}{)}^2.$$

$$\frac{{\mathbb{E}}_{{\mathbf{Q}}^{(k)}}[K_{ij}(T)|\mathbf{P}]}{2Z^2}(q_{ij}^{(k)}-q_{ij}^{(k+1)}{)}^2\ge \frac{P_{ij}^uϵ}{2h{q}_{ij}^{(k)}}(q_{ij}^{(k)}-q_{ij}^{(k+1)}{)}^2.$$

$$q_{ij}\rightarrow\text{sign}\big{(}\text{Re}(q_{ij})\big{)}\times|q_{ij}|\,,$$

$$\frac{{\partial }^{2}L(\Psi ;y)}{\partial {\Psi }^2}={\left[\frac{{\partial }^{2}R({\Psi }^{\prime };\Psi )}{\partial {\Psi }^{\prime 2}}+\frac{{\partial }^{2}R({\Psi }^{\prime };\Psi )}{\partial {\Psi }^{\prime }\partial \Psi }\right]}_{{\Psi }^{\prime }=\Psi }.$$

$$\frac{{\partial }^{2}R({\mathbf{Q}}^{\prime };\mathbf{Q})}{\partial q_{\alpha \beta }^{\prime }\partial q_{\mu \nu }^{\prime }}=\frac{-1}{q_{\mu \nu }^{\prime 2}}{\mathbb{E}}_{\mathbf{Q}}[K_{\mu \nu }(t)|y]{\delta }_{\alpha \mu }{\delta }_{\beta \nu },$$

$$\frac{{\partial }^{2}R({\mathbf{Q}}^{\prime };\mathbf{Q})}{\partial q_{\alpha \beta }\partial q_{\mu \nu }^{\prime }}=\frac{1}{q_{\mu \nu }^{\prime }}\frac{\partial }{\partial q_{\alpha \beta }}{\mathbb{E}}_{\mathbf{Q}}[K_{\mu \nu }(t)|y]-\frac{\partial }{\partial q_{\alpha \beta }}{\mathbb{E}}_{\mathbf{Q}}[S_{\mu }(t)|y],$$

$$\frac{{\partial }^{2}R({\mathbf{Q}}^{\prime };\mathbf{Q})}{\partial q_{\alpha \beta }\partial q_{\mu \nu }^{\prime }}=$$

$$\displaystyle\qquad\qquad\qquad+(e^{\mathbf{Q}(t_{s+1}-t_{s})})^{-1}_{y_{s},y_{s+1}}\left(e^{\mathbf{C}_{\psi}^{(\alpha\beta,\mu\nu)}(t_{s+1}-t_{s})}\right)_{y_{s},3h+y_{s+1}}\biggr{]}$$

$$\displaystyle\qquad-\biggl{[}-(e^{\mathbf{Q}(t_{s+1}-t_{s})})^{-2}_{y_{s},y_{s+1}}\left(e^{\mathbf{C}^{(\alpha\beta)}_{\eta}(t_{s+1}-t_{s})}\right)_{y_{s},h+y_{s+1}}(e^{\mathbf{C}^{(\mu)}_{\phi}(t_{s+1}-t_{s})})_{y_{s},h+y_{s+1}}$$

$$\displaystyle\qquad\qquad\qquad+\,(e^{\mathbf{Q}(t_{s+1}-t_{s})})^{-1}_{y_{s},y_{s+1}}\left(e^{\mathbf{C}_{\omega}^{(\alpha\beta,\mu)}(t_{s+1}-t_{s})}\right)_{y_{s},3h+y_{s+1}}\biggr{]}\,,$$

$$\mathbf{C}^{(\alpha\beta)}_{\gamma}=\left[{\begin{array}[]{cc}\mathbf{Q}&q_{\alpha\beta}\mathbf{e}_{\alpha}\mathbf{e}_{\beta}^{\intercal}\\ 0&\mathbf{Q}\end{array}}\right]\leavevmode\nobreak\ ,\leavevmode\nobreak\ \mathbf{C}^{(\alpha)}_{\phi}=\left[{\begin{array}[]{cc}\mathbf{Q}&\mathbf{e}_{\alpha}\mathbf{e}_{\alpha}^{\intercal}\\ 0&\mathbf{Q}\end{array}}\right]\leavevmode\nobreak\ ,\leavevmode\nobreak\ \mathbf{C}^{(\alpha\beta)}_{\eta}=\left[{\begin{array}[]{cc}\mathbf{Q}&\mathbf{e}_{\alpha}\mathbf{e}_{\beta}^{\intercal}-\mathbf{e}_{\alpha}\mathbf{e}_{\alpha}^{\intercal}\\ 0&\mathbf{Q}\end{array}}\right]\,,$$

Full text

This work’s contribution is threefold. Firstly, we provide directly computable closed form expressions for quantities appearing in the EM algorithm and associated information matrix, allowing to easily approximate confidence intervals. Previously, these quantities had to be estimated numerically and considerable computational speedups have been gained. Secondly, we prove convergence to a single set of parameters under very weak conditions (for the TPM problem). Finally, we provide a numerical benchmark of our results against other known algorithms, in particular, on several problems related to credit risk. The EM algorithm we propose, padded with the new formulas (and error criteria), outperforms other known algorithms in several metrics, in particular, with much less overestimation of probabilities of default in higher ratings than other statistical algorithms.

Keywords: Likelihood inference, Credit Risk, Transition Probability Matrices, EM algorithm, Markov Chain Monte Carlo

2010 AMS subject classifications: Primary: 62M05; Secondary: 60J22, 91G40

JEL subject classifications: C13, C63 and G32

Final author version. Article to appear in Quantitative Finance (submitted 2016 Dec 13; accepted 2017 Sep 15), doi:10.1080/14697688.2017.1383627

CRAN R-package ctmcd: Estimating the Parameters of a Continuous-Time Markov Chain from Discrete-Time Data – (https://CRAN.R-project.org/package=ctmcd

1 Introduction

Credit ratings play a key role not just in the calculation of a bank’s capital charge (amount of capital a bank must hold) but also are typically a requirement for corporations wishing to issue bonds. There are different agencies which provide firms with a rating and the credit rating the agency gives a company determines in some respect the financial health of the company. Typically ratings are of the form $AAA$, $AA$, $A$, $BBB$, $BB$, $B$, $C$, $D$ (although it varies between agencies) with, $AAA$ the best (safest), $C$ the worst (riskiest) and $D$ to imply the firm has defaulted. It is also standard for banks to use their own internal ratings system (see [YavinEtAl2014]). For an overview of the ‘science’ involved in the rating procedure see [Cantor2004].

The main object of interest in this work is the so-called annual Transition Probability Matrix (TPM), it is a stochastic matrix which shows the migration probabilities of different rated companies within a year. Rating agencies produce these annually. It is possible that such matrices are not initially stochastic due to company mergers or rounding for example. However, they can be renormalized by methods as described in [KreininSidelnikova2001] and, as argued in [BangiaDieboldKronimusEtAl2002], renormalizing non rated companies across all ratings is indeed the industry standard.

The main problem considered here is that a TPM encases transition probabilities over a $1$-year time frame and often in practice one needs a $3$ month or $10$ day transition matrix for which probabilities of default are lower than those in the TPM. Therefore one wants to accurately estimate the sub-annual matrix given the annual matrix. In the Basel proposals, Basel 3 [Basel2013]*p.3 a large part of the risk charge will be measured using ES (expected shortfall), which (as shown in [ContEtAl2010]) is extremely sensitive to a small shifts in probabilities. Therefore, accurate and consistent estimation is essential in the calculation. Moreover, with the perspective of the IFTR 9 regulation, one needs to better estimate the related probabilities of default (PD) and Markov chain’s generator since for a company whose risk profile changes significantly one needs to assess its risk throughout the bond’s lifetime. As PD are given by exponential functions, smaller initial errors will be compounded in a significant way, especially in long termed bonds. This is compounded by the known fact, corroborated by our numerical experiments, that certain algorithms overestimate the PD. Our methodology yields a way to obtain point-in-time (PIT) estimates of the Probability of Default (PD) from the through-the-cycle (TTC) estimates.

Credit rating models within the Markovian framework are handy both from a theoretical and numerical perspective. Evidence is given in [LandoSkodeberg2002] that such Markovian structure is not true in practice, nonetheless, within the Markovian structure, efficient implementation of apt Markovian credit risk models and related risk measuring estimations able to deal with massive portfolios is a challenging problem, see [BrigoMaiScherer2014, RutkowskiTarca2014]. There have been several models that produce non-Markovian effects such as mixing two generators [FrydmanSchuermann2008] or considering hidden Markov models [Korolkiewicz2012], see also [LongKeenanNeaguEtAl2011]. All non-Markovian models require in one way or another access to additional data for accurate calibration. This data is costly, needs to be updated over time and many companies opt to deal only with the TPMs. This work focuses on practitioners that do not have access to the data. The issue of rating momentum will be dealt with in forthcoming research.

The problem at hand.

We take the view of a financial agent who wishes to estimate probabilities of default or assess risk in his portfolio due to credit transitions but does not have access to (the expensive) individual credit rating transitions. The agent only has the annual TPM, say $\mathbf{P}(1)$, and uses a continuous time Markov chain (CTMC), say $(\hat{\mathbf{P}}(t))_{t\geq 0}$, with a finite state space to model the changes in rating over time. Under standard conditions the evolution of the CTMC can be written as $\hat{\mathbf{P}}(t)=e^{\mathbf{Q}t}$ where $\mathbf{Q}$ is the generator matrix. The problem is then to estimate $\mathbf{Q}$ given $\mathbf{P}(1)$. This estimation is non-trivial due to the so-called embeddability problem (not reviewed here). It is discussed in great detail by [IsraelRosenthalWei2001] and, for more of the mathematics and many of the existing results on the embeddability problem, we point the reader to [Lin2011].

Several approaches exist to tackle this estimation problem [KreininSidelnikova2001, IsraelRosenthalWei2001, TrueckOezturkmen2004, BladtSorensen2005, Inamura2006, BladtSorensen2009], either using deterministic algorithms (e.g. diagonal or weighted adjustment, Quasi-optimization of the generator) or statistical ones (Expectation-Maximization (EM), Markov chain Monte-Carlo (MCMC) ones), see Section 3. We focus on the Expectation-Maximization algorithm of [BladtSorensen2005] for CTMCs and allow for an absorbing states.

Contribution of this work.

We provide sufficient conditions to extend the convergence result of [BladtSorensen2005] to individual parameters rather than just convergence of likelihoods. The conditions presented are trivially satisfied in the context of the TPM problem. 2. 2.

We derive closed form expressions for the entries of the Hessian of the likelihood function used in the EM algorithm. This eliminated several instability issues appearing in other numerical implementations found in the literature and allows for computational speedups (comparatively). Moreover, the result provides a way to estimate the error of the estimation and assess the nature of the stationary point the algorithm has converged to. 3. 3.

We give a short overview of known methods and implement them with some modifications as to improve their performance. See Sections 3 & 4 for precise meanings: we apply the algorithms to certain credit risk problems and carry out a simulation study to check the impact in the computation of risk charges, namely IRC (Incremental Risk Charge) with VaR (Value at Risk), IDR (Incremental Default Risk) with VaR and IRC with ES (Expected Shortfall). We distinguish portfolio types (mixed, investment or speculative); the impact of different types of generators (stable vs unstable); dependence on the sample size and general convergence. We compare probabilities of default as maps of time across different algorithms and find interesting results.

For the study carried out, the implemented EM algorithm is very competitive. It is slightly slower than the deterministic algorithms but much faster than the MCMC algorithms. It embeds statistical properties like robustness that deterministic algorithms cannot capture.

Remark 1.1**.**

We focus purely on continuous over discrete time models. Continuous time algorithms yield robust estimators while the discrete ones do not, with robustness understood in the following sense: from $\mathbf{P}(1)$ estimate $\mathbf{P}(0.5)$ and $\mathbf{P}(0.25)$. From $\mathbf{P}(0.5)$ estimate $\mathbf{P}(0.25)$ again. Continuous algorithms yield the same $\mathbf{P}(0.25)$, discrete algorithms (in general) will not.

Remark 1.2** (Software availability).**

The findings and algorithms of this work are now part of an improved version of the CRAN R-package ctmcd: Estimating the Parameters of a Continuous-Time Markov Chain from Discrete-Time Data (see [Pfeuffer2017]) — https://CRAN.R-project.org/package=ctmcd

This work is organized as follows. In Section 2 we present the EM algorithm and we state our main theoretical findings. In Section 3 we briefly present other known algorithms and in Section 4 we present the benchmarking results.

2 The EM Algorithm

There exists extensive literature on the majority of the algorithms we present in this paper, therefore we only provide brief discussions and include references for additional information. Further, we will use the theory of Markov chains extensively. We do not provide details of the theory, however, interested readers can consult texts such as [Norris1998].

Preliminaries and standing convention. Throughout this manuscript we consider companies defined on a finite state space $\{1,\dots,h\}$, where each state corresponds to a rating. We denote $AAA$ as rating $1$ and $D$ (default) as rating $h$. We adopt the standard notation that $\mathbf{P}$ is an $h$-by-$h$ stochastic matrix, which will be the observed TPM (at, say, time $t=1$) and $\mathbf{Q}$ is an $h$-by-$h$ generator matrix. We further denote by $P_{ij}:=\left(\mathbf{P}\right)_{ij}$, by $q_{ij}:=\left(\mathbf{Q}\right)_{ij}$ and the intensity of state $i$ by $q_{i}=\sum_{j\neq i}q_{ij}$ where $i,j\in\{1,\dots,h\}$. A standard assumption used in credit risk modelling is that default is an absorbing state, hence $P_{hh}=1$. We work with infinitesimal generators of the following class.

Definition 2.1** (Stable-Conservative infinitesimal Generator matrix of a CTMC).**

We say a matrix $\mathbf{Q}$ is a generator matrix if the following properties are satisfied for all $i,j\in\{1,\dots,h\}$: i) $0\leq q_{ij}<\infty$ for $i\neq j$; ii) $q_{ii}\leq 0$; and iii) $\sum_{j=1}^{h}q_{ij}=0$ $\forall i$.

If a matrix $\mathbf{Q}$ satisfies the above properties, then for all $t\geq 0$ the matrix $\mathbf{P}(t):=e^{\mathbf{Q}t}$ is a stochastic matrix, where $e^{\mathbf{A}}$ is the matrix exponential of matrix $\mathbf{A}$ [Norris1998]*p.63. The goal of the algorithms presented is to calculate a generator matrix $\mathbf{Q}$ such that $e^{\mathbf{Q}t}$ is the best fit to the observed TPM, where $t$ denotes the length of time between the rating updates (typically one year).

Throughout let $(X(t))_{t\geq 0}$ denote a CTMC over the finite state space $\{1,\dots,h\}$ with a generator $\mathbf{Q}$ of the above class. Associated to $X(t)$ is, for $i,j$ in the state space, $K_{ij}(t)$ the number of jumps from $i$ to $j$ in the interval $[0,t]$ and by $S_{i}(t)$ the holding time of state $i$ in the interval $[0,t]$.

Remark 2.2**.**

If a matrix $\mathbf{P}$ is embeddable333In this setting a stochastic matrix $\mathbf{P}$ is embeddable if there exists a generator $\mathbf{Q}$ such that $\mathbf{P}=e^{\mathbf{Q}}$., the algorithms below are pointless and one can easily tackle the problem through eigenvalue decomposition etc. Or in the case where the exact timing of rating transitions are known one can use the standard maximum likelihood estimator as in [JarrowLandoTurnbull1997]. In our examples the only data given is a set of yearly TPMs which in general are not embeddable and the methods just mentioned do not yield useful results.

2.1 The Algorithm

Many methods have been developed in statistics in order to calculate the maximum likelihood estimate, but many methods break in the presence of missing data. Mathematically, we are interested in some set $\mathcal{X}$, but we are only able to observe $\mathcal{Y}$, with the assumption there is a many-to-one mapping from $\mathcal{X}$ to $\mathcal{Y}$. That is, $\mathcal{X}$ is a much richer set than $\mathcal{Y}$. When dealing with such a case, the Expectation Maximization (EM) algorithm often offers a robust solution to the problem. [McLachlanKrishnan2007] provide a complete overview of the algorithm.

The basis of algorithm is, we observe data $y$ which is a realisation (element) of $\mathcal{Y}$. We know $y$ has density function $g$ (sometimes referred to as a sampling density) depending on parameters $\Psi$ in some space $\Lambda$, but we want the density (likelihood) of $\mathcal{X}(y)$. Hence, postulate some family of densities $f$, dependent on $\Psi$, where $f$ corresponds to the density of the complete data set $\mathcal{X}(y)$ (the set of points $x\in\mathcal{X}$ which are in the pre-image of $y\in\mathcal{Y}$). The relation between $f$ and $g$ is,

[TABLE]

The idea is, the EM algorithm maximizes $g$ w.r.t. $\Psi$, but we force it to do so by using the density $f$. Further, define,

[TABLE]

where $\mathbb{E}_{\Psi}[\cdot|y]$ is the conditional expectation, conditional on $y$ under parameters $\Psi$. We assume $R$ to exist for all pairs $(\Psi^{\prime},\Psi)$, in particular we assume $f(x;\Psi)>0$ almost everywhere in $\mathcal{X}$ for all $\Psi$ (otherwise the logarithm is infinite). Let us clarify, $f$ is calculated using $\Psi^{\prime}$, but the expectation is calculated using $\Psi$. The EM algorithm is then the following iterative procedure.

Choose an initial $\Psi^{(1)}$ and take $p=1$. 2. 2.

E-step: Compute $R(\Psi;\Psi^{(p)})$. 3. 3.

M-step: Choose $\Psi^{(p+1)}$ to be the value of $\Psi\in\Lambda$ that maximizes $R(\Psi;\Psi^{(p)})$. 4. 4.

Check if the predefined convergence criteria is met, if not, take $p=p+1$ and return to (ii).

2.1.1 The particular problem of generator estimation

For our problem the observed process is a discrete time Markov chain (DTMC), the unobserved process to estimate is a continuous time Markov chain (CTMC). Therefore, the observed data is the discrete transitions and the parameters we wish to estimate are the entries in the generator. The likelihood of a continuous time fully observed Markov chain with generator $\mathbf{Q}$ is given by the following expression (see [KuchlerSorensen1997]*Chapter 3.4),

[TABLE]

with $K$ and $S$ the same as before (immediately before Remark 2.2). Hence given two generators $\mathbf{Q}^{\prime},\mathbf{Q}$, the function $R$ in (2.1) is,

[TABLE]

where $y$ denotes the discrete time observations. Maximizing for $q^{\prime}_{ij}$ in $R(\mathbf{Q}^{\prime};\mathbf{Q})$ yields

[TABLE]

The difficult step is the calculation of $\mathbb{E}_{\mathbf{Q}}[K_{ij}(t)|y]$ and $\mathbb{E}_{\mathbf{Q}}[S_{i}(t)|y]$. We follow an approach similar to [BladtSorensen2005] (see also [BladtEtAl2002]) but express the result in a framework more suited to the problem of generator estimation from TPMs, rather than the estimation from individual movements. Furthermore, the result derived in [BladtEtAl2002] is for irreducible Markov chains making it not applicable to our case (CTMC with absorbing states), accounted for in Proposition 2.4.

Consider the following functions (see [BladtEtAl2002]), for $1\leq i,j\leq h$

[TABLE]

where we denote by $\mathbf{c}=(c_{1},\cdots,c_{h})\in\mathbb{R}^{h}$ and $\mathbf{Z}\in\mathbb{R}^{h\times h}$ with $Z_{ii}=1$ for $i\in\{1,\cdots,h\}$. Observe that $V^{*}_{ij}$ is the Laplace-Stieltjes transform of the holding times $S$ and the probability generating function of the jumps $K$, with initial and final states $X(0)=i$ and $X(t)=j$ respectively. Denoting by $\mathbf{V}^{*}(\mathbf{c},\mathbf{Z};t)$ the $h$-by-$h$ matrix of elements $V^{*}_{ij}(\mathbf{c},\mathbf{Z};t)$. This allows us to give the main theorem (similar version) in [BladtEtAl2002].

Theorem 2.3**.**

For $t\geq 0$, the matrix $\mathbf{V}^{*}(\mathbf{c},\mathbf{Z};t)$ is given by,

[TABLE]

where $\bullet$ is the Schur (Hadamard) product444The Shur product of two $h\times h$ matrices $A$ and $B$ is the $h\times h$ matrix with elements $A_{ij}B_{ij}$. of matrices, $\mathbf{Q}$ is the generator matrix from the CTMC and $\Delta(\mathbf{c})$ is the diagonal matrix with entries $c_{i}$ at position $ii$ for $i=1,\dots,h$.

A closed form expression for the expectation terms in (2.3) follows from a result in [VanLoan1978] (sketched also in [HobolthJensen2011]).

Proposition 2.4**.**

Let $\textbf{e}_{i}$ be the column vector of length $h$ which is one at entry $i$ and zero elsewhere, further let us define the $2h$-by-$2h$ matrices $\mathbf{C}^{(\alpha\beta)}_{\gamma}$ and $\mathbf{C}^{(\alpha)}_{\phi}$ as,

[TABLE]

Consider a CTMC $X$ that we observe at $n$ time points $0\leq t_{1}<t_{2}<\dots<t_{n}$ and denote by $y_{s}$ the state of the Markov chain at time $t_{s}$, i.e. $y_{s}:=X(t_{s})$. Then, the expected jumps and holding times over the observations are,

[TABLE]

Proof.

Observe that $\mathbf{V}^{*}$ in (2.4) satisfies the differential equation in the statement of Theorem 2.3 (see [BladtSorensen2005]). The proof is omitted as one can solve the equation by employing the methods in [VanLoan1978]. ∎

Thus we obtain closed form expressions for the two key quantities appearing in (2.3). This approach differs from [BladtSorensen2005] where they describe numerical schemes to solve the differential equations, namely Runge-Kutta and uniformization. These techniques can yield good results at this level, but our closed form expression will pay dividends when it comes to error estimation.

This yields the relation we desire, however, in our example we have an observed TPM (or sequence of TPMs), $\mathbf{P}$, in the case of equal observation windows, $t$ in the interval $[0,T]$ (although it is trivial to generalize) the expectation can be expressed as,

[TABLE]

where $N=T/t$ (the number of observations) and $\mathbf{P}^{u}$ is the TPM of the $u$-th observation.

Remark 2.5** (The reducible case).**

Previously, we only had observed transitions, hence they must have a non-zero probability of occurring. Here we can sum $s$ and $r$ over the full range because $\mathbf{P}_{sr}(t)$ acts as an indicator of possible transitions, that is, if $P_{sr}(t)=0$ we set the $s,\,r$ component as [math]. Clearly, if $P_{sr}(t)>0$, but $(\exp\{\mathbf{Q}t\})_{sr}=0$, $\mathbf{Q}$ is misspecified.

Likelihood Convergence of the EM algorithm

In the case of this problem [BladtSorensen2005] provide a proof that the likelihood function converges with one small caveat in order to keep the parameter space compact. Namely, they use the following constrained parameter space, $\mathcal{Q}_{\epsilon}$, which can be achieved by setting, $\mathcal{Q}_{\epsilon}=\{\mathbf{Q}\in\mathcal{Q}|\text{det}[\exp(\mathbf{Q})]\geq\epsilon\}$ ($\mathcal{Q}$ is the parameter space from Definition 2.1) for some $\epsilon>0$. Theorem 4 in [BladtSorensen2005] states that the algorithm will converge to a stationary point of the likelihood or hit the boundary of the parameter space they have induced. It is accepted this is a crude approach to solving the problem and further analysis is needed when $\text{det}[\exp(\mathbf{Q})]=\epsilon$. An alternative approach would be to use a penalized likelihood as discussed in [McLachlanKrishnan2007]*p.214.

Parameter convergence criteria

The above convergence is sufficient for one to conclude convergence of the likelihood. However, it is not sufficient for convergence of the parameters as one cannot state that the series of iterates $\mathbf{Q}^{(k)}$ converge ($\|\mathbf{Q}^{(k+1)}-\mathbf{Q}^{(k)}\|\to 0$ as $k\to\infty$). From a theoretical standpoint this may not be as important as convergence of the likelihood itself, nonetheless, it is of key importance for applications. For instance, without convergence of the parameters the risk charge different financial agents obtain from the same data may vary wildly, even under very strict convergence conditions. Before proving convergence we require two important points.

Remark 2.6**.**

With (2.1.1) in mind we assume that for any $s\neq r$ such that $P^{u}_{sr}(t)=0$ for all $u$, we take the starting point $q_{sr}^{(0)}:=(\mathbf{Q}^{(0)})_{sr}=0$. As discussed in [BladtSorensen2005], any point set to zero will stay at zero for all iterations. Note, we are not changing the problem since these terms will converge to zero under the EM algorithm.

Assumption 2.7** (Element constraint).**

Similar to [BladtSorensen2005], we will use a manual space constraint to obtain the convergence. Take $1>\epsilon>0$, such that $\forall\leavevmode\nobreak\ i\neq j$, $q_{ij}<1/\epsilon$. Moreover, we assume adjacent mixing, namely, for $i\in\{2,\dots,h-1\}$, $q_{i,i\pm 1}>\epsilon$ and $q_{1,2}>\epsilon$.

We denote the space of generator matrices which satisfy this condition as $\Lambda_{\epsilon}$.

The above assumption ensures non-zero entries in the tri-diagonal band and also only finite entries as one can take $\epsilon$ as small as we wish. In the case of TPMs associated to credit ratings, such an assumption is trivially satisfied as one generally has diagonally dominant matrices and companies can always be upgraded or downgraded by one, thus $P^{u}_{i,i\pm 1}$ are typically non-zero. Diagonal dominance is sufficient for the generator to be identifiable and therefore entries do not blow up, we discuss the notion of identifiability in Section 2.2.

Proving the parameters converge is more challenging than the likelihoods, however, [Wu1983] provide a sufficient condition for this to occur, namely a sufficient condition for $\|\mathbf{Q}^{(k+1)}-\mathbf{Q}^{(k)}\|\to 0$ as $k\to\infty$ is, there exists a forcing function555A forcing function is defined as any function $\sigma:[0,\infty)\rightarrow[0,\infty)$ such that for any sequence ${t_{k}}$ defined in $[0,\infty)$, $\lim_{k\rightarrow\infty}\sigma(t_{k})=0$ implies $\lim_{k\rightarrow\infty}t_{k}=0$. $\sigma$ such that,

[TABLE]

An example of a forcing function is $\sigma(t)=\lambda t^{2}$ where $\lambda>0$. We require the following bounds on the expected values to show convergence.

Lemma 2.8**.**

Let $N$ and $\mathbf{P}^{u}$ be as defined in (2.1.1) and assume for $i\neq j$ there exists a $u\in\{1,\dots,N\}$ such that $P^{u}_{ij}>0$ (we observe a movement from $i$ to $j$ in observation window $u$). Then we obtain the following bounds on the expected number of jumps:

[TABLE]

Moreover, assuming there exists a $u\in\{1,\dots,N\}$ such that $P^{u}_{ii}>0$, we obtain the following bound on the expected holding time,

[TABLE]

To maintain the flow of the text we state immediately our main convergence result, and defer the proof of the Lemma to Appendix A.1.

Theorem 2.9**.**

Under Assumption 2.7, then, there exists a $\lambda>0$ such that for all EM iterations $k\in\mathbb{N}$,

[TABLE]

where $||\cdot||$ is the Euclidean norm.

Proof.

Writing out the difference in the $R$ terms we obtain,

[TABLE]

Due to the form of the Euclidean norm squared and the function $R$, it is sufficient to show the inequality holds for all $i\neq j$. Namely, it is sufficient to show the existence of a $\lambda>0$ such that,

[TABLE]

for all $i\neq j$. We tackle the log terms first. It is well known that we can express any $C^{\infty}$-function using Taylor expansion to a finite number of terms with some error (remainder) term. Moreover, the error term has a known form and hence, using an exact Taylor expansion to second order, there exists a $Z\in[\min(q^{(k)}_{ij},q^{(k+1)}_{ij}),\max(q^{(k)}_{ij},q^{(k+1)}_{ij})]$ such that,

[TABLE]

where we have expanded $q^{(k)}_{ij}$ around $q^{(k+1)}_{ij}$. Substituting (2.3) into the LHS of (2.1), the condition simplifies to,

[TABLE]

In order to show this bound we need to get a handle on $Z$. Clearly, there are two options between iterations, either $q^{(k)}_{ij}>q^{(k+1)}_{ij}$ or $q^{(k)}_{ij}\leq q^{(k+1)}_{ij}$. In the latter case we obtain,

[TABLE]

Since we can element wise bound $\mathbf{Q}^{(k)}$, using Lemma 2.8 and Assumption 2.7 we can bound the term $\mathbb{E}_{\mathbf{Q}^{(k)}}[K_{ij}(T)|\mathbf{P}]$ from above and $\mathbb{E}_{\mathbf{Q}^{(k)}}[S_{i}(T)|\mathbf{P}]$ from below by constants (depending on $\epsilon$). Hence, we can choose a $\lambda$ independent of $k$ such that the condition is satisfied.

The second case $q^{(k)}_{ij}>q^{(k+1)}_{ij}$, follows a similar argument. Again, we can set $Z$ as the larger of the two values, thus we obtain the following inequality,

[TABLE]

Using Lemma 2.8, we can reduce this inequality to,

[TABLE]

Since $P^{u}_{ij}>0$ and we can bound each $q_{ij}$ from above, again we choose a $\lambda$ independent of $k$. ∎

Starting value for the EM algorithm

The final point to discuss, is the choice of the initial matrix $\mathbf{Q}$. It is useful from a computational point of view to start in a good place. Here we choose $\mathbf{Q}$ based on a generalization of the QOG algorithm (described in Section 3.1) that allows for complex inputs. For each entry $q_{ij}$ we define the input as,

[TABLE]

where $|q_{ij}|$ is the magnitude of $q_{ij}$ and Re$(q_{ij})$, is the real component of $q_{ij}$. With the newly defined $\mathbf{Q}$ we apply the QOG algorithm. We take any zero entries not in the final row to be a small number ($10^{-5}$, say) unless there are zero observed transitions. This defines our initial choice of $\mathbf{Q}$. We define the EM algorithm steps as,

Take an initial intensity matrix $\mathbf{Q}$ and positive value $\epsilon$. 2. 2.

While the convergence criteria is not met and all entries of $\mathbf{Q}$ are within the boundaries

1. (1)

E-step: calculate $\mathbb{E}_{\mathbf{Q}}[K_{ij}(T)|\mathbf{P}]$ and $\mathbb{E}_{\mathbf{Q}}[S_{i}(T)|\mathbf{P}]$. 2. (2)

M-step: set $q^{\prime}_{ij}=\mathbb{E}_{\mathbf{Q}}[K_{ij}(T)|\mathbf{P}]/\mathbb{E}_{\mathbf{Q}}[S_{i}(T)|\mathbf{P}]$, for all $i\neq j$ and set $q_{ii}$ appropriately. 3. (3)

Set $\mathbf{Q}=\mathbf{Q}^{\prime}$ (where $\mathbf{Q}^{\prime}$ is the matrix of $q^{\prime}$s) and return to E-step. 3. 3.

End while and return $\mathbf{Q}$.

This leads to the following theorem for convergence of the EM.

Theorem 2.10** (Convergence of the EM).**

Assume that our initial point is in the parameter space $\Lambda_{\epsilon}$: is a true generator and satisfies Assumption 2.7. Then either the sequence of points $\{\mathbf{Q}^{(k)}\}_{k}$ converges to a single point in $\Lambda_{\epsilon}$ which is also a stationary point of the likelihood, or the entries go to the boundary (blow up or some tri-diagonal elements in an non-absorbing row go to zero).

A proof of Theorem 2.10 follows directly from Theorem 4 in [BladtSorensen2005] and our Theorem 2.9.

Remark 2.11** (The unique maximizer of the Likelihood).**

The natural question one may ask is does this stronger form of convergence imply convergence to the global maximum? The problem of existence and uniqueness of maximum likelihoods in this setting is a very challenging problem with a long history. [BladtSorensen2005] give a wonderful overview on the subject, Theorem 1 in [BladtSorensen2005] also provides results on existence and uniqueness of the maximum. Unfortunately, one cannot say more than this, if one can derive conditions under which a unique maximum existed (for non-embeddable TPMs) then the above convergence result is sufficient to conclude the EM will converge to the MLE.

During the final revision of this manuscript, [PfeufferMoestelFischer2017] came to our attention. There, the authors propose two algorithms to mitigate the effect of the EM’s possible convergence towards a local but not necessarily the global maximum of the likelihood function (no proofs are given). Our Theorem 2.9 is handy in this context as it shows that once the EM lands “near” the global maximum the iteration will converge to it.

2.2 Variance Estimation

In this section we derive an expression for the Hessian of the likelihood. We use a result in [Oakes1999] and follow [BladtSorensen2009], however, unlike [BladtSorensen2009], we provide a closed form expression for the Hessian. This result eliminates the stability problems observed in the numerical simulation case when the entries in $\mathbf{Q}$ are small. The Hessian provides a way to estimate the standard error of the maximum likelihood estimates and further allows us to assess the nature of the converged stationary point (this is further discussed in Section 4.4).

We point the reader to [BladtSorensen2005]*Theorem 1 for results on the existence and uniqueness of maximum likelihood estimators with respect to this problem. Further, for discussions on consistency and asymptotic normality related to this problem one should consult [KremerWeissbach2013], [KremerWeissbach2014]. [KremerWeissbach2013], provide sufficient conditions for consistency, the key assumption relies on so-called model identifiability666In our setting a model is identifiable if there are no two generators $\mathbf{Q}\neq\mathbf{Q}^{\prime}$ such that $\exp\{\mathbf{Q}t\}=\exp\{\mathbf{Q}^{\prime}t\}$.. [KremerWeissbach2013] prove identifiability under conditions which are too restrictive for our purpose; [BladtSorensen2005, DehayYao2007] discusses the problem of identifiability in detail. From [Cuthbert1973], [BladtSorensen2005] for the model to be identifiable it is sufficient (though very crude) to have $\min_{i}(\exp\{\mathbf{Q}t\})_{ii}>1/2$, [Culver1966] gives a requirement for general matrices based on the eigenvalues which one can always aposteriori verify after a $Q$ is deduced. The crucial assumption in [KremerWeissbach2014] to obtain asymptotic normality, is that the Hessian must be invertible at the true value, we can of course verify invertibility a posteriori.

We recall a result from [Oakes1999] for calculating the Hessian of the likelihood.

Lemma 2.12**.**

The second derivative of the likelihood with parameter $\Psi$ and observed information $y$ is related to the EM function $R$ by

[TABLE]

Injecting (2.2) in the above we obtain,

[TABLE]

[TABLE]

where $\delta_{ab}$ is the Kronecker delta. From our previous work, (2.10) is easy to obtain, however, (2.11) involves derivatives of the expected jumps and holding times and is thus challenging. [BladtSorensen2009] opt for a simple numerical scheme to compute these derivatives and found unstable results, although the authors do remark that more sophisticated numerical schemes could yield improved results at greater computational expense.

It is worth noting we have made no comment on the allowed values of $\alpha,\beta,\mu$ and $\nu$, other than the clear restriction that they belong to $\{1,\dots,h\}$. Let us now state the following definition.

Definition 2.13** (Allowed pairs).**

We say that the pair $\alpha,\beta$ is allowed if $\alpha\neq\beta$ (not in the diagonal) and $q_{\alpha\beta}$ is not converging to zero under the EM algorithm.

For practical applications, one can imagine the set of allowed values, as the set of $\alpha,\beta$ such that $q_{\alpha\beta}>\epsilon$, where $\epsilon$ is some cut-off point ($10^{-8}$, say). The reason we must exclude small parameters is, this analysis only holds in the large data limit, since we do not have an infinite amount of data we cannot for certain rule out some jump, however, if $q_{\alpha\beta}$ is converging to zero, it implies that this parameter is either zero or extremely close to zero and therefore we can bound it above by a small number. Moreover, from a mathematical point of view a parameter which does tend to zero (or even becomes zero) lies on the boundary, where the notion of differentiability is not clear. Therefore, we can think of the “allowed pairs” as the variables when solving the problem in a restricted parameter space. We now present the following theorem.

Theorem 2.14**.**

Let $\mu,\nu,\alpha,\beta\in\{1,\dots,h\}$, and $\mathbf{Q}$, $\mathbf{Q}^{\prime}$ be two generator matrices ($\in\Lambda_{\epsilon}$ for some $\epsilon>0$). For any two allowed pairs $\alpha,\beta$ and $\mu,\nu$, the derivative in (2.11) is,

[TABLE]

where the $2h$-by-$2h$ matrices, $\mathbf{C}^{(\alpha\beta)}_{\gamma},\mathbf{C}^{(\alpha)}_{\phi},\mathbf{C}^{(\alpha\beta)}_{\eta}$, are defined as,

[TABLE]

and the $4h$-by-$4h$ matrices $\mathbf{C}_{\psi}^{(\alpha\beta,\mu\nu)},\mathbf{C}_{\omega}^{(\alpha\beta,\mu)}$ are defined

[TABLE]

The proof of this uses similar techniques to Proposition 2.4 along with differentiation properties of matrix-exponentials, and is deferred to Appendix A.2.

Remark 2.15**.**

In the above representation for the derivative of $R$, we use subscripts of the form $h+y_{s+1}$ and $3h+y_{s+1}$, this is simply a consequence of the result in [VanLoan1978]. Namely, we are not interested in all the entries of the matrix, only an $h$-by-$h$ segment. We therefore need to adjust the indexing to only take elements at this specific segment.

Using Theorem 2.14 and Lemma 2.12, we can write the elements of the Hessian corresponding to the $q_{\alpha\beta}q_{\mu\nu}$ derivative as,

[TABLE]

A similar transform to (2.1.1) can be applied here to obtain the Hessian from TPMs. When using the result to estimate the error, some knowledge of the number of companies per rating is required.

2.2.1 Computation of the Error

Due to the Hessian only being defined for parameters $q_{\alpha\beta}>0$, some parameters are not included in the Hessian, thus the matrix is smaller than $(h-1)^{2}$-by-$(h-1)^{2}$. We compute the Hessian as follows,

• •

Let $N_{a}$ be the number of allowed points in the estimated $\mathbf{Q}$ returned from the EM algorithm.

• •

Define an $N_{a}$-by-$2$ matrix $\mathbf{V}_{\mathbf{Q}}$ as the matrix which records the allowed (in the sense of Definition 2.13) components of $\mathbf{Q}$. The $ij^{th}$ component of the Hessian is then the differential,

[TABLE]

• •

If we denote the Hessian by the $N_{a}$-by-$N_{a}$ matrix $\mathbf{H}(\cdot)$, then the information matrix is given by $-\mathbf{H}(\cdot)$. The estimated variance of the allowed parameter $q_{ab}$ is then the $i^{th}$ diagonal element of $-\mathbf{H}(\cdot)^{-1}$, where $\mathbf{V}_{\mathbf{Q}}(i,1)=a$ and $\mathbf{V}_{\mathbf{Q}}(i,2)=b$.

• •

Following standard statistics, the normal based $95\%$ confidence interval for $q_{ab}$ is $q_{ab}\pm 1.96\sqrt{Var(q_{ab})}$.

3 Competitor Algorithms

3.1 Deterministic algorithms

Diagonal Adjustment (DA)

The first method to discuss is diagonal adjustment, see [IsraelRosenthalWei2001]. Given a TPM, $\mathbf{P}$, one calculates the matrix logarithm directly. However, due to the embeddability problem, the logarithm may not be a valid generator. To solve this problem [IsraelRosenthalWei2001] suggest setting for $i\neq j$,

[TABLE]

and adjusting (re-balancing) the diagonal element correspondingly, $q^{DA}_{ii}=\sum_{j\neq i}-q_{ij}$ for $i\in\{1,\cdots,h\}$. Hence forcing the corresponding matrix $\mathbf{Q}^{DA}$ to satisfy the properties of a generator.

Weighted Adjustment (WA)

Weighted adjustment is also suggested in [IsraelRosenthalWei2001]. It follows diagonal adjustment except, one re-balances across the entire row. Again, calculate the logarithm of the TPM to find $q$’s, then compute

[TABLE]

The entries corresponding to weighted adjustment are defined as,

[TABLE]

Quasi-Optimisation of the Generator (QOG)

The above two methods are unfortunately not optimal in any sense. The QOG (Quasi-Optimisation of the Generator), method suggested in [KreininSidelnikova2001] relies on optimisation and is therefore an improvement on the diagonal and weighted adjustment methods. QOG involve solves the minimization problem $\min_{\mathbf{Q}\in\mathcal{Q}}\|\mathbf{Q}-\log(\mathbf{P})\|$, where $\mathcal{Q}$ is the space of stable generator matrices and $||\cdot||$ is the Euclidean norm. Further, [KreininSidelnikova2001] provide an efficient algorithm to obtain $\mathbf{Q}$.

3.2 Statistical algorithm: Markov Chain Monte Carlo

An alternative statistical algorithm one can adopt is MCMC (Markov Chain Monte Carlo). For the reader’s convenience we have included a summary of MCMC in Appendix B. It should be noted that all MCMC algorithms presented here use a so-called auxiliary variable technique, by introducing the fully observed Markov chain, $X$ as a random variable. Moreover, the prior for Q is $\Gamma(\alpha,1/\beta)$ (shape and scale), which is conjugate for the likelihood of a CTMC.

3.2.1 Gibbs Sampling - Bladt & Sorensen 2005

To simulate the Markov process, $X$, [BladtSorensen2005] suggest a rejection sampling method. As is stated in [BladtSorensen2005], such a sampling method runs into difficulties when considering low probability events since the rejection rate will be high (e.g. default for high rated bonds). The MCMC algorithm is summarised as follows, [Inamura2006],

Construct an initial generator $\mathbf{Q}$ using the prior distribution ($\Gamma(\alpha_{ij},1/\beta_{i})$). 2. 2.

For some specified number of runs

1. (1)

Simulate $X$ for each observation from $Y$, with law according to $\mathbf{Q}$. 2. (2)

Calculate the quantities of interest $K$ and $S$, from $X$. 3. (3)

Construct a new $\mathbf{Q}$ by drawing samples from $\Gamma(K_{ij}(t)+\alpha_{ij},1/(S_{i}(t)+\beta_{i}))$. 4. (4)

Save this $\mathbf{Q}$ and use it in the next simulation. 3. 3.

From the list of $\mathbf{Q}$s, drop some proportion (burn in), then take the mean of the remainder.

The issues with this method are the choice of $\alpha$ and $\beta$ and the number of runs required before we know that the sample has converged (burn in). Both of these are critical in obtaining accurate answers from MCMC and although [BladtSorensen2005] suggested taking $\alpha_{ij}$ and $\beta_{i}$ to be 1, they observe MCMC overestimating entries in the generator when true entries were small. Furthermore, here we are required to use the TPM indirectly through inferring company transitions. That is, we consider $M$ companies in each rating and define the number of companies to make the transition $i$ to $j$ as $M\times P_{ij}$, this of course need not be an integer, but we can always normalise the entries. The reason we cannot use the TPM directly as we did in the EM algorithm is due to the fact that MCMC becomes very sensitive to the values in the prior. The burn in for MCMC will be of little concern to us here as will become apparent when carrying out analysis on the algorithms.

3.2.2 Importance Sampling - Bladt & Sorensen 2009

[BladtSorensen2009] address some of the issues in [BladtSorensen2005] by running the same algorithm as previous combined with an importance sampling scheme based on the Metropolis-Hastings algorithm (in its essence a single component Metropolis-Hastings algorithm). The proposal distribution suggested is a Markov chain with generator given by the ‘neutral matrix’ $\mathbf{Q}^{*}$, which takes the following form,

[TABLE]

where $\mathbf{1}_{h}$ and $\textbf{I}_{h}$ is the $h$-by-$h$ matrix of ones and identity matrix respectively and $W$ is a scaling factor set to match the intensities in the true generator matrix $\mathbf{Q}$. [BladtSorensen2009] note that if entries in $\mathbf{Q}$ are known to be zero, then the corresponding element in $\mathbf{Q}^{*}$ should also be set to zero and the diagonal modified accordingly. Thus transitions rarely produced by the generated Markov chain will occur much more frequently under $\mathbf{Q}^{*}$. Thus we have solved (at least partially) one of the problems faced in MCMC. The importance sampling weights for a chain $X$ are,

[TABLE]

where $L$ is the CTMC likelihood. For the priors, [BladtSorensen2009] do not suggest any significant improvement on their earlier work. The authors use $\alpha=1$ and $\beta=5$, which they claim gives better results than the suggestion in [BladtSorensen2005]. However, it still provides a problem when dealing with entries in $\mathbf{Q}$ which are close to zero. The problem stems from the fact that very little information is known (rarely observed) for certain transitions, therefore the output for these entries is mostly based on our prior beliefs.

3.2.3 MCMC Mode Algorithm

[Inamura2006] presented an alternative algorithm to the original MCMC algorithm presented in [BladtSorensen2005] whereby one calculates the mode rather than the mean. The author claims that this gives extremely accurate results and outperforms other algorithms. The reasoning presented is that the standard MCMC overestimates in the small probability cases due to the gamma distribution being ‘skewed’, therefore the mode is a better estimate. [Inamura2006] approximates the mode of $\{q_{ij}^{(k)}\}$ by kernel smoothing over the estimates (after taking the log transform to ensure all results are positive).

Remark 3.1** (Other MCMC based estimators).**

Many extensions and different MCMC methods to solve this problem are possible (e.g. priors as hyperparameters or sequential Monte Carlo techniques). Here, we consider less complex MCMC algorithms which already set the tone for a comparative study.

4 Benchmarking the Algorithms

Due to the diversity of investments bank’s make, one cannot assess an algorithms’ performance with a single test. With this in mind we consider a host of tests on different portfolios and matrices. The computations were carried out on a Dell PowerEdge R430 with four Intel Xeon E5-2680 processors. During the review process of our work we found [Pfeuffer2017] with an R-implementation of some of the algorithms covered in the previous section. The performance tests of [Pfeuffer2017] are a just subset of those we present next and independently confirm (where there is overlap) our findings, in particular the timing of the MCMC algorithms versus the EM. A version of our algorithms will appear in the mentioned R-package (see Remark 1.2).

The first observation we make is, transition matrices can vary substantially depending on the financial climate (see [ChristensenEtAl2004] and [Cantor2004]). Therefore we consider two different generator matrices which can be thought of as the generator in financial stress and the generator in financial calm. In order to keep these matrices ‘reasonable’ we start off with the generator given in [ChristensenEtAl2004] built using a large amount of data (see also [Inamura2006]) and consider a generator which has in general higher transition rates and one with lower transition rates. Through considering more than one generator this provides a more detailed assessment of the performance of the various algorithms than other comparative reviews, such as [Inamura2006]. The generators we consider are shown in Table 4.1 and Table 4.2. We observe that Table 4.1, has more non-zero entries and larger entries than that of Table 4.2.

Throughout the analysis we refer to the multiple MCMC algorithms introduced in Section 3 which we label in the following way: MCMC BS05 is [BladtSorensen2005]’s algorithm of Section 3.2.1; MCMC BS09 is [BladtSorensen2009]’s algorithm of Section 3.2.2; and MCMC Mode is [Inamura2006]’s algorithm in Section 3.2.3.

4.1 Sample Size Inference

The first test we consider is an extension to a test in [Inamura2006], where the author considers a true underlying generator and masks it by using it to simulate TPMs, which we view as observations, then applying the algorithms to each observation. The key point here is, [Inamura2006] only simulates $100$ companies per rating and hence the outputted TPM is non-embeddable (has [math] entries for accessible jumps). This is an extremely useful test because it provides a fair and intuitive way to assess the performance of each algorithm, however, [Inamura2006] only considers one true generator and only one level of information i.e. $100$ companies per rating. Alongside the two different generators we also consider a range of companies per rating to determine its effect on convergence for each algorithm. Furthermore, [Inamura2006] uses seven years worth of data, although one would likely have access to multiple years worth of TPM data, it is unlikely that we would have seven years of transitions from the same generator. Hence we consider four years, which is more consistent with time homogeneity estimates for generators (see [ChristensenEtAl2004]). We calculate our estimates for the generator as follows.

Take a range of obligors per rating, $[100,200,300,500,750,1000]$ and $10$ random seeds. 2. 2.

For each true generator simulate four one year TPMs for each seed and for each obligor per rating. Hence we have (#Years$\times$#Obligors categories$\times$#Random Seeds$\times$#True generators), simulated TPMs. 3. 3.

For each set of four simulated TPM we estimate the generator for each algorithm. MCMC may take a long time to run, therefore we consider the time taken to carry out the first 10 runs and the total time taken, if these exceed 180 or 18 000 seconds respectively, the algorithm is deemed to be too slow and no result is returned. Note, MCMC algorithms use 3000 runs with a burn in of 300. This is smaller than [Inamura2006] for example, however, [Inamura2006] shows apparent convergence to the stationary distribution in a small number of iterations and we observe a similar result. 4. 4.

Therefore, for each algorithm we have (# Obligors categories $\times$ # Random Seeds $\times$ # True generators) estimated generators to analyze.

We analyze the estimated generators by considering, distance between estimated generator and true generator in Euclidean norm and difference in one year probability of default. All results presented have been obtained by analyzing the estimated generator for each seed, then averaging. This gives a better picture of the average performance.

4.1.1 Convergence in Euclidean Norm

Our goal in this analysis is to consider the empirical rate of improvement of each algorithm as our ‘information’ about the true generator increases. For each obligor category we calculate the natural log of the distance (measured by the Euclidean norm) between the estimate and the true. The results are shown in figures  4.1  and  4.2.

Note the $x$-axis is on a logarithmic scale. We observe similarities between the two figures, most notably in the case of low information all algorithms have very similar convergence results, however, as we increase the information there is substantial variation in improvement, MCMC BS09 algorithm does not improve as well as the other algorithms. Missing points stem from an algorithm failing the acceptance times.

The MCMC algorithms have a potentially increased error due to the Monte Carlo simulation, lowering it requires a larger computational expense to the already most expensive algorithm being tested here. For the [BladtSorensen2009] algorithm, the neutral matrix approximation may give poor mixing, thus the additional error.

4.1.2 Error in Probability of Default

Although overall error is important, it does not provide details on the small probability scale. This is extremely important in banking, since estimation of the probability of default is crucial. Using the same estimated generators as previous we calculate the corresponding one year TPM, that is, we calculate $\exp\{\mathbf{Q}_{\text{estimate}}\}$ (using the expm function in MATLAB) for each seed then take the average. The averaged TPM default probabilities are compared to the true ones. To keep the numbers in the comparisons meaningful we plot the log of the relative error, where we define,

[TABLE]

The results of which are given in Figures 4.3 and 4.4.

Unlike the overall error, there appears to be far greater volatility in the error estimation w.r.t. the probability of default. Moreover, there appears to be no general downward trend in error for the investment grade ratings. A likely cause for this is, even with 1000 companies there are still no/few investment grade defaults. Of the algorithms MCMC BS09 performs the worst. The EM algorithm though has consistently one of the smallest errors and is clearly the best in the investment grades.We have only shown the results for the unstable generator, the stable generator was similar.

4.2 Time Dependent Probability of Default

A key question that has not been addressed in the literature is how do the probabilities of default change in time among the several algorithms. For ths we only consider EM, QOG, WA and the MCMC Mode algorithm from [Inamura2006], since these algorithms gave the best probability of default estimates.

We consider a non-embeddable TPM, then estimate the generator matrix $\mathbf{Q}$, from $\mathbf{Q}$ we can easily calculate the probability of a company with some initial rating defaulting in time $t>0$. The goal here is to assess how that probability changes with time. The TPM is given in Table 4.4, for the MCMC algorithm we took this table to be generated with 250 obligors per rating.

The probability of default across ratings over the one year time horizon is found in Figure 4.5. The plots give a deeper understanding to the algorithms themselves. As the probability of default increases the algorithms converge, however, in the case of less defaults we observe a much larger discrepancy. This can be thought of as the algorithm’s ability to deal with missing data, in the lower grades we observe defaults and thus have a handle on the probability, however, in the case of AAA ratings we observe no defaults and therefore it is an approximation by the algorithm. This shows the difference between the methods, shows the potential prior dependence in the MCMC algorithm. What is also extremely interesting is that QOG set the jump in the generator from AA to C as zero (even though the TPM has a non zero entry there), this implies QOG may in some places under estimate the risk for the investment ratings, this can be seen by the fact QOG puts a smaller probability of default on AAA.

There is a clear ovestimation of the probability of default at higher grades by the WA and MCMC algorithms.

4.3 Risk Charge

The previous tests have been rather theoretical, we now consider a practical test to asses the performance of these algorithms in calculating risk charges. We do not give much discussion to the calculation of these risk charges for more technical details readers should consult texts such as [SkoglundChen2011]. Here we consider multiple stylized portfolios to represent the risk appetites of different banks. To best of our knowledge analysis into how different risk measures react to different portfolio types has not been considered in the literature. The risk charges we consider are IRC (VaR at 99.9% with a 3 months liquidity horizon including mark to market loss), IDR (VaR at 99.9% over one year only considering default) and a theoretical risk charge which is IRC but measured using Expected Shortfall (ES) at 97.5%. The final risk charge is included due to the Basel committee showing an increasing interest in ES. We consider 4 years worth of simulated data, and to keep the analysis realistic we consider 200 companies per rating. We consider 3 different portfolios corresponding to risk adverse (all investment grade), a speculative portfolio (all speculative grades) and finally a mixed portfolio. The portfolios considered are given in Tables 4.7, 4.7 and 4.7. The tables show the values and ratings of the various bonds in each portfolio.

Alongside these portfolios we calculate the risk charges using the following information,

• •

The interest rates we receive for a bond in each rating are

[TABLE]

These figures are based on interest rates from Moody’s and can be found in Section 4.1 of [SkoglundChen2011]. Although these interest rates do not technically match the generators we are using for the TPMs they provide reasonable interest rates for our toy example.

• •

We assume that all money is lost in the case of default (zero recovery rate).

• •

We calculate credit migration using the one factor777This is technically not the true regulation for the calculation of IDR which requires a two factor model, however our goal here is only to use these calculations as a method for comparing algorithms. credit metrics model ([GuptonEtAl1997]), i.e. normalised asset returns follow,

[TABLE]

where $X$ is the systematic risk, $\epsilon_{i}$ is the idiosyncratic risk both standard normally distributed and $\beta_{i}$ is the correlation to the systematic risk, defined in [Basel2003]*p.50,

[TABLE]

where $P^{D}_{i}$ is the probability of default of asset $i$. Consequently we see that the higher $P_{i}^{D}$ the lower the value of $\beta$.

• •

Although more sophisticated methods are available for calculation of VaR and ES (see [Fermanian2014]), we calculate the risk charges using Monte Carlo. This is sufficient here since the portfolios are small relative to a typical bank portfolio, therefore we can obtain accurate estimates using a reasonable number of simulations.

• •

Again, we calculate 10 realizations of the TPMs and estimate a generator for each.

We consider 15$\times 10^{5}$ simulations for each portfolio, to assess whether this was sufficient we calculated VaR and ES using 7.5$\times 10^{5}$, 10$\times 10^{5}$, 12.5$\times 10^{5}$ and 15$\times 10^{5}$ simulations and found the difference between 7.5$\times 10^{5}$ and 15$\times 10^{5}$ to be $<5\%$ for all cases. Hence were are confident that 15$\times 10^{5}$ gives sufficiently accurate results for our purposes.

With respect to the risk charge calculation, similar to the previous analysis, we calculate the risk charges for every set of TPMs, then average over all the seeds to obtain the risk charge. The risk charges as set by the true generators are given in Table 4.8.

To asses the performance of each algorithm we measure the error by the following,

[TABLE]

where $\text{Risk Charge Estimate}(i)$ is the $i^{th}$ realization of the risk charge and $N$ is the number of TPM sets (10 here). The results obtained by the algorithms are shown in Table 4.9.

It should be noted, in the stable IDR some algorithms produce a non-zero value for the investment portfolio, therefore we have inserted the money value. The first observation we make is, all algorithms overestimate the risk for the investment portfolio. This is down to two key feature, one is the ‘step like’ nature of VaR, where in a small portfolio, small probability changes can make a large difference. The other is because we are averaging over multiple Monte Carlo simulations, thus having one default in one of those realizations will change the overall average dramatically. In terms of a typical bank portfolio this type of error should not be a problem since we would be dealing with a far larger number of assets and hence one would obtain multiple defaults. However, the results do still give a useful comparison between the algorithms. Although the MCMC algorithms can outperform the deterministic algorithms for the speculative grades, remarkably in all categories the EM produces the best results. From the tests we have considered we conclude the EM to be the superior algorithm for this problem.

4.4 Error estimation of the EM algorithm

A major advantage of the statistical algorithms over their deterministic counterparts is that one can derive error estimates (confidence intervals) without the brute force (slightly ad-hoc) method of bootstrapping. For MCMC this comes by looking at the posterior distribution, which we get for free. However, as we have seen MCMC is computationally expensive and we have derived a relatively cheap formula to calculate the confidence intervals from the EM algorithm.

In a similar fashion to the analysis we have carried out previously we now test the error estimate given by the EM. Again, we mask the true generator by using simulated TPMs, however, here we only consider the scenario of 300 obligors per rating, but the number of years worth of data is varied. That is, we simulate 50 years worth of TPMs and then apply the EM algorithm using 1 years worth then 2 years etc up to 50 years. This analysis shows both the estimated error for the parameter and also how the error changes when more information is added. It should also be noted that we replace companies who have defaulted to the rating they were pre default. This keeps the number of companies in the system constant and can be thought of as the flow of new companies being rated. Moreover, this is only one realisation of the data, hence the parameter estimate and confidence intervals are not particularly smooth.

The transitions shown in Figures 4.6 and 4.7 were chosen to show a spectrum of the magnitudes in the generators, the other entries not shown are similar. The first point to make is that the true value of the parameter is almost always within the confidence interval and the confidence interval shrinks as the number of years increases. One of most important features though is that the confidence interval is only small when the EM is stable and close to the true parameter, hence the EM is not “over confident” but is also providing reasonable estimates on its own error. The final point to make is, although some confidence intervals go below zero by a small amount, this is only true in the case where the parameter is extremely close to zero initially and further, once more data is considered, all parameters have confidence interval which are strictly positive.

Remark 4.1** (Confidence intervals in practice and regulation).**

Figures 4.6 and 4.7 show a very important feature. Namely, how much the estimate can vary with data, especially when the parameter is reasonably small. Such analysis exposes the variance in the estimate and how much data is actually required before the estimate levels out. From this example though the confidence intervals calculated from the information matrix appear to be able to capture this error. In the view of future regulation it may be prudent to take such confidence intervals into account when considering risk charges.

Connection to the Global Maximum

A previous problem with the EM was one could not be sure of the nature of the stationary point. However, we know the form of the Hessian, and therefore we can easily check if this point is a maximum by assessing the eigenvalues of this matrix. Clearly, if we were not at a maximum, then it would be worth perturbing the outputted generator and rerunning the algorithm. As discussed in Remark 2.11 the question of a global maximum is very difficult in this setting.

Remark 4.2**.**

One way that has been suggested to improve the chances of the EM converging to the global maximum is, to start from multiple points. Here we can consider creating starting points by setting for each $i\neq j$, $q_{ij}\sim\text{Exp}(\lambda)$ for an appropriate $\lambda$ then setting $q_{ii}$ appropriately.

We tested the EM according to the above remark and found in every case considered the EM always returns the same generator.

5 Conclusions and future research

In this manuscript we built upon the closed form expressions for the expected number of jumps and holding times of a CTMC with an absorbing state, over given observations and used the results to derive a closed form expression for the Hessian of the likelihood. This coupled with stronger convergence has elevated the EM algorithm to be the optimal algorithm to tackle this problem.

Across the battery of tests carried out, the EM algorithm outperforms competing algorithms. The EM is a tractable algorithm, slower than the deterministic algorithms but still several orders of magnitude faster than the Markov-Chain Monte-Carlo alternatives (Table 4.2). The statistical algorithms (EM and MCMC) embed a strong robustness property for the estimator contrary to the deterministic algorithms, i.e. the likelihood is far less sensitive to small changes in the underlying TPM. In terms of estimating risk charges, the EM algorithm has superior results in all scenarios.

On the more practical side, Figure 4.5 highlights that for lower ratings algorithms produce essentially the same estimates for the probabilities of default while a palpable difference emerges at higher ratings. Moreover, the error estimates in the EM may provide a sensible way to test in effect model risk.

Lastly, non-Markovian phenomena like rating momentum (see [LandoSkodeberg2002]) and appropriate models to tackle it will be addressed in forthcoming research.

Acknowledgement(s)

The authors would like to thank Dr. R. P. Jena at Nomura Bank plc London for the helpful comments. In addition, the authors would like to thank Ruth King (U. of Edinburgh), Ioannis Papastathopoulos (U. of Edinburgh) and Samuel Cohen (Oxford Uni.) for the helpful discussions.

We thank as well the two anonymous referees for their comments which led to improvements on the initial submission.

Funding

G. Smith was supported by The Maxwell Institute Graduate School in Analysis and its Applications, a Centre for Doctoral Training funded by the UK Engineering and Physical Sciences Research Council (grant [EP/L016508/01]), the Scottish Funding Council, Heriot-Watt University and the University of Edinburgh.

G. dos Reis acknowledges support from the Fundaç$\tilde{\text{a}}$o para a Ci$\hat{e}$ncia e a Tecnologia (Portuguese Foundation for Science and Technology) through the project [UID/MAT/00297/2013] (Centro de Matemática e Aplicaç$\tilde{\text{o}}$es CMA/FCT/UNL).

Appendix A Proofs

A.1 Proof of Lemma 2.8

We now provide the proof of Lemma 2.8, all terms used have the same definition as they did when the Lemma was stated. Throughout we assume $i\neq h$, thus from from Assumption 2.7 $\mathbb{P}_{\mathbf{Q}}(X(t)=j|X(0)=i)>0$ for all $j\in\{1,\dots,h\}$ and $t>0$. The first inequality we prove is the lower bound on the expected number of jumps. Following the assumptions in Lemma 2.8 and time homogeneity we make the observation

[TABLE]

The above inequality holds because we are only considering $X(0)=i$, $X(t)=j$ and not all possible combinations of start and end states, moreover, $\mathbb{P}_{\mathbf{Q}}(K_{ij}\geq 1|X(0)=i,X(t)=j)\leq\sum_{n=1}^{\infty}n\mathbb{P}_{\mathbf{Q}}(K_{ij}=n|X(0)=i,X(t)=j)$. We further observe,

[TABLE]

Thus the lower bound in inequality (2.7) can be easily obtained. We now prove the upper bound on the expected number of jumps. The first observation we make is for all $\nu\in\{1,\dots,h\}$,

[TABLE]

To see this, let $\mu\neq i$, then denote by $\tau_{i}$ the first time the process enters state $i$ (if $\mathbb{P}_{\mathbf{Q}}(X(t)=i|X(0)=\mu)=0$ for $t>0$, then the result is trivial), by the law of total probability we find,

[TABLE]

The second term is zero. Then, using the Markov property we obtain,

[TABLE]

Consequently from this observation and (2.1.1) we obtain,

[TABLE]

Observe that,

[TABLE]

The numerator is easy to bound by considering the expected number of jumps out of $i$,

[TABLE]

The denominator requires further analysis, firstly, let $n=|i-\nu|$, and therefore by Assumption 2.7 we can go from state $i$ to $\nu$ in $n$ jumps, w.l.o.g. let $i\geq\nu$ (it will be come clear that the ordering does not matter). Firstly, if $i=\nu$ then,

[TABLE]

For $i>\nu$, we use the Markov property to obtain,

[TABLE]

Conditioning on $X$ only making one jump in each increment we obtain,

[TABLE]

As $n\leq h$ and the terms are strictly smaller than $1$, the sought result follows (independent of $\nu\neq i$).

The last inequality to prove concerns the holding times. By taking for $P^{u}_{ii}>0$,

[TABLE]

where the final inequality follows by simply considering the case of no jumps. We can then apply the bounds from Assumption 2.7 to complete the inequality.

A.2 Proof of Theorem 2.14

We recall from [Wilcox1967], [TsaiChan2003] that for a square matrix $\mathbf{M}$ whose elements depend on a vector of parameters $\{\lambda_{1},\dots\lambda_{r}\}$ (for $r\in\mathbb{N}$), the following identity holds

[TABLE]

for all $i\in\{1,\dots,r\}$. Let $\mu,\nu,\alpha,\beta\in\{1,\dots,h\}$. Recalling Proposition 2.4, differentiating $\mathbb{E}_{\mathbf{Q}}[K_{\mu\nu}(t)|y]$ w.r.t. $q_{\alpha\beta}$ yields,

[TABLE]

Note that although the expected value of $K$ only depends on individual elements of the matrix and not the full matrix, we are still able to use the differentiation result since $A_{ij}=\mathbf{e}_{i}^{\intercal}\mathbf{A}\mathbf{e}_{j}$. Hence, from (A.1) we obtain,

[TABLE]

Clearly, since $q_{\alpha\beta}$ appears twice in $\mathbf{Q}$,

[TABLE]

Then, by [VanLoan1978] we can solve these integrals explicitly to obtain,

[TABLE]

again $\mathbf{C}^{(\alpha\beta)}_{\eta}$ and $\mathbf{C}_{\psi}^{(\alpha\beta,\mu\nu)}$ are as defined in the Theorem’s statement.

Therefore, we have a closed form expression for the derivative of expected jumps w.r.t. $q_{\alpha\beta}$. Applying a similar argument for the expected holding time we obtain,

[TABLE]

where $\mathbf{C}_{\omega}^{(\alpha\beta,\mu)}$ is as defined in the Theorem. Combining these yields the required result.

Appendix B Overview of Markov Chain Monte Carlo algorithm

For details on the Markov Chain Monte Carlo (MCMC) theory we refer the reader to [GilksRichardsonEtAl1996]. Algorithms for implementing MCMC to estimate a generator, from discrete observations are discussed in [BladtSorensen2005] and [BladtSorensen2009]. MCMC differs from the EM in the sense that EM estimates the set of parameters which maximises the likelihood function, while MCMC samples from the posterior distribution. Namely, given some data $D$, the posterior distribution of parameters $\theta$ is $\pi(\theta|D)$, which by Bayes’ theorem is,

[TABLE]

with $\pi(D|\theta)$ denoting the likelihood and $\pi(\theta)$ the prior distribution. MCMC obtains the best guess of $\theta$ by sampling from $\pi(\theta|D)$ and taking the Monte Carlo approximation of the expected value. The reason the expectation is our best guess is due to the fact we use both the data (likelihood) but also our experience on what the outcome should approximately be (the prior). Although the prior can be extremely useful in stopping ‘bad’ answers it is also a criticism of MCMC due to so-called prior sensitivity.

Remark B.1**.**

Here we purely discuss MCMC to sample from the posterior, algorithms which approximate the maximum likelihood in the presence of missing data do exist, but are more useful when for example one cannot explicitly write the E step in the EM algorithm (see [GelfandCarlin1993]).

Similar to the case of the EM algorithm the problem faced here is missing data. Namely we wish to consider the so-called posterior distribution of the generator matrix $\mathbf{Q}$, which we denote by $\pi(\mathbf{Q}|D)$ (although it is common to suppress the data and only write $\pi(\mathbf{Q})$). The difficulty is, in its current state this is an extremely hard distribution to evaluate so we augment with an auxiliary variable $X$ (see [GilksRichardsonEtAl1996]*p.105 and [BesagGreen1993]). In general $X$ need not require an interpretation, although here it will correspond to the full Markov chain. In order to generate realisations of $\pi(\mathbf{Q}|D)$, we specify the conditional distribution $\pi(X|\mathbf{Q},D)$ which provides the joint distribution $\pi(\mathbf{Q},X|D)=\pi(\mathbf{Q}|D)\pi(X|\mathbf{Q},D)$ and therefore the marginal distribution of $\mathbf{Q}$ is $\pi(\mathbf{Q}|D)$. One can then sample from the marginal distribution by using any sampling method that preserves the joint distribution $\pi(\mathbf{Q},X|D)$ (and by extension $\pi(\mathbf{Q}|D)$), such as Gibbs or Metropolis Hastings.

The method used in [BladtSorensen2005] and [BladtSorensen2009] is the data augmentation algorithm from [TannerWong1987] (see also [LittleRubin2002]*p.200). We specify the prior distribution $\pi(\mathbf{Q})$ and take a realisation from this distribution, $\mathbf{Q}^{(0)}$, we then construct a sequence $\{\mathbf{Q}^{(k)},X^{(k)}\}$, for $k=1,\dots,M$ by:

Draw, $X^{(k)}\sim\pi(X|\mathbf{Q}^{(k-1)},D)$. 2. 2.

Draw, $\mathbf{Q}^{(k)}\sim\pi(\mathbf{Q}|X^{(k)},D)=\pi(\mathbf{Q}|X^{(k)})$ (since $X^{(k)}$ is richer than $D$). 3. 3.

Save $\{\mathbf{Q}^{(k)},X^{(k)}\}$ and take $k=k+1$.

Under mild conditions (see [GilksRichardsonEtAl1996]*Chapter 4), after some burn-in $n$, the sequence $\{\mathbf{Q}^{(k)},X^{(k)}\}$ for $k\geq n$ has the same distribution as $\pi(\mathbf{Q},X|D)$. Moreover, the marginals also have the correct distribution, namely, $\{\mathbf{Q}^{(k)}\}\sim\pi(\mathbf{Q}|D)$ for $k\geq n$. Therefore we estimate the generator matrix by, $\frac{1}{M-n+1}\sum_{k=n}^{M}\mathbf{Q}^{(k)}$.

For the choice of prior, $\pi(\mathbf{Q})$, [BladtSorensen2005] suggest a prior from the gamma distribution with shape $\alpha_{ij}$ and scale $1/\beta_{i}$. Hence, $q_{ij}\sim\Gamma(\alpha_{ij},1/\beta_{i})$, where $\alpha_{ij},\beta_{i}\geq 0$, $\forall\,i\neq j\in\{1,\dots,h\}$. With this choice, the prior is a conjugate prior. Although this prior has some drawbacks, we note, by assuming the prior to follow a Gamma distribution we effectively bound the parameter space, therefore there is no need to make the space compact. Noting that the posterior distribution of $X$ is equivalent to the likelihood i.e. $\pi(X|\mathbf{Q})=L_{t}(X;\mathbf{Q})$, one has

[TABLE]

From the likelihood of a CTMC and the assumption on the prior we infer that,

[TABLE]

We do not have equality here since there is no normalisation term. We generate $q_{ij}$ with $i\neq j$ from the distribution \Gamma\big{(}K_{ij}(t)+\alpha_{ij},1/(S_{i}(t)+\beta_{i})\big{)} (since each $q_{ij}$ is independent).