# Credit Portfolios, Credibility Theory, and Dynamic Empirical Bayes

Authors: Tze Leung Lai
Journal: ISRN Probability and Statistics
 (2012)

---

## Abstract

We begin with a review of (a) the pricing theory of multiname credit derivatives to hedge the credit risk of a portfolio of corporate bonds and (b) current approaches to modeling correlated default intensities. We then consider pricing of insurance contracts using credibility theory in actuarial science. After a brief discussion of the similarities and differences of both pricing theories, we propose a new unified approach, which uses recent advances in dynamic empirical Bayes modeling, to evolutionary credibility in insurance rate-making and default modeling of credit portfolios.


---

## Body

## 1. Introduction

Credit markets provide liquidity to the underlying financial system. They include not only banks that give loans to consumers and businesses, but also financial institutions and corporations that develop and trade securities and derivative contracts derived from these loans and credit-sensitive investments such as corporate bonds, and insurance companies and their subsidiaries that issue insurance and reinsurance contracts to protect against adverse credit events. The recent financial crisis in the USA was associated with unpreparedly high default rates of mortgage loans and defaultable bonds in credit markets, beginning with subprime mortgage loans in 2007 and culminating in the collapse of large financial institutions such as Bear Stearns and Lehman Brothers in 2008. In Sections2 and 3, we give a brief review of the conventional approaches to modeling default risk and of widely used derivative contracts to hedge credit risk, respectively. In particular, we point out in Section 2 the similarities and differences between the intensity models in credit risk and those in survival analysis.Parallel to the increasing volume of subprime mortgage loans whose value was estimated to be $1.3 trillion by March 2007, an important development in financial markets during the past decade was the rapid growth of credit derivatives, culminating in $32 trillion worth of notional principal for outstanding credit derivatives by December 2009. These derivative contracts are used to hedge against credit loss of either a single one or a portfolio of corporate bonds. A “credit default swap” (CDS) is a contract between a buyer, who obtains the right to sell bonds issued by a company or country (a “single name” that is called thereference entity), and a seller (of the CDS contract) who agrees to buy the bonds for their face values when a credit event occurs. The face value of a coupon-bearing bond is the principal that the issuer repays at maturity if it does not default. The total face value of the bonds is the notional principal of the CDS. The buyer of the CDS makes periodic payments (usually quarterly) to the seller until the expiration date T of the CDS or until a credit event occurs. If there is no credit event up to time T, the buyer receives no payoff. If there is a credit event at time τ<T, the buyer's payoff is the difference between the face value and the mid-market value of the cheapest deliverable bond determined by an auction process several days after the credit event if the contract specifies cash settlement, or between the face value and what the buyer receives by selling the bond if the contract specifies physical settlement. The credit event is defined in the CDS contract and includes failure to pay according to schedule, agency downgrade, bankruptcy, and debt restructuring. The credit risk of a portfolio of corporate bonds can be mitigated by using a multiname credit derivative. Stochastic models of default intensities are used to price these derivatives. For multiname credit derivatives, joint modeling of the default intensities has been a challenging statistical problem. The review in Section 3.3 describes the correlated default models that have been widely used prior to the financial crisis and points out some of their limitations and deficiencies.In July 2007, Bear Stearns disclosed that two of its subprime hedge funds, the High-Grade Structured Credit Fund and the High-Grade Structured Credit Strategies Enhanced Leverage Master Fund, which were invested in CDOs, had lost nearly all their value following a rapid decline in the subprime mortgage market. Lawsuits were brought by investors against former managers of the hedge funds, and Standard & Poor's (S&P) downgraded the company's credit rating. In March 2008, the Federal Reserve Bank of New York initially agreed to provide a $25 billion collateralized 28-day loan to Bear Stearns, but subsequently changed the deal to make a $30 billion loan to J. P. Morgan Chase to purchase Bear Stearns at $2 a share, or less than 7% of Bear Stearns' market value just two days before. After a class action suit was filed on behalf of the shareholders, J. P. Morgan Chase raised the offer to $10 per share, which was approved by the shareholders in May 2008. Lehman Brothers also suffered unprecedented loss for its large positions in subprime and other lower-rated mortgage-backed securities in 2008. After attempts to sell it to Korea Development Bank and then to Bank of America and to Barclays failed, it filed for Chapter 11 bankruptcy protection on September 15, 2008, making the largest bankruptcy filing, with over $600 billion in assets, in US history.A day after Lehman's collapse, American International Group (AIG) needed bailout by the Federal Reserve Bank, which gave the insurance company a secured credit facility of up to $85 billion to enable it to meet collateral obligations after its credit ratings were downgraded below AA, in exchange for a stock warrant for 79.9% of its equity. AIG's London unit had sold credit protection in the form of CDS and CDO to insure $44 billion worth of securities originally rated AAA. As Lehman's stock price was plummeting, investors found that AIG had valued its subprime mortgage-backed securities at 1.7 to 2 times the values used by Lehman and lost confidence in AIG. Its share prices had fallen over 95% to just $1.25 by September 16, 2008, from a 52-week high of $70.13. The “contagion” phenomenon, from increased default probabilities of subprime mortgages to those of counterparties in credit derivative contracts whose values vary with credit ratings, was mostly neglected in the models of joint default intensities that were used to price CDOs and mortgage-backed securities. These models also failed to predict well the “frailty” traits of latent macroeconomic variables that underlie the mortgages and mortgage-backed securities. Section4 reviews recent works in the literature to incorporate frailty and contagion in modeling joint default intensities.The major business of AIG is insurance, and insurance contracts rather than credit derivatives contracts are what it sells primarily. These insurance contracts provide protection against loss due to automobile accidents and other unexpected events such as fire and floods in home insurance, just like credit derivatives that provide protection against loss due to default or other credit events. On the other hand, as will be explained in Section5, whereas pricing credit derivatives involves the risk-neutral stochastic dynamics of the credit events that are priced by the market, pricing insurance contracts involves real-world (rather than risk-neutral) statistical models of the events. Section 5 discusses the similarities and differences between both pricing theories and also reviews credibility theory in actuarial science, which aims at deriving the premium of an insurance contract that balances the experience of an individual risk with the class risk experience. In particular, pricing a CDO that provides CDS protection for a portfolio of companies and setting the premium of automobile insurance for drivers in the same risk class both involve predicting the probability of occurrence and the loss incurred for certain future events. Estimation of the probability and the loss can pool information from the individual (company or driver) and the collective (portfolio of companies or risk group) levels.Such pooling of information is the essence of empirical Bayes (EB) methodolgy in statistics, introduced by Robbins [1] and Stein [2]. The review of credibility theory in Section 5 also reviews EB and points out their connections. Lai et al. [3] have recently shown how the predictive performance of conventional EB methods in estimating the accident-proneness of a driver in a future period, or other insurance applications such as workers' compensation losses due to permanent partial disability, can be improved by using a dynamic EB approach to longitudinal data. Section 6 describes this new approach and relates it to generalized linear mixed models (GLMM) in biostatistics introduced by Breslow and Clayton [4]. It also reviews recent work applying the dynamic EB approach via GLMM to credibility theory, credit portfolios, and default probability forecasts of subprime mortgage loans. Thus the dynamic EB methodology provides a unified statistical approach to credit portfolios and credibility theory, and Section 7 provides further discussion and concluding remarks.


## 2. Intensity and Structural Models of Default Risk

### 2.1. Survival Analysis, Intensity Processes, and Intensity Models for Default Risk

We begin with a brief review of survival analysis in biostatistics and reliability theory and then show its connections with widely used intensity models for pricing corporate bonds and other credit-sensitive securities. To pinpoint these connections, our review focuses on the martingale approach to survival analysis, as described in [5, 6]. A major theme of the statistical problems in survival analysis is to study the distribution of a failure time τ, based on possibly censored data from either a homogeneous population or a regression model with covariates. The survival function of τ is S(t)=P(τ≥t)=1-F(t-), where F(s)=P(τ≤s) and F(t-)=lims→t-F(s)=P(τ<t). Clearly, S(t) is a nonincreasing left-continuous function of t with S(0)=1 and limt→∞S(t)=0. If τ is absolutely continuous, the probability density function of τ is f(t)=-dS(t)/dt, and the hazard (or intensity) function is
(2.1)λ(t)=limh→0+P(t≤τ<t+h∣τ≥t)h=f(t)S(t).
Integrating (2.1) with respect to t yields
(2.2)S(t)=exp[-∫0t‍λ(u)du]=exp[-Λ(t)],
where Λ(t)=∫0t‍λ(s)ds is called the cumulative hazard function. Without assuming F to be continuous, we can define the cumulative hazard function by the Stieltjes integral
(2.3)Λ(t)=∫0t‍dF(s)1-F(s-).
Regarding the indicator functions I{τ≤t},t≥0, as a stochastic process, an important property of Λ is that
(2.4){I{τ≤t}-Λ(τ∧t),t≥0}isamartingale,
where τ∧t=min(τ,t). Fleming and Harrington [5] and Andersen et al. [6] give overviews of continuous-time martingales and their applications, via (2.4) and its analogs, to survival analysis. Based on n independent observations τi sampled from F, the nonparametric maximum likelihood estimate of F is the empirical distribution function F^(t)=(1/n)∑i=1n‍I{τi≤t}. However, in applications of survival analysis to biomedical and econometric data, the τi may not be fully observable because of censoring. Some subjects (or firms) may not fail during the observation period or may have been lost in followup during the period; the data on these individuals are said to be right-censored.Right-censoring can be formulated by introducing censoring variablesci that indicate the time of loss to followup. The censoring variable ci can also indicate the length of the observation period. The observations, therefore, are (Ti,δi),i=1,…,n, where Ti=min(τi,ci) and δi=I{τi≤ci} is the censoring indicator that indicates whether Ti is an actual failure time or is censored. Subject i is “at risk” at time s if Ti≥s (i.e., has not failed and has not been lost to followup prior to s). Let
(2.5)Y(s)=∑i=1n‍I{Ti≥s},N(s)=∑i=1n‍I{Ti≤s,δi=1}.
Note that Y(s) is the risk set size and ΔN(s)=N(s)-N(s-) is the number of observed deaths at time s, and that ΔN(s)/Y(s) is the analog of p^j. We use ΔN(s) to denote the jump size of N at s (i.e., ΔN(s)=N(s)-N(s-)) and also use the convention 0/0=0. The cumulative hazard function Λ(t) can be estimated by
(2.6)Λ^(t)=∑s≤t‍ΔN(s)Y(s)=∫0t‍I{Y(s)>0}Y(s)dN(s)
for the censored data (Ti,δi),1≤i≤n. This is called the Nelson-Aalen estimator [7, 8]. Analogous to (2.4),
(2.7){∫0tU(s)[dN(s)-Y(s)dΛ(s)],t≥0}isamartingale,
for every left-continuous stochastic process U(s); note that Y(s) is left-continuous. From this and the martingale central limit theorem,
(2.8)Λ^(t)-Λ(t){∫0t‍(I{Y(s)>0}/Y2(s))dN(s)}1/2→𝒟N(0,1)
as n→∞, where →𝒟 denotes convergence in distribution; see [5, 6].Altman [9] and others have developed mortality tables for loans and bonds by using methods in actuarial science. These mortality tables (also called life tables) are used by actuaries to set premiums for life insurance policies. Partition time into disjoint intervals I1=(0,t1],I2=(t1,t2], and so forth. A life table summarizes the mortality results of a large cohort of n subjects as follows:nj= number of subjects alive at the beginning of Ij,dj= number of deaths during Ij, lj= number lost to followup during Ij.It estimates pj=P(dyingduringIj∣aliveatthebeginningofIj) by p^j=dj/(nj-lj) so that the actuarial (life-table) estimate of P(τ>tk) is the product
(2.9)∏j=1k‍(1-p^j)=∏j=1k‍(1-djnj-lj).
Without discretizing the failure times, Kaplan and Meier [10] introduced the product-limit estimator of S(t)=P(τ≥t):
(2.10)S^(t)=∏s<t‍(1-ΔN(s)Y(s)).
Since N(s) has at most n jumps, the product in (2.10) is finite. Note that S^(t)=∏s<t‍(1-ΔΛ^(s)) by (2.6) and (2.10). As shown in [6],
(2.11)S(t)=∏s<t‍(1-dΛ(s))∀tsuchthatΛ(t)<∞.
The product in (2.11) is called the product integral of the nondecreasing function Λ(·) on [0,t]; it is defined by the limit of ∏{1-[Λ(ti)-Λ(ti-1)]}, where 0=t0<t1<⋯<tm=t is a partition of [0,t] and the product is in the natural order from left to right, as the mesh size max1≤i≤m|ti-ti-1| approaches 0. Moreover, by making use of the martingale central limit theorem, it is shown in [5, 6] that
(2.12)S^(t)-S(t)[S^(t){∫0t‍(I{Y(s)>0}/Y2(s))dN(s)}1/2]→𝒟N(0,1).


#### 2.1.1. Regression Models for Hazard Functions with Covariates

In applications, one often wants to use a model forτ to predict future failures from a vector x(t) of predictors based on current and past observations; x(t) is called a time-varying covariate. When x(t)=x does not depend on t, it is called a time-independent (or baseline) covariate. In practice, some predictors may be time-independent, while other components of x(t) may be time-varying. Since prediction of future default from x(t) is relevant only if τ>t (i.e., if default has not occurred at or before t), one is interested in modeling the conditional distribution of τ given τ>t, for example, by relating the hazard function λ(t) to x(t). Noting that λ(t)dt is the probability that τ∈[t,t+dt], Cox [11] introduced the proportional hazards (or Cox regression) model
(2.13)λ(t)=λ0(t)exp(βTx(t)).
Putting x(t)=0 in (2.13) shows that λ0(·) is also a hazard function; it is called the baseline hazard. Cox regression fits (2.13) to the observed data (Ti,δi,xi(t)),t≤Ti,1≤i≤n. Lane et al. [12] apply Cox regression to predict bank default, using time-independent covariates.The parameters of (2.13) are β and the cumulative hazard function Λ0, which is an infinite-dimensional parameter. One way to estimate (β,Λ0) is to specify Λ0 via a parameter vector θ so that (β,θ) can be estimated by maximum likelihood. Since an observed failure contributes a factor fi(Ti) to the likelihood, where fi=λiSi is the density function of the τi while a censored failure contributes a factor Si(Ti) to the likelihood, and since logSi(t)=-Λi(t) and Λi(t)=Λ0(t)exp(βTxi(t)), the log-likelihood function can be written as
(2.14)l(β,θ)=∑i=1n{δilogfi(Ti)+(1-δi)logSi(Ti)}=∑i=1n{δi[logλ0(Ti;θ)+βTxi(Ti)]-Λ0(Ti;θ)eβTxi(Ti)};
see [13]. Commonly used parametric models for failure times are the exponential density θe-θt that has constant hazard θ>0 and the Weibull density θtθ-1exp(-tθ) that has hazard function θtθ-1 and cumulative hazard function tθ, with θ>0. Lee and Urrutia [14] use a Weibull distribution of default times to predict insolvency in the property-liability insurance industry. Hazard models with time-varying covariates have been used by McDonald and Van de Gucht [15] to model the timing of high-yield bond default and of exercising call options. Related survival analysis to predict bankruptcy has been carried out in [16–19].Instead of assuming a parametric model to estimate the baseline hazard functionλ0(·), Cox [11, 20] introduced a semiparametric method to estimate the finite-dimensional parameter β in the presence of an infinite-dimensional nuisance parameter λ0(·); it is semiparametric in the sense of being nonparametric in λ0 but parametric in β (in terms of βTxt). Cox's partial likelihood method decomposes the likelihood function into two factors, with one involving only β and the other involving both β and Λ. It estimates β by maximizing the following partial likelihood, which is a factor that only involves β. Order the observed censored failure times as τ(1)<⋯<τ(m), with m≤n. Let Cj denote the set of censored Ti's in the interval [τ(j-1),τ(j)), and let (j) denote the individual failing at τ(j), noting that with probability 1 there is only one failure at τ(j) because the failure time distributions have density functions. Let R(j)={i:Ti≥τ(j)} denote the risk set at τ(j). Cox's regression estimator β^ is the maximizer of the partial log-likelihood
(2.15)l(β)=∑j=1m{βTx(j)(τ(j))-log(∑i∈R(j)exp(βTx(i)(τ(j))))},
or equivalently, is the solution of (∂/∂β)l(β)=0. Martingale theory can be used to show that
(2.16)(-l¨(β))1/2(β^-β0)hasalimitingstandardnormaldistributionasn→∞,
where l¨(β) is the Hessian matrix of second partial derivatives (∂2/∂βk∂βh)l(β). One, therefore, can perform usual likelihood inference, treating the partial likelihood as a likelihood function. Moreover, even though β^ is based on partial likelihood, it has been shown to be asymptotically efficient; see [6]. Modifying Y(s) in (2.5) to
(2.17)Y(s)=∑i∈R(j)exp(βTxi(s))ats=τ(j)
leads to the Breslow estimator [21] of Λ0, which is again given by (2.6) but with Y(s) defined by (2.17).


#### 2.1.2. Intensity-Based Models for Pricing Default Risk

Besides the important roles they have played in the statistical analysis of failure-time data as explained in the preceding section, intensity modeling has also become a standard approach to pricing corporate bonds and other credit-sensitive securities since the 1990s. Intensity models provide an elegant way to combine term-structure modeling for default-free bonds with the default risk and recovery features of corporate bonds. Chapter 10 of [22] gives an overview of interest rate markets and the term structure of interest rates, statistical methods to estimate the zero-coupon yield curve from the current prices of default-free Treasury bonds, arbitrage-free pricing theory of interest rate derivatives, short rate and other models under the risk-neutral (instead of the real-world) measure entailed by arbitrage-free pricing, and calibration of these models to market data involving actively traded interest rate derivatives. For corporate bonds, besides the interest rate processes, we also include an additional stochastic process for the default intensity of the bond. All the stochastic processes are defined on the same probability space in which the probability measure Q is a risk-neutral (or martingale) measure associated with arbitrage-free pricing. Expectations taken under this measure will be denoted by E^ or EQ, while E is used to denote expectation under the real-world (or physical) measure P.ACox process (also called stochastic Poisson process) is often used to model the default intensity of the bond's issuer. The default intensity λt is assumed to be governed by an exogenous stochastic process Xt,t≥0, so that λt=λ(Xt) and the stochastic dynamics of λt is specified through Xt. Let τ be the default time and Λ(t)=∫0t‍λ(Xs)ds. Then Λ(τ) is distributed as an exponential random variable ϵ1 with mean 1 that is independent of {Xs,s≥0} since P{Λ(τ)>t} = P{τ>Λ-1(t)} = S(Λ-1(t)) = exp{-Λ(Λ-1(t))} = e-t. Hence we can use λ(Xs) to generate τ by τ=inf{t:∫0t‍λ(Xs)ds≥ϵ1}. Consider a zero-coupon bond, with maturity date T and face value 1, issued by a firm at time 0. Assume that there is a short-rate process r(Xs) under the risk-neutral measure Q such that the default-free bond price is given by
(2.18)p(0,T)=E^{exp(-∫0Tr(Xs)ds)};
see [22, page 257]. Assuming λ(Xs) to be the intensity process for the default time τ of the firm and assuming zero recovery at default, the price of the defaultable bond at time 0 is
(2.19)π(0,T)=E^{I{τ>T}exp(-∫0Tr(Xs)ds)}=E^{E^[⋯∣Xs,s≤T]}=E^{exp(-∫0Tr(Xs)ds)E^[I{τ>T}∣Xs,s≤T]}=E^{exp(-∫0T(r+λ)(Xs)ds)},
since E^[I{τ>T}∣Xs,s≤T]=Q{∫0T‍λ(Xs)<ϵ1∣Xs,s≤T}=exp(-∫0T‍λ(Xs)ds), noting that P(ϵ1>x)=e-x and that ϵ1 is independent of {Xs,s≥0}. Thus (2.19) replaces the short rate in the default-free bond pricing formula by the default-adjusted short rate (r+λ)(Xs). Moreover, closed-form expressions are available for π(0,T) if one uses affine models for both r and λ; see [23, 24]. The short rate process and the intensity process have unknown parameters which can be encoded in a vector θ and can be estimated from the bond prices by using least squares in the same way as that for Treasury bonds; see [22]. Note that unlike survival analysis in biostatistics that uses censored failure-time data to estimate the intensity in (2.13) or its integral in (2.3), estimation of the parameters of the intensity process is based not on observed defaults but rather on the bond prices that involve the risk-neutral (rather than the real-world) measure.In practice, recovery rates also need to be considered and the vectorXs may include predictors that make r(Xs) and λ(Xs) much more complex than the affine model that gives an explicit formula for the integral in (2.19). Monte Carlo methods are used instead to evaluate (2.19) in these more complex situations. Since corporate bonds typically promise cash flows, which mimic those of Treasury bonds, it is natural to consider yields when one compares the effects of different recovery assumptions at default. The yield at date t of a zero-coupon bond with price Bt, face value A, and maturity T is defined by
(2.20)y(t,T)=1T-tlogABt,
and the difference between the yield of a defaultable bond and a corresponding Treasury bond is referred to as the credit spread or yield spread. The first type of recovery assumptions is recovery of market value. This measures the change in market prices before and after the default, representing the loss in the bond's value associated with the default. More generally, consider the price process of a defaultable contingent claim V promising a payoff of f(XT) at T. The claim is said to have a fractional recovery of market value of δ at default time τ if the amount recovered in the event of default is equal to
(2.21)h(τ)=δV(τ-)forτ≤T,
where V(τ-) is the value of the claim just prior to default and δ∈[0,1). With this recovery assumption, the price at t of the defaultable contingent claim (if there has not been a default up to date t) is
(2.22)V(t)=E^[f(XT)exp(-∫tT(r+(1-δ)λ)(Xs)ds)∣Xs,s≤T],
which is an extension of (2.19). A variant of recovery of market value is recovery of face value, which is the closest to legal practice in the sense that debt with the same priority is assigned a fractional recovery depending on the outstanding notional amount but not on maturity or coupon. It is also the measure typically used in rating-agency studies. The third type of recovery assumptions is recovery of treasury, under which the corporate bond in default is replaced with a treasury bond with the same maturity but a reduced payment. Unlike recovery of face value, it tries to correct for the fact that amounts of principal with long maturity should be discounted more than principal payments with short maturity.


### 2.2. Econometric Modeling of Default Risk

The intensity-based models in the preceding section are sometimes calledreduced-form models, using a terminology in econometrics that refers to models of the relationship of an economic variable to exogenous variables. In contrast, structural models attempt to derive the relationship from some underlying economic theory. Structural models for the default of a firm assume that default is triggered when the firm's asset and debt values enter some region whose boundary is called a “default barrier.”


#### 2.2.1. Merton's Model of Debt and Equity Values and the KMV Credit Monitor Model

As in the Black-Scholes model described in Section 8.1 of [22], Merton [25] assumes a market with continuous trading, no transaction costs, unlimited short selling and perfectly divisible assets, and with a risk-free asset that has constant interest rate r. The asset value At of the firm is assumed to follow a geometric Brownian motion (GBM) dAt/At=μdt+σdwt. Suppose at time 0 the firm has issued two kinds of contingent claims: equity and debt, in which debt is assumed to be a zero-coupon bond with face value D and maturity date T. We can regard the firm as being run by equity owners. At maturity date T of the bond, equity owners pay D to retain ownership if At≥D, and would not pay D if At<D since they have limited liability. In the latter case, bond holders receive a recovery of AT instead of the promised payment D. Therefore, the equity value at time t is given by the Black-Scholes formula
(2.23)Ve(t)=EQ{e-r(T-t)(AT-D)+∣At}=AtΦ(d1)-De-r(T-t)Φ(d2),
where Q is the risk-neutral measure under which the drift of the GBM is r instead of μ, and
(2.24)d1=log(At/D)+(r+σ2/2)(T-t)σT-t,d2=d1-σT-t;
see [22, page 186]. Similarly, the debt value at time t is
(2.25)Vd(t)=EQ{e-r(T-t)min(AT,D)∣At}=At-Ve(t),
since min(At,D)=At-(At-D)+.In the 1990s, the KMV Corporation of San Francisco (now merged into Moody's KMV) introduced its Credit Monitor Model for default prediction using Merton's framework. To begin with, the GBM for the asset value can be expressed in terms of the Brownian motionwt by At=A0exp{(μ-(1/2)σ2)t+σwt}; see [22, page 66]. Therefore, the probability of default on the debt at maturity is
(2.26)P(AT<D)=P{(μ-σ22)T+σwT<log(DA0)}=Φ(-logA0+(μ-σ2/2)T-logDσT)=Φ(-E(logAt)-logDσT).
This leads to KMV's distance to default defined by
(2.27)DD={E(logAT)-logD}σT
so that the predicted probability of default under Merton's model is Φ(-DD). Since σT is the standard deviation of logAT, DD can be regarded as the number of standard deviations of the log asset value from the log liability. Instead of using the normal distribution for logAT to predict the default probability as in (2.26), the KMV Credit Monitor Model uses a large proprietary database of firms and firm defaults to construct an EDF (expected default frequency) score for a firm with κ standard deviations away from the default (i.e., DD =κ) at the beginning of a future period. The idea is to stratify the firms based on their DD's and estimate empirically the default probability of the firm from its associated stratum. The firm's DD is evaluated from its balance sheet and option models for its debt.


#### 2.2.2. Black-Cox Framework of Default Barriers and Other Extensions of Merton's Model

Black and Cox [26] extend the Merton model by allowing defaults to occur prior to the maturity of the bond; specifically, when the asset value At hits a nonrandom, continuous boundary g(t)≤D. They also allow the equity owners to be paid dividends, at rate q, so that under the risk-neutral measure Q the GBM for At has drift r-q; see [22, page 186]. Let τ=inf{0≤t≤T:At=g(t)}. The payoff to bond holders is
(2.28)min(AT,D)I{τ>T}+g(τ)I{τ≤T},
in which the first summand reflects that the bond holders receive min(AT,D) at the bond's maturity if default has not occurred up to that time, while the second summand reflects that the bond holders receive the recovery of Aτ=g(τ) upon default. Similarly, the payoff to equity holders is (AT-D)+I{τ>T}.Stochastic models for the default-free interest rate dynamics (see [22, Chapter 10]) have been used to replace the constant r in Merton's model. Another extension is to allow coupon payments to the bond holders. More general jump diffusion models have been introduced to replace GBM in Merton's model. Lando [24] provides an overview of these extensions, which also have been considered in conjunction with default barriers.


### 2.3. Credit Ratings and Credit Value at Risk

Credit scoring has been used in the banking industry since the 1940s to assess the creditworthiness of loan applicants. The basic idea is to assign scores to loan applications based on their attributes such as gender, marital status, education level, income, profession, and credit history. Such scoring systems, calledscorecards, are usually developed by in-house credit departments and/or consulting firms. One of the most widely used credit scoring systems was developed by Fair Isaac & Company (FICO), which built credit scoring systems for US banks and retail stores in the 1960s and 1970s and then expanded those systems to meet the needs of other industries and to evaluate the credit of borrowers. Credit scoring systems have now been developed for virtually all types of credit analysis, from consumer credit to commercial loans. The idea is to preidentify certain key factors that determine the probability of default and combine or weight them into a quantitative score.Credit ratings of the creditworthiness of corporate bonds are provided by rating agencies such as Moody's, S&P, and Fitch. Using Moody's system, the ratings are Aaa, Aa, A, Baa, Ba, B, Caa, Ca, and C and “investment grade” ratings are Baa or above. The corresponding S&P and Fitch ratings are AAA, AA, A, BBB, BB, B, CCC, CC, and C. For finer ratings, Moody's divides the Aa rating category into Aa1, Aa2, Aa3, the A category into A1, A2, A3, and so forth, and S&P and Fitch divide AA into AA+, AA, and AA−, A into A+, A, and A−, and so forth. Although methodologies and standards differ from one rating agency to another, regulators generally do not make distinctions among the agencies. Furthermore, although there is a strong similarity between the rating systems of Moody's and S&P, different agencies might assign slightly different ratings for the same bond. The ratings given by rating agencies are used by several credit risk software models such as CreditMetrics of J. P. Morgan, a system that evaluates risks individually or across an entire portfolio. Most corporations approach rating agencies to request a rating prior to sale or registration of a debt issue. S&P assigns and publishes ratings for all public corporate debt issues over $50 million, with or without a request from the issuer; but in all instances, S&P analytical staff first contact the issuer to call for meetings.


#### 2.3.1. CreditMetrics, Credit Ratings Migration, and Product-Limit Estimators

CreditMetrics was introduced in 1997 by J. P. Morgan and its cosponsors (including KMV, UBS, Bank of America, and others). The system measures the risk and potential losses of loans and portfolios of bonds and loans over a period (e.g., a year) by simulation and analytic methods. In particular, it uses Monte Carlo simulations to estimate credit rating changes of all counterparties. Historical data provided by rating agencies on credit rating changes can be used to estimate the transition probability matrix for credit rating changes. In the Basel Accord for bank regulation, banks are required to construct these transition probability matrices based on their own data to stress test their portfolios. Some risk management systems such as CreditMetrics use these credit rating transition matrices for Monte Carlo simulations of credit losses due to credit ratings migration in a future period. A Markov chain model of ratings migration assumes a finite state space𝒮={1,…,K} consisting of different rating classes, in which state 1 denotes the best credit rating class and state K represents the default case. The K×K probability transition matrix of the Markov chain is
(2.29)P=(p11p12⋯p1Kp21p22⋯p2K⋮⋮⋮⋮pK-1,1pK-1,2⋯pK-1,K00⋯1),
in which pii=1-∑j≠i‍pij for all i and pij represents the actual probability of moving to state j from initial rating state i in one time step. Note that the default state K is an absorbing state. The estimates of credit transition matrices published by rating agencies use this Markov chain framework.Suppose there areNi firms in a given rating category i at the beginning of the period and that Nij firms migrate to the category j at the end of the period. An obvious estimate of pij is p^ij=Nij/Ni for j≠i. However, there are serious difficulties with this estimate in practice. First, if no transition from rating class i to j is observed in the period, then the estimated transition probability p^ij becomes 0. Second, the assumption of homogeneity over time is too strong for credit transition matrices. The product-limit estimator (2.10) for the survival function has been generalized by Aalen [27], Aalen and Johansen [28], and Fleming [29] to estimate the transition probability matrix P(s,t) of a nonhomogeneous continuous-time Markov process on a finite state-space. Analogous to (2.11), the transition probability matrix from time s to t can be expressed as the product integral P(s,t)=∏(s,t]‍(I+dA(u)), and the components Aij of the cumulative transition intensity matrix can be estimated by the Nelson-Aalen estimator of the form (2.6); see Section  IV.4 of [6].To simulate the distribution of credit losses due to credit rating changes of two counterparties, say one with initial Aaa rating (state 1) and the other with initial B rating, over a 1-year period, we can use the estimateP^ of the transition matrix (2.29) to sample the new ratings for both counterparties in the following way. Without assuming the rating changes of the two counterparties to be independent, a Gaussian copula model, which will be described in Section 3.3, is typically used to simulate the joint probability distribution of the rating changes. Therefore, we first generate two standard normal random variables ZA and ZB with correlation ρ. Noting that Φ(ZA) is uniform, the new rating of the Aaa company remains Aaa if Φ(ZA)<p^11, moves to states 2 (i.e., Aa) if p^11≤Φ(ZA)<p^11+p^12, and so forth. This is tantamount to generating the state transitions of the Aaa company by the first row of the probability transition matrix P^. While the state transitions of the B company are likewise generated by the corresponding row of P^, the state transitions of the two companies are not independent and their dependence structure is inherited from the correlation of ZA and ZB. From the simulated credit rating changes of the two counterparties we can revalue their outstanding contracts to determine the total credit losses during the year. From a large number of such simulations, a probability distribution for credit losses is obtained, which can be used to evaluate credit VaR (Values at Risk) for risk management.


#### 2.3.2. Credit VaR and a One-Factor Default Model

A100(1-α)% credit VaR over a period of T years is the (1-α)th quantile of the credit loss within the period. J. P. Morgan's CreditMetrics defines credit losses as those arising not only from defaults but also from credit ratings downgrade. This involves estimating the probability distribution of credit losses by Monte Carlo simulation of defaults and credit rating changes of all counterparties. Instead of Monte Carlo, a closed-form approximation, which does not involve downgrade in credit ratings, to credit VaR for a large portfolio of similar loans is used in the minimum capital requirement formulas of the Basel Committee on Banking Supervision. Suppose the portfolio consists of M loans in the same risk class. Letting τi denote the time to default of the ith loan, we can express the proportion of loans that default during the period as ∑i=1M‍I{τi≤T}/M. Since the loans belong to the same risk class, their τi are identically distributed, with common distribution function G that is assumed to be continuous. Then G(τi) is uniform and Zi=Φ-1(G(τi)) is standard normal, where Φ is the standard normal distribution function. Vasicek [30, 31] and Schonbucher [32] assume the one-factor Gaussian model
(2.30)Zi=ρX+1-ρϵi,i=1,…,M,
in which X and ϵi are i.i.d. standard normal. Since τi≤T if and only if Zi≤Φ-1(G(T)), it follows that
(2.31)P(τi≤T∣X)=P{ϵi≤Φ-1(G(T))-ρX1-ρ∣X}=Φ(Φ-1(G(T))-ρX1-ρ).
Conditional on X, I{ϵi≤[Φ-1(G(T))-ρX]/1-ρ} are i.i.d. and therefore by the strong law of large numbers,
(2.32)M-1∑i=1M‍I{τi≤T}=M-1∑i=1M‍I{ϵi≤[Φ-1(G(T))-ρX]/1-ρ}→E[I{ϵi≤[Φ-1(G(T))-ρX]/1-ρ}∣X]=Φ(Φ-1(G(T))-ρX1-ρ)
as M→∞, with probability 1.LetL be the size of the loan portfolio and R the recovery rate. The proportional loss given default is 1-R and therefore the loss due to loan default in the period is L(1-R)∑i=1M‍I{τi≤T}/M. Using the approximation (2.32), the 100(1-α)% credit VaR over the period is
(2.33)L(1-R)Φ(Φ-1(G(T))+ρΦ-1(α)1-ρ),
since the (1-α)th quantile of X is Φ-1(1-α)=-Φ-1(α).


## 2.1. Survival Analysis, Intensity Processes, and Intensity Models for Default Risk

We begin with a brief review of survival analysis in biostatistics and reliability theory and then show its connections with widely used intensity models for pricing corporate bonds and other credit-sensitive securities. To pinpoint these connections, our review focuses on the martingale approach to survival analysis, as described in [5, 6]. A major theme of the statistical problems in survival analysis is to study the distribution of a failure time τ, based on possibly censored data from either a homogeneous population or a regression model with covariates. The survival function of τ is S(t)=P(τ≥t)=1-F(t-), where F(s)=P(τ≤s) and F(t-)=lims→t-F(s)=P(τ<t). Clearly, S(t) is a nonincreasing left-continuous function of t with S(0)=1 and limt→∞S(t)=0. If τ is absolutely continuous, the probability density function of τ is f(t)=-dS(t)/dt, and the hazard (or intensity) function is
(2.1)λ(t)=limh→0+P(t≤τ<t+h∣τ≥t)h=f(t)S(t).
Integrating (2.1) with respect to t yields
(2.2)S(t)=exp[-∫0t‍λ(u)du]=exp[-Λ(t)],
where Λ(t)=∫0t‍λ(s)ds is called the cumulative hazard function. Without assuming F to be continuous, we can define the cumulative hazard function by the Stieltjes integral
(2.3)Λ(t)=∫0t‍dF(s)1-F(s-).
Regarding the indicator functions I{τ≤t},t≥0, as a stochastic process, an important property of Λ is that
(2.4){I{τ≤t}-Λ(τ∧t),t≥0}isamartingale,
where τ∧t=min(τ,t). Fleming and Harrington [5] and Andersen et al. [6] give overviews of continuous-time martingales and their applications, via (2.4) and its analogs, to survival analysis. Based on n independent observations τi sampled from F, the nonparametric maximum likelihood estimate of F is the empirical distribution function F^(t)=(1/n)∑i=1n‍I{τi≤t}. However, in applications of survival analysis to biomedical and econometric data, the τi may not be fully observable because of censoring. Some subjects (or firms) may not fail during the observation period or may have been lost in followup during the period; the data on these individuals are said to be right-censored.Right-censoring can be formulated by introducing censoring variablesci that indicate the time of loss to followup. The censoring variable ci can also indicate the length of the observation period. The observations, therefore, are (Ti,δi),i=1,…,n, where Ti=min(τi,ci) and δi=I{τi≤ci} is the censoring indicator that indicates whether Ti is an actual failure time or is censored. Subject i is “at risk” at time s if Ti≥s (i.e., has not failed and has not been lost to followup prior to s). Let
(2.5)Y(s)=∑i=1n‍I{Ti≥s},N(s)=∑i=1n‍I{Ti≤s,δi=1}.
Note that Y(s) is the risk set size and ΔN(s)=N(s)-N(s-) is the number of observed deaths at time s, and that ΔN(s)/Y(s) is the analog of p^j. We use ΔN(s) to denote the jump size of N at s (i.e., ΔN(s)=N(s)-N(s-)) and also use the convention 0/0=0. The cumulative hazard function Λ(t) can be estimated by
(2.6)Λ^(t)=∑s≤t‍ΔN(s)Y(s)=∫0t‍I{Y(s)>0}Y(s)dN(s)
for the censored data (Ti,δi),1≤i≤n. This is called the Nelson-Aalen estimator [7, 8]. Analogous to (2.4),
(2.7){∫0tU(s)[dN(s)-Y(s)dΛ(s)],t≥0}isamartingale,
for every left-continuous stochastic process U(s); note that Y(s) is left-continuous. From this and the martingale central limit theorem,
(2.8)Λ^(t)-Λ(t){∫0t‍(I{Y(s)>0}/Y2(s))dN(s)}1/2→𝒟N(0,1)
as n→∞, where →𝒟 denotes convergence in distribution; see [5, 6].Altman [9] and others have developed mortality tables for loans and bonds by using methods in actuarial science. These mortality tables (also called life tables) are used by actuaries to set premiums for life insurance policies. Partition time into disjoint intervals I1=(0,t1],I2=(t1,t2], and so forth. A life table summarizes the mortality results of a large cohort of n subjects as follows:nj= number of subjects alive at the beginning of Ij,dj= number of deaths during Ij, lj= number lost to followup during Ij.It estimates pj=P(dyingduringIj∣aliveatthebeginningofIj) by p^j=dj/(nj-lj) so that the actuarial (life-table) estimate of P(τ>tk) is the product
(2.9)∏j=1k‍(1-p^j)=∏j=1k‍(1-djnj-lj).
Without discretizing the failure times, Kaplan and Meier [10] introduced the product-limit estimator of S(t)=P(τ≥t):
(2.10)S^(t)=∏s<t‍(1-ΔN(s)Y(s)).
Since N(s) has at most n jumps, the product in (2.10) is finite. Note that S^(t)=∏s<t‍(1-ΔΛ^(s)) by (2.6) and (2.10). As shown in [6],
(2.11)S(t)=∏s<t‍(1-dΛ(s))∀tsuchthatΛ(t)<∞.
The product in (2.11) is called the product integral of the nondecreasing function Λ(·) on [0,t]; it is defined by the limit of ∏{1-[Λ(ti)-Λ(ti-1)]}, where 0=t0<t1<⋯<tm=t is a partition of [0,t] and the product is in the natural order from left to right, as the mesh size max1≤i≤m|ti-ti-1| approaches 0. Moreover, by making use of the martingale central limit theorem, it is shown in [5, 6] that
(2.12)S^(t)-S(t)[S^(t){∫0t‍(I{Y(s)>0}/Y2(s))dN(s)}1/2]→𝒟N(0,1).


### 2.1.1. Regression Models for Hazard Functions with Covariates

In applications, one often wants to use a model forτ to predict future failures from a vector x(t) of predictors based on current and past observations; x(t) is called a time-varying covariate. When x(t)=x does not depend on t, it is called a time-independent (or baseline) covariate. In practice, some predictors may be time-independent, while other components of x(t) may be time-varying. Since prediction of future default from x(t) is relevant only if τ>t (i.e., if default has not occurred at or before t), one is interested in modeling the conditional distribution of τ given τ>t, for example, by relating the hazard function λ(t) to x(t). Noting that λ(t)dt is the probability that τ∈[t,t+dt], Cox [11] introduced the proportional hazards (or Cox regression) model
(2.13)λ(t)=λ0(t)exp(βTx(t)).
Putting x(t)=0 in (2.13) shows that λ0(·) is also a hazard function; it is called the baseline hazard. Cox regression fits (2.13) to the observed data (Ti,δi,xi(t)),t≤Ti,1≤i≤n. Lane et al. [12] apply Cox regression to predict bank default, using time-independent covariates.The parameters of (2.13) are β and the cumulative hazard function Λ0, which is an infinite-dimensional parameter. One way to estimate (β,Λ0) is to specify Λ0 via a parameter vector θ so that (β,θ) can be estimated by maximum likelihood. Since an observed failure contributes a factor fi(Ti) to the likelihood, where fi=λiSi is the density function of the τi while a censored failure contributes a factor Si(Ti) to the likelihood, and since logSi(t)=-Λi(t) and Λi(t)=Λ0(t)exp(βTxi(t)), the log-likelihood function can be written as
(2.14)l(β,θ)=∑i=1n{δilogfi(Ti)+(1-δi)logSi(Ti)}=∑i=1n{δi[logλ0(Ti;θ)+βTxi(Ti)]-Λ0(Ti;θ)eβTxi(Ti)};
see [13]. Commonly used parametric models for failure times are the exponential density θe-θt that has constant hazard θ>0 and the Weibull density θtθ-1exp(-tθ) that has hazard function θtθ-1 and cumulative hazard function tθ, with θ>0. Lee and Urrutia [14] use a Weibull distribution of default times to predict insolvency in the property-liability insurance industry. Hazard models with time-varying covariates have been used by McDonald and Van de Gucht [15] to model the timing of high-yield bond default and of exercising call options. Related survival analysis to predict bankruptcy has been carried out in [16–19].Instead of assuming a parametric model to estimate the baseline hazard functionλ0(·), Cox [11, 20] introduced a semiparametric method to estimate the finite-dimensional parameter β in the presence of an infinite-dimensional nuisance parameter λ0(·); it is semiparametric in the sense of being nonparametric in λ0 but parametric in β (in terms of βTxt). Cox's partial likelihood method decomposes the likelihood function into two factors, with one involving only β and the other involving both β and Λ. It estimates β by maximizing the following partial likelihood, which is a factor that only involves β. Order the observed censored failure times as τ(1)<⋯<τ(m), with m≤n. Let Cj denote the set of censored Ti's in the interval [τ(j-1),τ(j)), and let (j) denote the individual failing at τ(j), noting that with probability 1 there is only one failure at τ(j) because the failure time distributions have density functions. Let R(j)={i:Ti≥τ(j)} denote the risk set at τ(j). Cox's regression estimator β^ is the maximizer of the partial log-likelihood
(2.15)l(β)=∑j=1m{βTx(j)(τ(j))-log(∑i∈R(j)exp(βTx(i)(τ(j))))},
or equivalently, is the solution of (∂/∂β)l(β)=0. Martingale theory can be used to show that
(2.16)(-l¨(β))1/2(β^-β0)hasalimitingstandardnormaldistributionasn→∞,
where l¨(β) is the Hessian matrix of second partial derivatives (∂2/∂βk∂βh)l(β). One, therefore, can perform usual likelihood inference, treating the partial likelihood as a likelihood function. Moreover, even though β^ is based on partial likelihood, it has been shown to be asymptotically efficient; see [6]. Modifying Y(s) in (2.5) to
(2.17)Y(s)=∑i∈R(j)exp(βTxi(s))ats=τ(j)
leads to the Breslow estimator [21] of Λ0, which is again given by (2.6) but with Y(s) defined by (2.17).


### 2.1.2. Intensity-Based Models for Pricing Default Risk

Besides the important roles they have played in the statistical analysis of failure-time data as explained in the preceding section, intensity modeling has also become a standard approach to pricing corporate bonds and other credit-sensitive securities since the 1990s. Intensity models provide an elegant way to combine term-structure modeling for default-free bonds with the default risk and recovery features of corporate bonds. Chapter 10 of [22] gives an overview of interest rate markets and the term structure of interest rates, statistical methods to estimate the zero-coupon yield curve from the current prices of default-free Treasury bonds, arbitrage-free pricing theory of interest rate derivatives, short rate and other models under the risk-neutral (instead of the real-world) measure entailed by arbitrage-free pricing, and calibration of these models to market data involving actively traded interest rate derivatives. For corporate bonds, besides the interest rate processes, we also include an additional stochastic process for the default intensity of the bond. All the stochastic processes are defined on the same probability space in which the probability measure Q is a risk-neutral (or martingale) measure associated with arbitrage-free pricing. Expectations taken under this measure will be denoted by E^ or EQ, while E is used to denote expectation under the real-world (or physical) measure P.ACox process (also called stochastic Poisson process) is often used to model the default intensity of the bond's issuer. The default intensity λt is assumed to be governed by an exogenous stochastic process Xt,t≥0, so that λt=λ(Xt) and the stochastic dynamics of λt is specified through Xt. Let τ be the default time and Λ(t)=∫0t‍λ(Xs)ds. Then Λ(τ) is distributed as an exponential random variable ϵ1 with mean 1 that is independent of {Xs,s≥0} since P{Λ(τ)>t} = P{τ>Λ-1(t)} = S(Λ-1(t)) = exp{-Λ(Λ-1(t))} = e-t. Hence we can use λ(Xs) to generate τ by τ=inf{t:∫0t‍λ(Xs)ds≥ϵ1}. Consider a zero-coupon bond, with maturity date T and face value 1, issued by a firm at time 0. Assume that there is a short-rate process r(Xs) under the risk-neutral measure Q such that the default-free bond price is given by
(2.18)p(0,T)=E^{exp(-∫0Tr(Xs)ds)};
see [22, page 257]. Assuming λ(Xs) to be the intensity process for the default time τ of the firm and assuming zero recovery at default, the price of the defaultable bond at time 0 is
(2.19)π(0,T)=E^{I{τ>T}exp(-∫0Tr(Xs)ds)}=E^{E^[⋯∣Xs,s≤T]}=E^{exp(-∫0Tr(Xs)ds)E^[I{τ>T}∣Xs,s≤T]}=E^{exp(-∫0T(r+λ)(Xs)ds)},
since E^[I{τ>T}∣Xs,s≤T]=Q{∫0T‍λ(Xs)<ϵ1∣Xs,s≤T}=exp(-∫0T‍λ(Xs)ds), noting that P(ϵ1>x)=e-x and that ϵ1 is independent of {Xs,s≥0}. Thus (2.19) replaces the short rate in the default-free bond pricing formula by the default-adjusted short rate (r+λ)(Xs). Moreover, closed-form expressions are available for π(0,T) if one uses affine models for both r and λ; see [23, 24]. The short rate process and the intensity process have unknown parameters which can be encoded in a vector θ and can be estimated from the bond prices by using least squares in the same way as that for Treasury bonds; see [22]. Note that unlike survival analysis in biostatistics that uses censored failure-time data to estimate the intensity in (2.13) or its integral in (2.3), estimation of the parameters of the intensity process is based not on observed defaults but rather on the bond prices that involve the risk-neutral (rather than the real-world) measure.In practice, recovery rates also need to be considered and the vectorXs may include predictors that make r(Xs) and λ(Xs) much more complex than the affine model that gives an explicit formula for the integral in (2.19). Monte Carlo methods are used instead to evaluate (2.19) in these more complex situations. Since corporate bonds typically promise cash flows, which mimic those of Treasury bonds, it is natural to consider yields when one compares the effects of different recovery assumptions at default. The yield at date t of a zero-coupon bond with price Bt, face value A, and maturity T is defined by
(2.20)y(t,T)=1T-tlogABt,
and the difference between the yield of a defaultable bond and a corresponding Treasury bond is referred to as the credit spread or yield spread. The first type of recovery assumptions is recovery of market value. This measures the change in market prices before and after the default, representing the loss in the bond's value associated with the default. More generally, consider the price process of a defaultable contingent claim V promising a payoff of f(XT) at T. The claim is said to have a fractional recovery of market value of δ at default time τ if the amount recovered in the event of default is equal to
(2.21)h(τ)=δV(τ-)forτ≤T,
where V(τ-) is the value of the claim just prior to default and δ∈[0,1). With this recovery assumption, the price at t of the defaultable contingent claim (if there has not been a default up to date t) is
(2.22)V(t)=E^[f(XT)exp(-∫tT(r+(1-δ)λ)(Xs)ds)∣Xs,s≤T],
which is an extension of (2.19). A variant of recovery of market value is recovery of face value, which is the closest to legal practice in the sense that debt with the same priority is assigned a fractional recovery depending on the outstanding notional amount but not on maturity or coupon. It is also the measure typically used in rating-agency studies. The third type of recovery assumptions is recovery of treasury, under which the corporate bond in default is replaced with a treasury bond with the same maturity but a reduced payment. Unlike recovery of face value, it tries to correct for the fact that amounts of principal with long maturity should be discounted more than principal payments with short maturity.


## 2.1.1. Regression Models for Hazard Functions with Covariates

In applications, one often wants to use a model forτ to predict future failures from a vector x(t) of predictors based on current and past observations; x(t) is called a time-varying covariate. When x(t)=x does not depend on t, it is called a time-independent (or baseline) covariate. In practice, some predictors may be time-independent, while other components of x(t) may be time-varying. Since prediction of future default from x(t) is relevant only if τ>t (i.e., if default has not occurred at or before t), one is interested in modeling the conditional distribution of τ given τ>t, for example, by relating the hazard function λ(t) to x(t). Noting that λ(t)dt is the probability that τ∈[t,t+dt], Cox [11] introduced the proportional hazards (or Cox regression) model
(2.13)λ(t)=λ0(t)exp(βTx(t)).
Putting x(t)=0 in (2.13) shows that λ0(·) is also a hazard function; it is called the baseline hazard. Cox regression fits (2.13) to the observed data (Ti,δi,xi(t)),t≤Ti,1≤i≤n. Lane et al. [12] apply Cox regression to predict bank default, using time-independent covariates.The parameters of (2.13) are β and the cumulative hazard function Λ0, which is an infinite-dimensional parameter. One way to estimate (β,Λ0) is to specify Λ0 via a parameter vector θ so that (β,θ) can be estimated by maximum likelihood. Since an observed failure contributes a factor fi(Ti) to the likelihood, where fi=λiSi is the density function of the τi while a censored failure contributes a factor Si(Ti) to the likelihood, and since logSi(t)=-Λi(t) and Λi(t)=Λ0(t)exp(βTxi(t)), the log-likelihood function can be written as
(2.14)l(β,θ)=∑i=1n{δilogfi(Ti)+(1-δi)logSi(Ti)}=∑i=1n{δi[logλ0(Ti;θ)+βTxi(Ti)]-Λ0(Ti;θ)eβTxi(Ti)};
see [13]. Commonly used parametric models for failure times are the exponential density θe-θt that has constant hazard θ>0 and the Weibull density θtθ-1exp(-tθ) that has hazard function θtθ-1 and cumulative hazard function tθ, with θ>0. Lee and Urrutia [14] use a Weibull distribution of default times to predict insolvency in the property-liability insurance industry. Hazard models with time-varying covariates have been used by McDonald and Van de Gucht [15] to model the timing of high-yield bond default and of exercising call options. Related survival analysis to predict bankruptcy has been carried out in [16–19].Instead of assuming a parametric model to estimate the baseline hazard functionλ0(·), Cox [11, 20] introduced a semiparametric method to estimate the finite-dimensional parameter β in the presence of an infinite-dimensional nuisance parameter λ0(·); it is semiparametric in the sense of being nonparametric in λ0 but parametric in β (in terms of βTxt). Cox's partial likelihood method decomposes the likelihood function into two factors, with one involving only β and the other involving both β and Λ. It estimates β by maximizing the following partial likelihood, which is a factor that only involves β. Order the observed censored failure times as τ(1)<⋯<τ(m), with m≤n. Let Cj denote the set of censored Ti's in the interval [τ(j-1),τ(j)), and let (j) denote the individual failing at τ(j), noting that with probability 1 there is only one failure at τ(j) because the failure time distributions have density functions. Let R(j)={i:Ti≥τ(j)} denote the risk set at τ(j). Cox's regression estimator β^ is the maximizer of the partial log-likelihood
(2.15)l(β)=∑j=1m{βTx(j)(τ(j))-log(∑i∈R(j)exp(βTx(i)(τ(j))))},
or equivalently, is the solution of (∂/∂β)l(β)=0. Martingale theory can be used to show that
(2.16)(-l¨(β))1/2(β^-β0)hasalimitingstandardnormaldistributionasn→∞,
where l¨(β) is the Hessian matrix of second partial derivatives (∂2/∂βk∂βh)l(β). One, therefore, can perform usual likelihood inference, treating the partial likelihood as a likelihood function. Moreover, even though β^ is based on partial likelihood, it has been shown to be asymptotically efficient; see [6]. Modifying Y(s) in (2.5) to
(2.17)Y(s)=∑i∈R(j)exp(βTxi(s))ats=τ(j)
leads to the Breslow estimator [21] of Λ0, which is again given by (2.6) but with Y(s) defined by (2.17).


## 2.1.2. Intensity-Based Models for Pricing Default Risk

Besides the important roles they have played in the statistical analysis of failure-time data as explained in the preceding section, intensity modeling has also become a standard approach to pricing corporate bonds and other credit-sensitive securities since the 1990s. Intensity models provide an elegant way to combine term-structure modeling for default-free bonds with the default risk and recovery features of corporate bonds. Chapter 10 of [22] gives an overview of interest rate markets and the term structure of interest rates, statistical methods to estimate the zero-coupon yield curve from the current prices of default-free Treasury bonds, arbitrage-free pricing theory of interest rate derivatives, short rate and other models under the risk-neutral (instead of the real-world) measure entailed by arbitrage-free pricing, and calibration of these models to market data involving actively traded interest rate derivatives. For corporate bonds, besides the interest rate processes, we also include an additional stochastic process for the default intensity of the bond. All the stochastic processes are defined on the same probability space in which the probability measure Q is a risk-neutral (or martingale) measure associated with arbitrage-free pricing. Expectations taken under this measure will be denoted by E^ or EQ, while E is used to denote expectation under the real-world (or physical) measure P.ACox process (also called stochastic Poisson process) is often used to model the default intensity of the bond's issuer. The default intensity λt is assumed to be governed by an exogenous stochastic process Xt,t≥0, so that λt=λ(Xt) and the stochastic dynamics of λt is specified through Xt. Let τ be the default time and Λ(t)=∫0t‍λ(Xs)ds. Then Λ(τ) is distributed as an exponential random variable ϵ1 with mean 1 that is independent of {Xs,s≥0} since P{Λ(τ)>t} = P{τ>Λ-1(t)} = S(Λ-1(t)) = exp{-Λ(Λ-1(t))} = e-t. Hence we can use λ(Xs) to generate τ by τ=inf{t:∫0t‍λ(Xs)ds≥ϵ1}. Consider a zero-coupon bond, with maturity date T and face value 1, issued by a firm at time 0. Assume that there is a short-rate process r(Xs) under the risk-neutral measure Q such that the default-free bond price is given by
(2.18)p(0,T)=E^{exp(-∫0Tr(Xs)ds)};
see [22, page 257]. Assuming λ(Xs) to be the intensity process for the default time τ of the firm and assuming zero recovery at default, the price of the defaultable bond at time 0 is
(2.19)π(0,T)=E^{I{τ>T}exp(-∫0Tr(Xs)ds)}=E^{E^[⋯∣Xs,s≤T]}=E^{exp(-∫0Tr(Xs)ds)E^[I{τ>T}∣Xs,s≤T]}=E^{exp(-∫0T(r+λ)(Xs)ds)},
since E^[I{τ>T}∣Xs,s≤T]=Q{∫0T‍λ(Xs)<ϵ1∣Xs,s≤T}=exp(-∫0T‍λ(Xs)ds), noting that P(ϵ1>x)=e-x and that ϵ1 is independent of {Xs,s≥0}. Thus (2.19) replaces the short rate in the default-free bond pricing formula by the default-adjusted short rate (r+λ)(Xs). Moreover, closed-form expressions are available for π(0,T) if one uses affine models for both r and λ; see [23, 24]. The short rate process and the intensity process have unknown parameters which can be encoded in a vector θ and can be estimated from the bond prices by using least squares in the same way as that for Treasury bonds; see [22]. Note that unlike survival analysis in biostatistics that uses censored failure-time data to estimate the intensity in (2.13) or its integral in (2.3), estimation of the parameters of the intensity process is based not on observed defaults but rather on the bond prices that involve the risk-neutral (rather than the real-world) measure.In practice, recovery rates also need to be considered and the vectorXs may include predictors that make r(Xs) and λ(Xs) much more complex than the affine model that gives an explicit formula for the integral in (2.19). Monte Carlo methods are used instead to evaluate (2.19) in these more complex situations. Since corporate bonds typically promise cash flows, which mimic those of Treasury bonds, it is natural to consider yields when one compares the effects of different recovery assumptions at default. The yield at date t of a zero-coupon bond with price Bt, face value A, and maturity T is defined by
(2.20)y(t,T)=1T-tlogABt,
and the difference between the yield of a defaultable bond and a corresponding Treasury bond is referred to as the credit spread or yield spread. The first type of recovery assumptions is recovery of market value. This measures the change in market prices before and after the default, representing the loss in the bond's value associated with the default. More generally, consider the price process of a defaultable contingent claim V promising a payoff of f(XT) at T. The claim is said to have a fractional recovery of market value of δ at default time τ if the amount recovered in the event of default is equal to
(2.21)h(τ)=δV(τ-)forτ≤T,
where V(τ-) is the value of the claim just prior to default and δ∈[0,1). With this recovery assumption, the price at t of the defaultable contingent claim (if there has not been a default up to date t) is
(2.22)V(t)=E^[f(XT)exp(-∫tT(r+(1-δ)λ)(Xs)ds)∣Xs,s≤T],
which is an extension of (2.19). A variant of recovery of market value is recovery of face value, which is the closest to legal practice in the sense that debt with the same priority is assigned a fractional recovery depending on the outstanding notional amount but not on maturity or coupon. It is also the measure typically used in rating-agency studies. The third type of recovery assumptions is recovery of treasury, under which the corporate bond in default is replaced with a treasury bond with the same maturity but a reduced payment. Unlike recovery of face value, it tries to correct for the fact that amounts of principal with long maturity should be discounted more than principal payments with short maturity.


## 2.2. Econometric Modeling of Default Risk

The intensity-based models in the preceding section are sometimes calledreduced-form models, using a terminology in econometrics that refers to models of the relationship of an economic variable to exogenous variables. In contrast, structural models attempt to derive the relationship from some underlying economic theory. Structural models for the default of a firm assume that default is triggered when the firm's asset and debt values enter some region whose boundary is called a “default barrier.”


### 2.2.1. Merton's Model of Debt and Equity Values and the KMV Credit Monitor Model

As in the Black-Scholes model described in Section 8.1 of [22], Merton [25] assumes a market with continuous trading, no transaction costs, unlimited short selling and perfectly divisible assets, and with a risk-free asset that has constant interest rate r. The asset value At of the firm is assumed to follow a geometric Brownian motion (GBM) dAt/At=μdt+σdwt. Suppose at time 0 the firm has issued two kinds of contingent claims: equity and debt, in which debt is assumed to be a zero-coupon bond with face value D and maturity date T. We can regard the firm as being run by equity owners. At maturity date T of the bond, equity owners pay D to retain ownership if At≥D, and would not pay D if At<D since they have limited liability. In the latter case, bond holders receive a recovery of AT instead of the promised payment D. Therefore, the equity value at time t is given by the Black-Scholes formula
(2.23)Ve(t)=EQ{e-r(T-t)(AT-D)+∣At}=AtΦ(d1)-De-r(T-t)Φ(d2),
where Q is the risk-neutral measure under which the drift of the GBM is r instead of μ, and
(2.24)d1=log(At/D)+(r+σ2/2)(T-t)σT-t,d2=d1-σT-t;
see [22, page 186]. Similarly, the debt value at time t is
(2.25)Vd(t)=EQ{e-r(T-t)min(AT,D)∣At}=At-Ve(t),
since min(At,D)=At-(At-D)+.In the 1990s, the KMV Corporation of San Francisco (now merged into Moody's KMV) introduced its Credit Monitor Model for default prediction using Merton's framework. To begin with, the GBM for the asset value can be expressed in terms of the Brownian motionwt by At=A0exp{(μ-(1/2)σ2)t+σwt}; see [22, page 66]. Therefore, the probability of default on the debt at maturity is
(2.26)P(AT<D)=P{(μ-σ22)T+σwT<log(DA0)}=Φ(-logA0+(μ-σ2/2)T-logDσT)=Φ(-E(logAt)-logDσT).
This leads to KMV's distance to default defined by
(2.27)DD={E(logAT)-logD}σT
so that the predicted probability of default under Merton's model is Φ(-DD). Since σT is the standard deviation of logAT, DD can be regarded as the number of standard deviations of the log asset value from the log liability. Instead of using the normal distribution for logAT to predict the default probability as in (2.26), the KMV Credit Monitor Model uses a large proprietary database of firms and firm defaults to construct an EDF (expected default frequency) score for a firm with κ standard deviations away from the default (i.e., DD =κ) at the beginning of a future period. The idea is to stratify the firms based on their DD's and estimate empirically the default probability of the firm from its associated stratum. The firm's DD is evaluated from its balance sheet and option models for its debt.


### 2.2.2. Black-Cox Framework of Default Barriers and Other Extensions of Merton's Model

Black and Cox [26] extend the Merton model by allowing defaults to occur prior to the maturity of the bond; specifically, when the asset value At hits a nonrandom, continuous boundary g(t)≤D. They also allow the equity owners to be paid dividends, at rate q, so that under the risk-neutral measure Q the GBM for At has drift r-q; see [22, page 186]. Let τ=inf{0≤t≤T:At=g(t)}. The payoff to bond holders is
(2.28)min(AT,D)I{τ>T}+g(τ)I{τ≤T},
in which the first summand reflects that the bond holders receive min(AT,D) at the bond's maturity if default has not occurred up to that time, while the second summand reflects that the bond holders receive the recovery of Aτ=g(τ) upon default. Similarly, the payoff to equity holders is (AT-D)+I{τ>T}.Stochastic models for the default-free interest rate dynamics (see [22, Chapter 10]) have been used to replace the constant r in Merton's model. Another extension is to allow coupon payments to the bond holders. More general jump diffusion models have been introduced to replace GBM in Merton's model. Lando [24] provides an overview of these extensions, which also have been considered in conjunction with default barriers.


## 2.2.1. Merton's Model of Debt and Equity Values and the KMV Credit Monitor Model

As in the Black-Scholes model described in Section 8.1 of [22], Merton [25] assumes a market with continuous trading, no transaction costs, unlimited short selling and perfectly divisible assets, and with a risk-free asset that has constant interest rate r. The asset value At of the firm is assumed to follow a geometric Brownian motion (GBM) dAt/At=μdt+σdwt. Suppose at time 0 the firm has issued two kinds of contingent claims: equity and debt, in which debt is assumed to be a zero-coupon bond with face value D and maturity date T. We can regard the firm as being run by equity owners. At maturity date T of the bond, equity owners pay D to retain ownership if At≥D, and would not pay D if At<D since they have limited liability. In the latter case, bond holders receive a recovery of AT instead of the promised payment D. Therefore, the equity value at time t is given by the Black-Scholes formula
(2.23)Ve(t)=EQ{e-r(T-t)(AT-D)+∣At}=AtΦ(d1)-De-r(T-t)Φ(d2),
where Q is the risk-neutral measure under which the drift of the GBM is r instead of μ, and
(2.24)d1=log(At/D)+(r+σ2/2)(T-t)σT-t,d2=d1-σT-t;
see [22, page 186]. Similarly, the debt value at time t is
(2.25)Vd(t)=EQ{e-r(T-t)min(AT,D)∣At}=At-Ve(t),
since min(At,D)=At-(At-D)+.In the 1990s, the KMV Corporation of San Francisco (now merged into Moody's KMV) introduced its Credit Monitor Model for default prediction using Merton's framework. To begin with, the GBM for the asset value can be expressed in terms of the Brownian motionwt by At=A0exp{(μ-(1/2)σ2)t+σwt}; see [22, page 66]. Therefore, the probability of default on the debt at maturity is
(2.26)P(AT<D)=P{(μ-σ22)T+σwT<log(DA0)}=Φ(-logA0+(μ-σ2/2)T-logDσT)=Φ(-E(logAt)-logDσT).
This leads to KMV's distance to default defined by
(2.27)DD={E(logAT)-logD}σT
so that the predicted probability of default under Merton's model is Φ(-DD). Since σT is the standard deviation of logAT, DD can be regarded as the number of standard deviations of the log asset value from the log liability. Instead of using the normal distribution for logAT to predict the default probability as in (2.26), the KMV Credit Monitor Model uses a large proprietary database of firms and firm defaults to construct an EDF (expected default frequency) score for a firm with κ standard deviations away from the default (i.e., DD =κ) at the beginning of a future period. The idea is to stratify the firms based on their DD's and estimate empirically the default probability of the firm from its associated stratum. The firm's DD is evaluated from its balance sheet and option models for its debt.


## 2.2.2. Black-Cox Framework of Default Barriers and Other Extensions of Merton's Model

Black and Cox [26] extend the Merton model by allowing defaults to occur prior to the maturity of the bond; specifically, when the asset value At hits a nonrandom, continuous boundary g(t)≤D. They also allow the equity owners to be paid dividends, at rate q, so that under the risk-neutral measure Q the GBM for At has drift r-q; see [22, page 186]. Let τ=inf{0≤t≤T:At=g(t)}. The payoff to bond holders is
(2.28)min(AT,D)I{τ>T}+g(τ)I{τ≤T},
in which the first summand reflects that the bond holders receive min(AT,D) at the bond's maturity if default has not occurred up to that time, while the second summand reflects that the bond holders receive the recovery of Aτ=g(τ) upon default. Similarly, the payoff to equity holders is (AT-D)+I{τ>T}.Stochastic models for the default-free interest rate dynamics (see [22, Chapter 10]) have been used to replace the constant r in Merton's model. Another extension is to allow coupon payments to the bond holders. More general jump diffusion models have been introduced to replace GBM in Merton's model. Lando [24] provides an overview of these extensions, which also have been considered in conjunction with default barriers.


## 2.3. Credit Ratings and Credit Value at Risk

Credit scoring has been used in the banking industry since the 1940s to assess the creditworthiness of loan applicants. The basic idea is to assign scores to loan applications based on their attributes such as gender, marital status, education level, income, profession, and credit history. Such scoring systems, calledscorecards, are usually developed by in-house credit departments and/or consulting firms. One of the most widely used credit scoring systems was developed by Fair Isaac & Company (FICO), which built credit scoring systems for US banks and retail stores in the 1960s and 1970s and then expanded those systems to meet the needs of other industries and to evaluate the credit of borrowers. Credit scoring systems have now been developed for virtually all types of credit analysis, from consumer credit to commercial loans. The idea is to preidentify certain key factors that determine the probability of default and combine or weight them into a quantitative score.Credit ratings of the creditworthiness of corporate bonds are provided by rating agencies such as Moody's, S&P, and Fitch. Using Moody's system, the ratings are Aaa, Aa, A, Baa, Ba, B, Caa, Ca, and C and “investment grade” ratings are Baa or above. The corresponding S&P and Fitch ratings are AAA, AA, A, BBB, BB, B, CCC, CC, and C. For finer ratings, Moody's divides the Aa rating category into Aa1, Aa2, Aa3, the A category into A1, A2, A3, and so forth, and S&P and Fitch divide AA into AA+, AA, and AA−, A into A+, A, and A−, and so forth. Although methodologies and standards differ from one rating agency to another, regulators generally do not make distinctions among the agencies. Furthermore, although there is a strong similarity between the rating systems of Moody's and S&P, different agencies might assign slightly different ratings for the same bond. The ratings given by rating agencies are used by several credit risk software models such as CreditMetrics of J. P. Morgan, a system that evaluates risks individually or across an entire portfolio. Most corporations approach rating agencies to request a rating prior to sale or registration of a debt issue. S&P assigns and publishes ratings for all public corporate debt issues over $50 million, with or without a request from the issuer; but in all instances, S&P analytical staff first contact the issuer to call for meetings.


### 2.3.1. CreditMetrics, Credit Ratings Migration, and Product-Limit Estimators

CreditMetrics was introduced in 1997 by J. P. Morgan and its cosponsors (including KMV, UBS, Bank of America, and others). The system measures the risk and potential losses of loans and portfolios of bonds and loans over a period (e.g., a year) by simulation and analytic methods. In particular, it uses Monte Carlo simulations to estimate credit rating changes of all counterparties. Historical data provided by rating agencies on credit rating changes can be used to estimate the transition probability matrix for credit rating changes. In the Basel Accord for bank regulation, banks are required to construct these transition probability matrices based on their own data to stress test their portfolios. Some risk management systems such as CreditMetrics use these credit rating transition matrices for Monte Carlo simulations of credit losses due to credit ratings migration in a future period. A Markov chain model of ratings migration assumes a finite state space𝒮={1,…,K} consisting of different rating classes, in which state 1 denotes the best credit rating class and state K represents the default case. The K×K probability transition matrix of the Markov chain is
(2.29)P=(p11p12⋯p1Kp21p22⋯p2K⋮⋮⋮⋮pK-1,1pK-1,2⋯pK-1,K00⋯1),
in which pii=1-∑j≠i‍pij for all i and pij represents the actual probability of moving to state j from initial rating state i in one time step. Note that the default state K is an absorbing state. The estimates of credit transition matrices published by rating agencies use this Markov chain framework.Suppose there areNi firms in a given rating category i at the beginning of the period and that Nij firms migrate to the category j at the end of the period. An obvious estimate of pij is p^ij=Nij/Ni for j≠i. However, there are serious difficulties with this estimate in practice. First, if no transition from rating class i to j is observed in the period, then the estimated transition probability p^ij becomes 0. Second, the assumption of homogeneity over time is too strong for credit transition matrices. The product-limit estimator (2.10) for the survival function has been generalized by Aalen [27], Aalen and Johansen [28], and Fleming [29] to estimate the transition probability matrix P(s,t) of a nonhomogeneous continuous-time Markov process on a finite state-space. Analogous to (2.11), the transition probability matrix from time s to t can be expressed as the product integral P(s,t)=∏(s,t]‍(I+dA(u)), and the components Aij of the cumulative transition intensity matrix can be estimated by the Nelson-Aalen estimator of the form (2.6); see Section  IV.4 of [6].To simulate the distribution of credit losses due to credit rating changes of two counterparties, say one with initial Aaa rating (state 1) and the other with initial B rating, over a 1-year period, we can use the estimateP^ of the transition matrix (2.29) to sample the new ratings for both counterparties in the following way. Without assuming the rating changes of the two counterparties to be independent, a Gaussian copula model, which will be described in Section 3.3, is typically used to simulate the joint probability distribution of the rating changes. Therefore, we first generate two standard normal random variables ZA and ZB with correlation ρ. Noting that Φ(ZA) is uniform, the new rating of the Aaa company remains Aaa if Φ(ZA)<p^11, moves to states 2 (i.e., Aa) if p^11≤Φ(ZA)<p^11+p^12, and so forth. This is tantamount to generating the state transitions of the Aaa company by the first row of the probability transition matrix P^. While the state transitions of the B company are likewise generated by the corresponding row of P^, the state transitions of the two companies are not independent and their dependence structure is inherited from the correlation of ZA and ZB. From the simulated credit rating changes of the two counterparties we can revalue their outstanding contracts to determine the total credit losses during the year. From a large number of such simulations, a probability distribution for credit losses is obtained, which can be used to evaluate credit VaR (Values at Risk) for risk management.


### 2.3.2. Credit VaR and a One-Factor Default Model

A100(1-α)% credit VaR over a period of T years is the (1-α)th quantile of the credit loss within the period. J. P. Morgan's CreditMetrics defines credit losses as those arising not only from defaults but also from credit ratings downgrade. This involves estimating the probability distribution of credit losses by Monte Carlo simulation of defaults and credit rating changes of all counterparties. Instead of Monte Carlo, a closed-form approximation, which does not involve downgrade in credit ratings, to credit VaR for a large portfolio of similar loans is used in the minimum capital requirement formulas of the Basel Committee on Banking Supervision. Suppose the portfolio consists of M loans in the same risk class. Letting τi denote the time to default of the ith loan, we can express the proportion of loans that default during the period as ∑i=1M‍I{τi≤T}/M. Since the loans belong to the same risk class, their τi are identically distributed, with common distribution function G that is assumed to be continuous. Then G(τi) is uniform and Zi=Φ-1(G(τi)) is standard normal, where Φ is the standard normal distribution function. Vasicek [30, 31] and Schonbucher [32] assume the one-factor Gaussian model
(2.30)Zi=ρX+1-ρϵi,i=1,…,M,
in which X and ϵi are i.i.d. standard normal. Since τi≤T if and only if Zi≤Φ-1(G(T)), it follows that
(2.31)P(τi≤T∣X)=P{ϵi≤Φ-1(G(T))-ρX1-ρ∣X}=Φ(Φ-1(G(T))-ρX1-ρ).
Conditional on X, I{ϵi≤[Φ-1(G(T))-ρX]/1-ρ} are i.i.d. and therefore by the strong law of large numbers,
(2.32)M-1∑i=1M‍I{τi≤T}=M-1∑i=1M‍I{ϵi≤[Φ-1(G(T))-ρX]/1-ρ}→E[I{ϵi≤[Φ-1(G(T))-ρX]/1-ρ}∣X]=Φ(Φ-1(G(T))-ρX1-ρ)
as M→∞, with probability 1.LetL be the size of the loan portfolio and R the recovery rate. The proportional loss given default is 1-R and therefore the loss due to loan default in the period is L(1-R)∑i=1M‍I{τi≤T}/M. Using the approximation (2.32), the 100(1-α)% credit VaR over the period is
(2.33)L(1-R)Φ(Φ-1(G(T))+ρΦ-1(α)1-ρ),
since the (1-α)th quantile of X is Φ-1(1-α)=-Φ-1(α).


## 2.3.1. CreditMetrics, Credit Ratings Migration, and Product-Limit Estimators

CreditMetrics was introduced in 1997 by J. P. Morgan and its cosponsors (including KMV, UBS, Bank of America, and others). The system measures the risk and potential losses of loans and portfolios of bonds and loans over a period (e.g., a year) by simulation and analytic methods. In particular, it uses Monte Carlo simulations to estimate credit rating changes of all counterparties. Historical data provided by rating agencies on credit rating changes can be used to estimate the transition probability matrix for credit rating changes. In the Basel Accord for bank regulation, banks are required to construct these transition probability matrices based on their own data to stress test their portfolios. Some risk management systems such as CreditMetrics use these credit rating transition matrices for Monte Carlo simulations of credit losses due to credit ratings migration in a future period. A Markov chain model of ratings migration assumes a finite state space𝒮={1,…,K} consisting of different rating classes, in which state 1 denotes the best credit rating class and state K represents the default case. The K×K probability transition matrix of the Markov chain is
(2.29)P=(p11p12⋯p1Kp21p22⋯p2K⋮⋮⋮⋮pK-1,1pK-1,2⋯pK-1,K00⋯1),
in which pii=1-∑j≠i‍pij for all i and pij represents the actual probability of moving to state j from initial rating state i in one time step. Note that the default state K is an absorbing state. The estimates of credit transition matrices published by rating agencies use this Markov chain framework.Suppose there areNi firms in a given rating category i at the beginning of the period and that Nij firms migrate to the category j at the end of the period. An obvious estimate of pij is p^ij=Nij/Ni for j≠i. However, there are serious difficulties with this estimate in practice. First, if no transition from rating class i to j is observed in the period, then the estimated transition probability p^ij becomes 0. Second, the assumption of homogeneity over time is too strong for credit transition matrices. The product-limit estimator (2.10) for the survival function has been generalized by Aalen [27], Aalen and Johansen [28], and Fleming [29] to estimate the transition probability matrix P(s,t) of a nonhomogeneous continuous-time Markov process on a finite state-space. Analogous to (2.11), the transition probability matrix from time s to t can be expressed as the product integral P(s,t)=∏(s,t]‍(I+dA(u)), and the components Aij of the cumulative transition intensity matrix can be estimated by the Nelson-Aalen estimator of the form (2.6); see Section  IV.4 of [6].To simulate the distribution of credit losses due to credit rating changes of two counterparties, say one with initial Aaa rating (state 1) and the other with initial B rating, over a 1-year period, we can use the estimateP^ of the transition matrix (2.29) to sample the new ratings for both counterparties in the following way. Without assuming the rating changes of the two counterparties to be independent, a Gaussian copula model, which will be described in Section 3.3, is typically used to simulate the joint probability distribution of the rating changes. Therefore, we first generate two standard normal random variables ZA and ZB with correlation ρ. Noting that Φ(ZA) is uniform, the new rating of the Aaa company remains Aaa if Φ(ZA)<p^11, moves to states 2 (i.e., Aa) if p^11≤Φ(ZA)<p^11+p^12, and so forth. This is tantamount to generating the state transitions of the Aaa company by the first row of the probability transition matrix P^. While the state transitions of the B company are likewise generated by the corresponding row of P^, the state transitions of the two companies are not independent and their dependence structure is inherited from the correlation of ZA and ZB. From the simulated credit rating changes of the two counterparties we can revalue their outstanding contracts to determine the total credit losses during the year. From a large number of such simulations, a probability distribution for credit losses is obtained, which can be used to evaluate credit VaR (Values at Risk) for risk management.


## 2.3.2. Credit VaR and a One-Factor Default Model

A100(1-α)% credit VaR over a period of T years is the (1-α)th quantile of the credit loss within the period. J. P. Morgan's CreditMetrics defines credit losses as those arising not only from defaults but also from credit ratings downgrade. This involves estimating the probability distribution of credit losses by Monte Carlo simulation of defaults and credit rating changes of all counterparties. Instead of Monte Carlo, a closed-form approximation, which does not involve downgrade in credit ratings, to credit VaR for a large portfolio of similar loans is used in the minimum capital requirement formulas of the Basel Committee on Banking Supervision. Suppose the portfolio consists of M loans in the same risk class. Letting τi denote the time to default of the ith loan, we can express the proportion of loans that default during the period as ∑i=1M‍I{τi≤T}/M. Since the loans belong to the same risk class, their τi are identically distributed, with common distribution function G that is assumed to be continuous. Then G(τi) is uniform and Zi=Φ-1(G(τi)) is standard normal, where Φ is the standard normal distribution function. Vasicek [30, 31] and Schonbucher [32] assume the one-factor Gaussian model
(2.30)Zi=ρX+1-ρϵi,i=1,…,M,
in which X and ϵi are i.i.d. standard normal. Since τi≤T if and only if Zi≤Φ-1(G(T)), it follows that
(2.31)P(τi≤T∣X)=P{ϵi≤Φ-1(G(T))-ρX1-ρ∣X}=Φ(Φ-1(G(T))-ρX1-ρ).
Conditional on X, I{ϵi≤[Φ-1(G(T))-ρX]/1-ρ} are i.i.d. and therefore by the strong law of large numbers,
(2.32)M-1∑i=1M‍I{τi≤T}=M-1∑i=1M‍I{ϵi≤[Φ-1(G(T))-ρX]/1-ρ}→E[I{ϵi≤[Φ-1(G(T))-ρX]/1-ρ}∣X]=Φ(Φ-1(G(T))-ρX1-ρ)
as M→∞, with probability 1.LetL be the size of the loan portfolio and R the recovery rate. The proportional loss given default is 1-R and therefore the loss due to loan default in the period is L(1-R)∑i=1M‍I{τi≤T}/M. Using the approximation (2.32), the 100(1-α)% credit VaR over the period is
(2.33)L(1-R)Φ(Φ-1(G(T))+ρΦ-1(α)1-ρ),
since the (1-α)th quantile of X is Φ-1(1-α)=-Φ-1(α).


## 3. Credit Derivatives: An Overview

### 3.1. Pricing and Hedging with Credit Default Swaps

Credit default swaps (CDSs) have been described in the second paragraph of Section1. The total amount paid per year by the buyer of the CDS, as a percentage of the notional principal per annum, is called the CDS spread. Suppose the reference bond has coupon dates t1,…,tM=T and the protection buyer pays a constant premium c at these dates prior to τ, the time of the credit event. Assume that the face value of the bond is 1 and that the protection seller pays 1-δ in the case of a credit event. Then the value vb of the premium leg of the swap is
(3.1)vb=cE^{∑i=1Mexp(-∫0tirsds)I{τ>ti}}=cE^{∑i=1Mexp(-∫0ti(rs+λs)ds)}=c∑i=1Mπ(0,ti),
where rs is the short rate process, λs is the default intensity process, and π(0,ti) is the price of the defaultable bond with zero recovery and maturity ti, using the same argument as that in (2.19). Similarly the value vs of the default leg of the swap, to be paid by the protection seller, is
(3.2)vs=(1-δ)E^{exp(-∫0τrsds)I{τ≤T}}=(1-δ)∫0TE^{exp(-∫0τrsds)}f^(t)dt,
where f^ is the density function of τ under the risk-neutral measure Q, assuming that the short rate process and the default intensity process are independent under Q. Note that the survival function of τ under Q is
(3.3)S^(t)=Q(τ≥t)=E^exp(-∫0tλsds).
Letting λ^=f^/S^ be the hazard function of τ under Q, it then follows from (3.2) that
(3.4)vs=(1-δ)∫0Tp(0,t)λ^(t)S^(t)dt,
where p(0,t) is the price of the default-free bond; see (2.18). Note that under this independence assumption on the short rate and default intensity processes, we can further simplify π(0,ti) in (3.1) to
(3.5)π(0,ti)=E^[exp(-∫0tirsds)]E^[exp(-∫0tiλsds)]=p(0,ti)S^(ti).
In the absence of arbitrage, the premium cds of the default swap is the premium c for which vb=vs. Therefore, combining (3.1) with (3.4) yields
(3.6)cds=(1-δ)∫0Tp(0,t)S^(t)λ^(t)dt∑i=1Mp(0,ti)S^(ti).


### 3.2. Multiname Credit Derivatives

The contingent claim of a multiname credit derivative is on credit loss of the portfolio involvingK reference entities, with respective default times τ1,…,τK. The default process Nt counts the number of defaults up to time t, while the loss process Lt records the cumulative financial loss due to default up to time t, for example,
(3.7)Nt=∑k=1KI{τk≤t},Lt=∑k=1KlkI{τk≤t},
where lk is the loss due to default of the kth reference entity. Modeling the dependence of the intensity processes for the τi will be discussed in Section 3.3. The risk-neutral measure is used to price the derivative and the physical measure to predict future losses. Once the model is specified, one can use Monte Carlo simulations to evaluate the derivative price or the future expected loss.


#### 3.2.1. Basket CDS

A basket CDS is an over-the-counter derivative involving a portfolio ofn (typically between 5 and 15) reference entities. The protection buyer pays a premium, called the basket spread and expressed as a fraction of the notional principal per annum, until the kth credit event or maturity T, whichever is first; this is the premium leg of the basket CDS, similar to that of a regular CDS. The protection seller compensates the protection buyer for the credit loss associated with the kth credit event, constituting the default leg of the basket CDS. The fair basket spread equates the values of the premium and default legs, similar to a regular CDS. For k=1, the basket swap is called a first-to-default CDS. For k=n, it is called an add-up basket CDS. More generally, we can have 1≤k≤n for a kth-to-default CDS.


#### 3.2.2. Credit Indices

Indices to track CDS spreads in credit markets have been developed by different providers who reached an agreement in 2004 to consolidate them. Two index portfolios that have become standard are CDX NA IG, which is a portfolio of 125 investment-grade companies in North America, and iTraxx, which is a portfolio of 125 investment-grade companies in Europe. An index portfolio consists of CDSs referenced on all companies in the portfolio, with common maturityT, common premium dates, and common notional principal. The principal leg consists of payments by the protection buyer that are proportional to the total notional principal on the surviving companies at the premium dates. The default leg is a stream of payments by the protection seller that cover the portfolio losses up to T as they occur.


#### 3.2.3. Synthetic or Cash CDO

In a collateralized debt obligation (CDO), the sponsor chooses a portfolio of companies and a maturityT. The sponsor of a synthetic CDO sells CDS protection for each company in the portfolio with the CDS maturities equal to T. The synthetic CDO principal is the total of the notional principals underlying the CDSs. The CDO sponsor issues a collection of tranches to investors and sets different rules for determining the cash inflows and outflows for different tranches. The tranche holders sell protection to the CDO sponsor and thereby earn the specified spreads on the outstanding tranche principal, but are responsible for the sponsor's payouts to the CDS protection buyers when the companies in the portfolio default. As an illustration, suppose there are three tranches: equity, mezzanine, and senior. The equity tranche earns the largest spread of the three tranches, say 1000 basis points per year, on the outstanding tranche principal. They are responsible for all payouts on the CDS up to a certain percentage, say 5%, of the CDO principal. Thus, if the equity, mezzanine, and senior tranches have initial capital of $5 million, $15 million, and $80 million, respectively, and the defaults by companies in the portfolio lead to payouts of $2 million after 1 year, then the equity tranche holders are responsible for these payouts and the equity tranche capital reduces to $3 million.Acash CDO is created from a bond portfolio, and the sponsor issues a collection of tranches and requires an initial investment from the tranche holders to finance the underlying corporate bonds. Rules are designed to ensure that if one tranche is more senior than another it is more likely to receive the promised interest payments and repayments of principal, recognizing that there may be credit losses in these defaultable bonds.


### 3.3. Correlated Default Intensities of Multiple Obligors

For a multiname credit derivative involvingk firms, it is important to model not only the individual default intensity process but also the joint distribution of these processes. Finding tractable models that can capture the key features of the interrelationships of the firms' default intensities has been an active area of research. The following are some traditional modeling approaches described by [24].


#### 3.3.1. Gaussian Copula Approach

LetX=(X1,…,XmT) be a continuous random vector with given marginal distribution functions Fi for Xi. Since Ui=Fi(Xi) is uniformly distributed, we can specify the joint distribution of X by a copula function, which is the joint distribution function of m uniform random variables. Specifically, the copula function C(u1,…,um) defined on 0≤ui≤1  (1≤i≤m) can be used to provide a joint distribution function of X via
(3.8)C(F1(xi),…,Fm(xm))=P{F1(X1)≤F1(x1),…,Fm(Xm)≤Fm(xm)}=P{X1≤x1,…,Xm≤xm}.
Using the representation
(3.9)τi=inf{t:exp(-∫0tλi(Xs)ds)≤Ui},
in which Ui=e-ϵi is uniformly distributed on [0,1] (since ϵi is exponential with mean 1) and is independent of {Xs,s≥0}, we can use a copula function on [0,1]M to model the joint distribution of U1,…,UM and thereby also of τ1,…,τM via (3.9).LetGi be the distribution function of default time τi for the ith obligor, 1≤i≤M. Then Zi=Φ-1(Gi(τi)) is standard normal. A Gaussian copula model assumes that (Z1,…,ZM) is multivariate normal and specifies its correlation matrix Γ. Li [33] introduced it in 2000 to model default correlations, and it quickly became a widely used tool to price CDOs and other multiname credit derivatives that were previously too complex to price. A convenient way is to use the correlations of the asset returns to specify Γ. A better way is to use a factor model
(3.10)Zi=βiTX+σiϵi
to generate Γ, in which X is a multivariate normal vector of common factors affecting defaults and ϵi is standard normal independent of X. A special case is the one-factor model (2.30).


#### 3.3.2. Dependence among Default Intensities

Whereas the factor model (3.10) is applied to the Gaussian Zi=Φ-1(Gi(τi)) that involves the distribution function Gi of τi, a more direct approach is to decompose the default intensity process λti of the ith firm as λti=μt+νti, in which μt is the default intensity of a common factor and νti is that of an idiosyncratic component such that μt,νt1,…,νtM are independent processes. Instead of the preceding additive decomposition of λit, an alternative is to use a common multiplicative factor ξ so that λti=ξνti, in which ξ is a positive random variable (e.g., exponential) that is independent of the independent processes νt1,…,νtM.


#### 3.3.3. Mixture Models for Default Indicators

Vasicek [34] introduced a mixed binomial model to study the distribution of the number of defaults in a loan portfolio over a time horizon T. It played an important role in the Basel II rules and inspired subsequent work on using mixture models and factor variables; see [35, 36]. The mixed binomial model assumes that the indicator variables Xi=I{τi≤T} are conditionally independent Bernoulli (pi) given p1,…,pM and that the pi are random variables having the same distribution Π with mean p and such that Cov(pi,pj)=τ for i≠j. Then
(3.11)E(Xi)=p,Var(Xi)=p(1-p),Cov(Xi,Xj)=τfori≠j.
Therefore, the mixing distribution induces correlations among the M Bernoulli random variables that are conditionally independent given the unknown probabilities pi that have common prior distribution Π. The mean and variance of the number of defaults up to time T, denoted by #, are
(3.12)E(#)=Mp,Var(#)=Mp(1-p)+M(M-1)τ.


## 3.1. Pricing and Hedging with Credit Default Swaps

Credit default swaps (CDSs) have been described in the second paragraph of Section1. The total amount paid per year by the buyer of the CDS, as a percentage of the notional principal per annum, is called the CDS spread. Suppose the reference bond has coupon dates t1,…,tM=T and the protection buyer pays a constant premium c at these dates prior to τ, the time of the credit event. Assume that the face value of the bond is 1 and that the protection seller pays 1-δ in the case of a credit event. Then the value vb of the premium leg of the swap is
(3.1)vb=cE^{∑i=1Mexp(-∫0tirsds)I{τ>ti}}=cE^{∑i=1Mexp(-∫0ti(rs+λs)ds)}=c∑i=1Mπ(0,ti),
where rs is the short rate process, λs is the default intensity process, and π(0,ti) is the price of the defaultable bond with zero recovery and maturity ti, using the same argument as that in (2.19). Similarly the value vs of the default leg of the swap, to be paid by the protection seller, is
(3.2)vs=(1-δ)E^{exp(-∫0τrsds)I{τ≤T}}=(1-δ)∫0TE^{exp(-∫0τrsds)}f^(t)dt,
where f^ is the density function of τ under the risk-neutral measure Q, assuming that the short rate process and the default intensity process are independent under Q. Note that the survival function of τ under Q is
(3.3)S^(t)=Q(τ≥t)=E^exp(-∫0tλsds).
Letting λ^=f^/S^ be the hazard function of τ under Q, it then follows from (3.2) that
(3.4)vs=(1-δ)∫0Tp(0,t)λ^(t)S^(t)dt,
where p(0,t) is the price of the default-free bond; see (2.18). Note that under this independence assumption on the short rate and default intensity processes, we can further simplify π(0,ti) in (3.1) to
(3.5)π(0,ti)=E^[exp(-∫0tirsds)]E^[exp(-∫0tiλsds)]=p(0,ti)S^(ti).
In the absence of arbitrage, the premium cds of the default swap is the premium c for which vb=vs. Therefore, combining (3.1) with (3.4) yields
(3.6)cds=(1-δ)∫0Tp(0,t)S^(t)λ^(t)dt∑i=1Mp(0,ti)S^(ti).


## 3.2. Multiname Credit Derivatives

The contingent claim of a multiname credit derivative is on credit loss of the portfolio involvingK reference entities, with respective default times τ1,…,τK. The default process Nt counts the number of defaults up to time t, while the loss process Lt records the cumulative financial loss due to default up to time t, for example,
(3.7)Nt=∑k=1KI{τk≤t},Lt=∑k=1KlkI{τk≤t},
where lk is the loss due to default of the kth reference entity. Modeling the dependence of the intensity processes for the τi will be discussed in Section 3.3. The risk-neutral measure is used to price the derivative and the physical measure to predict future losses. Once the model is specified, one can use Monte Carlo simulations to evaluate the derivative price or the future expected loss.


### 3.2.1. Basket CDS

A basket CDS is an over-the-counter derivative involving a portfolio ofn (typically between 5 and 15) reference entities. The protection buyer pays a premium, called the basket spread and expressed as a fraction of the notional principal per annum, until the kth credit event or maturity T, whichever is first; this is the premium leg of the basket CDS, similar to that of a regular CDS. The protection seller compensates the protection buyer for the credit loss associated with the kth credit event, constituting the default leg of the basket CDS. The fair basket spread equates the values of the premium and default legs, similar to a regular CDS. For k=1, the basket swap is called a first-to-default CDS. For k=n, it is called an add-up basket CDS. More generally, we can have 1≤k≤n for a kth-to-default CDS.


### 3.2.2. Credit Indices

Indices to track CDS spreads in credit markets have been developed by different providers who reached an agreement in 2004 to consolidate them. Two index portfolios that have become standard are CDX NA IG, which is a portfolio of 125 investment-grade companies in North America, and iTraxx, which is a portfolio of 125 investment-grade companies in Europe. An index portfolio consists of CDSs referenced on all companies in the portfolio, with common maturityT, common premium dates, and common notional principal. The principal leg consists of payments by the protection buyer that are proportional to the total notional principal on the surviving companies at the premium dates. The default leg is a stream of payments by the protection seller that cover the portfolio losses up to T as they occur.


### 3.2.3. Synthetic or Cash CDO

In a collateralized debt obligation (CDO), the sponsor chooses a portfolio of companies and a maturityT. The sponsor of a synthetic CDO sells CDS protection for each company in the portfolio with the CDS maturities equal to T. The synthetic CDO principal is the total of the notional principals underlying the CDSs. The CDO sponsor issues a collection of tranches to investors and sets different rules for determining the cash inflows and outflows for different tranches. The tranche holders sell protection to the CDO sponsor and thereby earn the specified spreads on the outstanding tranche principal, but are responsible for the sponsor's payouts to the CDS protection buyers when the companies in the portfolio default. As an illustration, suppose there are three tranches: equity, mezzanine, and senior. The equity tranche earns the largest spread of the three tranches, say 1000 basis points per year, on the outstanding tranche principal. They are responsible for all payouts on the CDS up to a certain percentage, say 5%, of the CDO principal. Thus, if the equity, mezzanine, and senior tranches have initial capital of $5 million, $15 million, and $80 million, respectively, and the defaults by companies in the portfolio lead to payouts of $2 million after 1 year, then the equity tranche holders are responsible for these payouts and the equity tranche capital reduces to $3 million.Acash CDO is created from a bond portfolio, and the sponsor issues a collection of tranches and requires an initial investment from the tranche holders to finance the underlying corporate bonds. Rules are designed to ensure that if one tranche is more senior than another it is more likely to receive the promised interest payments and repayments of principal, recognizing that there may be credit losses in these defaultable bonds.


## 3.2.1. Basket CDS

A basket CDS is an over-the-counter derivative involving a portfolio ofn (typically between 5 and 15) reference entities. The protection buyer pays a premium, called the basket spread and expressed as a fraction of the notional principal per annum, until the kth credit event or maturity T, whichever is first; this is the premium leg of the basket CDS, similar to that of a regular CDS. The protection seller compensates the protection buyer for the credit loss associated with the kth credit event, constituting the default leg of the basket CDS. The fair basket spread equates the values of the premium and default legs, similar to a regular CDS. For k=1, the basket swap is called a first-to-default CDS. For k=n, it is called an add-up basket CDS. More generally, we can have 1≤k≤n for a kth-to-default CDS.


## 3.2.2. Credit Indices

Indices to track CDS spreads in credit markets have been developed by different providers who reached an agreement in 2004 to consolidate them. Two index portfolios that have become standard are CDX NA IG, which is a portfolio of 125 investment-grade companies in North America, and iTraxx, which is a portfolio of 125 investment-grade companies in Europe. An index portfolio consists of CDSs referenced on all companies in the portfolio, with common maturityT, common premium dates, and common notional principal. The principal leg consists of payments by the protection buyer that are proportional to the total notional principal on the surviving companies at the premium dates. The default leg is a stream of payments by the protection seller that cover the portfolio losses up to T as they occur.


## 3.2.3. Synthetic or Cash CDO

In a collateralized debt obligation (CDO), the sponsor chooses a portfolio of companies and a maturityT. The sponsor of a synthetic CDO sells CDS protection for each company in the portfolio with the CDS maturities equal to T. The synthetic CDO principal is the total of the notional principals underlying the CDSs. The CDO sponsor issues a collection of tranches to investors and sets different rules for determining the cash inflows and outflows for different tranches. The tranche holders sell protection to the CDO sponsor and thereby earn the specified spreads on the outstanding tranche principal, but are responsible for the sponsor's payouts to the CDS protection buyers when the companies in the portfolio default. As an illustration, suppose there are three tranches: equity, mezzanine, and senior. The equity tranche earns the largest spread of the three tranches, say 1000 basis points per year, on the outstanding tranche principal. They are responsible for all payouts on the CDS up to a certain percentage, say 5%, of the CDO principal. Thus, if the equity, mezzanine, and senior tranches have initial capital of $5 million, $15 million, and $80 million, respectively, and the defaults by companies in the portfolio lead to payouts of $2 million after 1 year, then the equity tranche holders are responsible for these payouts and the equity tranche capital reduces to $3 million.Acash CDO is created from a bond portfolio, and the sponsor issues a collection of tranches and requires an initial investment from the tranche holders to finance the underlying corporate bonds. Rules are designed to ensure that if one tranche is more senior than another it is more likely to receive the promised interest payments and repayments of principal, recognizing that there may be credit losses in these defaultable bonds.


## 3.3. Correlated Default Intensities of Multiple Obligors

For a multiname credit derivative involvingk firms, it is important to model not only the individual default intensity process but also the joint distribution of these processes. Finding tractable models that can capture the key features of the interrelationships of the firms' default intensities has been an active area of research. The following are some traditional modeling approaches described by [24].


### 3.3.1. Gaussian Copula Approach

LetX=(X1,…,XmT) be a continuous random vector with given marginal distribution functions Fi for Xi. Since Ui=Fi(Xi) is uniformly distributed, we can specify the joint distribution of X by a copula function, which is the joint distribution function of m uniform random variables. Specifically, the copula function C(u1,…,um) defined on 0≤ui≤1  (1≤i≤m) can be used to provide a joint distribution function of X via
(3.8)C(F1(xi),…,Fm(xm))=P{F1(X1)≤F1(x1),…,Fm(Xm)≤Fm(xm)}=P{X1≤x1,…,Xm≤xm}.
Using the representation
(3.9)τi=inf{t:exp(-∫0tλi(Xs)ds)≤Ui},
in which Ui=e-ϵi is uniformly distributed on [0,1] (since ϵi is exponential with mean 1) and is independent of {Xs,s≥0}, we can use a copula function on [0,1]M to model the joint distribution of U1,…,UM and thereby also of τ1,…,τM via (3.9).LetGi be the distribution function of default time τi for the ith obligor, 1≤i≤M. Then Zi=Φ-1(Gi(τi)) is standard normal. A Gaussian copula model assumes that (Z1,…,ZM) is multivariate normal and specifies its correlation matrix Γ. Li [33] introduced it in 2000 to model default correlations, and it quickly became a widely used tool to price CDOs and other multiname credit derivatives that were previously too complex to price. A convenient way is to use the correlations of the asset returns to specify Γ. A better way is to use a factor model
(3.10)Zi=βiTX+σiϵi
to generate Γ, in which X is a multivariate normal vector of common factors affecting defaults and ϵi is standard normal independent of X. A special case is the one-factor model (2.30).


### 3.3.2. Dependence among Default Intensities

Whereas the factor model (3.10) is applied to the Gaussian Zi=Φ-1(Gi(τi)) that involves the distribution function Gi of τi, a more direct approach is to decompose the default intensity process λti of the ith firm as λti=μt+νti, in which μt is the default intensity of a common factor and νti is that of an idiosyncratic component such that μt,νt1,…,νtM are independent processes. Instead of the preceding additive decomposition of λit, an alternative is to use a common multiplicative factor ξ so that λti=ξνti, in which ξ is a positive random variable (e.g., exponential) that is independent of the independent processes νt1,…,νtM.


### 3.3.3. Mixture Models for Default Indicators

Vasicek [34] introduced a mixed binomial model to study the distribution of the number of defaults in a loan portfolio over a time horizon T. It played an important role in the Basel II rules and inspired subsequent work on using mixture models and factor variables; see [35, 36]. The mixed binomial model assumes that the indicator variables Xi=I{τi≤T} are conditionally independent Bernoulli (pi) given p1,…,pM and that the pi are random variables having the same distribution Π with mean p and such that Cov(pi,pj)=τ for i≠j. Then
(3.11)E(Xi)=p,Var(Xi)=p(1-p),Cov(Xi,Xj)=τfori≠j.
Therefore, the mixing distribution induces correlations among the M Bernoulli random variables that are conditionally independent given the unknown probabilities pi that have common prior distribution Π. The mean and variance of the number of defaults up to time T, denoted by #, are
(3.12)E(#)=Mp,Var(#)=Mp(1-p)+M(M-1)τ.


## 3.3.1. Gaussian Copula Approach

LetX=(X1,…,XmT) be a continuous random vector with given marginal distribution functions Fi for Xi. Since Ui=Fi(Xi) is uniformly distributed, we can specify the joint distribution of X by a copula function, which is the joint distribution function of m uniform random variables. Specifically, the copula function C(u1,…,um) defined on 0≤ui≤1  (1≤i≤m) can be used to provide a joint distribution function of X via
(3.8)C(F1(xi),…,Fm(xm))=P{F1(X1)≤F1(x1),…,Fm(Xm)≤Fm(xm)}=P{X1≤x1,…,Xm≤xm}.
Using the representation
(3.9)τi=inf{t:exp(-∫0tλi(Xs)ds)≤Ui},
in which Ui=e-ϵi is uniformly distributed on [0,1] (since ϵi is exponential with mean 1) and is independent of {Xs,s≥0}, we can use a copula function on [0,1]M to model the joint distribution of U1,…,UM and thereby also of τ1,…,τM via (3.9).LetGi be the distribution function of default time τi for the ith obligor, 1≤i≤M. Then Zi=Φ-1(Gi(τi)) is standard normal. A Gaussian copula model assumes that (Z1,…,ZM) is multivariate normal and specifies its correlation matrix Γ. Li [33] introduced it in 2000 to model default correlations, and it quickly became a widely used tool to price CDOs and other multiname credit derivatives that were previously too complex to price. A convenient way is to use the correlations of the asset returns to specify Γ. A better way is to use a factor model
(3.10)Zi=βiTX+σiϵi
to generate Γ, in which X is a multivariate normal vector of common factors affecting defaults and ϵi is standard normal independent of X. A special case is the one-factor model (2.30).


## 3.3.2. Dependence among Default Intensities

Whereas the factor model (3.10) is applied to the Gaussian Zi=Φ-1(Gi(τi)) that involves the distribution function Gi of τi, a more direct approach is to decompose the default intensity process λti of the ith firm as λti=μt+νti, in which μt is the default intensity of a common factor and νti is that of an idiosyncratic component such that μt,νt1,…,νtM are independent processes. Instead of the preceding additive decomposition of λit, an alternative is to use a common multiplicative factor ξ so that λti=ξνti, in which ξ is a positive random variable (e.g., exponential) that is independent of the independent processes νt1,…,νtM.


## 3.3.3. Mixture Models for Default Indicators

Vasicek [34] introduced a mixed binomial model to study the distribution of the number of defaults in a loan portfolio over a time horizon T. It played an important role in the Basel II rules and inspired subsequent work on using mixture models and factor variables; see [35, 36]. The mixed binomial model assumes that the indicator variables Xi=I{τi≤T} are conditionally independent Bernoulli (pi) given p1,…,pM and that the pi are random variables having the same distribution Π with mean p and such that Cov(pi,pj)=τ for i≠j. Then
(3.11)E(Xi)=p,Var(Xi)=p(1-p),Cov(Xi,Xj)=τfori≠j.
Therefore, the mixing distribution induces correlations among the M Bernoulli random variables that are conditionally independent given the unknown probabilities pi that have common prior distribution Π. The mean and variance of the number of defaults up to time T, denoted by #, are
(3.12)E(#)=Mp,Var(#)=Mp(1-p)+M(M-1)τ.


## 4. Frailty and Contagion: HMM and Point Process Models

As noted in Section3.3, the Gaussian copula was popular because of its ease of use. However, there is no convincing argument to connect the asset return correlations to the correlations of the normally distributed transformed default times. In a 2009 commentary on “the biggest financial meltdown since the Great Depression,” Salmon [37] mentioned that the Gaussian copula approach, which “looked like an unambiguously positive breakthrough,” “was adopted by everybody from bond investors and Wall Street banks to rating agencies and regulators” and “became so deeply entrenched—and was making people so much money—that warnings about its limitations were largely ignored.” In the wake of the financial crisis, it was recognized that better albeit less tractable models of correlated default intensities are needed for pricing CDOs and risk management of credit portfolios. It was also recognized that such models should include relevant firm-level and macroeconomic variables for default prediction and also incorporate frailty and contagion mentioned in Section 1.


### 4.1. Intensity Processes with Common Frailty

The mixture binomial model for the indicator variablesI{τi≤T} with unknown means pi, allowing “random effects” in the subject-specific probabilities pi, described in Section 3.3 has more general analogs in survival analysis under the rubric of “frailty.” In the case of the proportional hazards regression model (2.13), random effects can be introduced by including a positive random factor Y in λ(t)=Yλ0(t)exp(βTx(t)), or equivalently, by including i.i.d. random effects Zi in
(4.1)λi(t)=λ0(t)exp(Zi+βTxi(t)),1≤i≤n.
Such ideas first arose in demography, in which Vaupel et al. [38] and Manton et al. [39, 40] discussed the impact of heterogeneity on the dynamics of mortality and introduced the term “frailty” for the latent random variable Y=eZ to account for the heterogeneity, which they model with a gamma distribution. Subsequently, Heckman and Singer [41, 42] and Hougaard [43, 44] considered other distributions to model frailty and applications to econometrics and biostatistics. Clayton and Cuzick [45] use the frailty approach to give an explicit formulation of the joint survival function with given marginals and to derive multivariate generalizations of survival analysis.Takingλ0≡1, Duffie et al. [46] propose an enhancement of (4.1) by including another unobserved stochastic process zθ(t), leading to the dynamic frailty model(4.2)λi(t)=exp(Zi+zθ(t)+βTxi(t)),1≤i≤n,
in which zθ(t) is an affine process, with unknown parameter vector θ, of the type used in (2.19). Because the process Zi+zθ(t) undergoing Markovian dynamics is not directly observable and the observations are (Ti,δi,xi(t)),t≤Ti,1≤i≤n, as in Section 2.1, (4.2) is a hidden Markov model (HMM) for which parameter estimation and prediction of future outcomes are much harder than the likelihood or partial likelihood methods in (2.14) or (2.15).


#### 4.1.1. State and Parameter Estimation in HMM

Duffie et al. [46] use (4.2) to model default of US public nonfinancial firms, carrying out an empirical study of 2793 companies for the period from January 1979 to March 2004. Fitting the model involves estimation of not only the parameter vectors β and θ but also the unobservable states zθ(t) and Zi. Mergers and other exits of the companies are regarded as censoring variables. Estimation of the posterior distribution of λi(t) for 0≤t≤Ti given the observations is called the “smoothing” problem in HMM, in contrast to the “filtering” problem that involves observations up to time t; see [47]. Duffie et al. [46] use the stochastic EM algorithm [48, 49] to estimate θ and β by maximum likelihood, noting that the complete data likelihood involves Zi and zθ(t) that are not observed. Replacing θ and β by their maximum likelihood estimates, they use Markov chain Monte Carlo (MCMC) to evaluate the posterior distribution of λi(t),0≤t≤Ti,i=1,…,n. The MCMC scheme they use involves both Gibbs sampling and random walk Metropolis-Hastings steps [50, 51]. Because of the computational complexity in the stochastic EM and MCMC schemes, they choose a simple Ornstein-Uhlenbeck process dzθ(t)=-θzθ(t)dt+dw(t), with a univariate parameter θ and standard Brownian motion w(t), for the affine process. The covariate vector xi(t) in (4.2) is the same as used earlier by [52] without the frailty terms Zi and zθ(t). It includes the following firm-specific covariates xij(t) and macroeconomic covariates uk(t):(i)
the firm's distance to default,(ii)
the firm's trailing 1-year stock return,(iii)
the 3-month Treasury bill rate, which is a macroeconomic variable “consistent with the effect of a monetary policy that lowers short-term interest rates when the economy is likely to perform poorly,”(iv)
the trailing 1-year return on the S&P 500 index.


### 4.2. Contagion and Point Process Models

Davis and Lo [53] introduced a simple model to reflect default contagion among n homogeneous firms. Suppose prior to any default all firms have default intensity λ. After a default occurs, the default intensities of the surviving firms increase by a common amount μ. This is an example of the more general “top-down” approach to modeling portfolio loss and contagion. Consider a bond portfolio involving n companies and order their default times by τ(1)≤τ(2)≤⋯≤τ(n). The number of defaults up to time t is the point process Nt=∑i=1n‍I{τ(i)≤t}. The top-down approach models directly the intensity of the point process Nt. In contrast, the “bottom-up” approach models the intensities λi(t) of the individual default times τi and aggregates them into the point process Nt with intensity ∑i=1n‍λi(t)I{τi≥t}. One such contagion model is λi(t)=λi0(t)+βiNt, where Nt=(I{τ1≤t}⋯I{τn≤t})T, βij measures the increase in hazard for firm i after firm j defaults, and λi0 is the Cox-Ingersoll-Ross process for default intensity that is given by
(4.3)dλi0(t)=κi(θi-λi0(t))dt+σiλi0(t)dwi(t),
in which wi(t),t≥0, is standard Brownian motion; see [54].Additive regression models for hazard functions are often used to model contagion in biostatistics, generalizing the simple contagion model of [53] described above. Aalen [55] and Huffer and McKeague [56] provide the statistical counterparts of the multiplicative hazard model (2.13).


## 4.1. Intensity Processes with Common Frailty

The mixture binomial model for the indicator variablesI{τi≤T} with unknown means pi, allowing “random effects” in the subject-specific probabilities pi, described in Section 3.3 has more general analogs in survival analysis under the rubric of “frailty.” In the case of the proportional hazards regression model (2.13), random effects can be introduced by including a positive random factor Y in λ(t)=Yλ0(t)exp(βTx(t)), or equivalently, by including i.i.d. random effects Zi in
(4.1)λi(t)=λ0(t)exp(Zi+βTxi(t)),1≤i≤n.
Such ideas first arose in demography, in which Vaupel et al. [38] and Manton et al. [39, 40] discussed the impact of heterogeneity on the dynamics of mortality and introduced the term “frailty” for the latent random variable Y=eZ to account for the heterogeneity, which they model with a gamma distribution. Subsequently, Heckman and Singer [41, 42] and Hougaard [43, 44] considered other distributions to model frailty and applications to econometrics and biostatistics. Clayton and Cuzick [45] use the frailty approach to give an explicit formulation of the joint survival function with given marginals and to derive multivariate generalizations of survival analysis.Takingλ0≡1, Duffie et al. [46] propose an enhancement of (4.1) by including another unobserved stochastic process zθ(t), leading to the dynamic frailty model(4.2)λi(t)=exp(Zi+zθ(t)+βTxi(t)),1≤i≤n,
in which zθ(t) is an affine process, with unknown parameter vector θ, of the type used in (2.19). Because the process Zi+zθ(t) undergoing Markovian dynamics is not directly observable and the observations are (Ti,δi,xi(t)),t≤Ti,1≤i≤n, as in Section 2.1, (4.2) is a hidden Markov model (HMM) for which parameter estimation and prediction of future outcomes are much harder than the likelihood or partial likelihood methods in (2.14) or (2.15).


### 4.1.1. State and Parameter Estimation in HMM

Duffie et al. [46] use (4.2) to model default of US public nonfinancial firms, carrying out an empirical study of 2793 companies for the period from January 1979 to March 2004. Fitting the model involves estimation of not only the parameter vectors β and θ but also the unobservable states zθ(t) and Zi. Mergers and other exits of the companies are regarded as censoring variables. Estimation of the posterior distribution of λi(t) for 0≤t≤Ti given the observations is called the “smoothing” problem in HMM, in contrast to the “filtering” problem that involves observations up to time t; see [47]. Duffie et al. [46] use the stochastic EM algorithm [48, 49] to estimate θ and β by maximum likelihood, noting that the complete data likelihood involves Zi and zθ(t) that are not observed. Replacing θ and β by their maximum likelihood estimates, they use Markov chain Monte Carlo (MCMC) to evaluate the posterior distribution of λi(t),0≤t≤Ti,i=1,…,n. The MCMC scheme they use involves both Gibbs sampling and random walk Metropolis-Hastings steps [50, 51]. Because of the computational complexity in the stochastic EM and MCMC schemes, they choose a simple Ornstein-Uhlenbeck process dzθ(t)=-θzθ(t)dt+dw(t), with a univariate parameter θ and standard Brownian motion w(t), for the affine process. The covariate vector xi(t) in (4.2) is the same as used earlier by [52] without the frailty terms Zi and zθ(t). It includes the following firm-specific covariates xij(t) and macroeconomic covariates uk(t):(i)
the firm's distance to default,(ii)
the firm's trailing 1-year stock return,(iii)
the 3-month Treasury bill rate, which is a macroeconomic variable “consistent with the effect of a monetary policy that lowers short-term interest rates when the economy is likely to perform poorly,”(iv)
the trailing 1-year return on the S&P 500 index.


## 4.1.1. State and Parameter Estimation in HMM

Duffie et al. [46] use (4.2) to model default of US public nonfinancial firms, carrying out an empirical study of 2793 companies for the period from January 1979 to March 2004. Fitting the model involves estimation of not only the parameter vectors β and θ but also the unobservable states zθ(t) and Zi. Mergers and other exits of the companies are regarded as censoring variables. Estimation of the posterior distribution of λi(t) for 0≤t≤Ti given the observations is called the “smoothing” problem in HMM, in contrast to the “filtering” problem that involves observations up to time t; see [47]. Duffie et al. [46] use the stochastic EM algorithm [48, 49] to estimate θ and β by maximum likelihood, noting that the complete data likelihood involves Zi and zθ(t) that are not observed. Replacing θ and β by their maximum likelihood estimates, they use Markov chain Monte Carlo (MCMC) to evaluate the posterior distribution of λi(t),0≤t≤Ti,i=1,…,n. The MCMC scheme they use involves both Gibbs sampling and random walk Metropolis-Hastings steps [50, 51]. Because of the computational complexity in the stochastic EM and MCMC schemes, they choose a simple Ornstein-Uhlenbeck process dzθ(t)=-θzθ(t)dt+dw(t), with a univariate parameter θ and standard Brownian motion w(t), for the affine process. The covariate vector xi(t) in (4.2) is the same as used earlier by [52] without the frailty terms Zi and zθ(t). It includes the following firm-specific covariates xij(t) and macroeconomic covariates uk(t):(i)
the firm's distance to default,(ii)
the firm's trailing 1-year stock return,(iii)
the 3-month Treasury bill rate, which is a macroeconomic variable “consistent with the effect of a monetary policy that lowers short-term interest rates when the economy is likely to perform poorly,”(iv)
the trailing 1-year return on the S&P 500 index.


## 4.2. Contagion and Point Process Models

Davis and Lo [53] introduced a simple model to reflect default contagion among n homogeneous firms. Suppose prior to any default all firms have default intensity λ. After a default occurs, the default intensities of the surviving firms increase by a common amount μ. This is an example of the more general “top-down” approach to modeling portfolio loss and contagion. Consider a bond portfolio involving n companies and order their default times by τ(1)≤τ(2)≤⋯≤τ(n). The number of defaults up to time t is the point process Nt=∑i=1n‍I{τ(i)≤t}. The top-down approach models directly the intensity of the point process Nt. In contrast, the “bottom-up” approach models the intensities λi(t) of the individual default times τi and aggregates them into the point process Nt with intensity ∑i=1n‍λi(t)I{τi≥t}. One such contagion model is λi(t)=λi0(t)+βiNt, where Nt=(I{τ1≤t}⋯I{τn≤t})T, βij measures the increase in hazard for firm i after firm j defaults, and λi0 is the Cox-Ingersoll-Ross process for default intensity that is given by
(4.3)dλi0(t)=κi(θi-λi0(t))dt+σiλi0(t)dwi(t),
in which wi(t),t≥0, is standard Brownian motion; see [54].Additive regression models for hazard functions are often used to model contagion in biostatistics, generalizing the simple contagion model of [53] described above. Aalen [55] and Huffer and McKeague [56] provide the statistical counterparts of the multiplicative hazard model (2.13).


## 5. Insurance Contracts, Credibility Theory, and Linear EB

Credibility theory for pricing insurance contracts in actuarial science uses historical data from an individual contract and the risk class to which it belongs to arrive at acredibility premium of the form
(5.1)αY-+(1-α)μ,with0<α<1,
where Y- is the observed mean claim amount for the individual contract and μ is the overall mean claim per contract in the class during the previous period; see [57]. It can be derived as a Bayes estimate of θ under the assumption that Y- is normal (resp., gamma, binomial, or Poisson) with mean θ and that θ has a prior normal (resp., inverse gamma, beta, or gamma) distribution with mean μ, for which the form of α is explicitly determined by the distributional parameters. Without assuming a prior distribution with specified parameters, Bühlmann [58, 59] considered the best linear estimator of the form αY-+b based on observed claims Y1,…,Yn. The best linear estimator involves the mean and variance of the Yi, which can be estimated consistently by the method of moments, as will be explained below.


### 5.1. Linear Empirical Bayes and Credibility Factors

The empirical Bayes (EB) methodology introduced by Robbins [1] for n independent and structurally similar problems of statistical inference on unknown parameters θi based on observed data Yi (i=1,…,n), where Yi has probability density f(y∣θi) and θi and Yi may be vectors. The θi are assumed to have a common prior distribution G that has unspecified hyperparameters. Letting dG(y) be the Bayes decision rule (with respect to some loss function and assuming known hyperparameters) when Yi=y is observed, the basic principle underlying empirical Bayes is that dG can often be consistently estimated from Y1,…,Yn, leading to the empirical Bayes decision rule dG^. Thus, the n structurally similar problems can be pooled to provide information about unspecified hyperparameters in the prior distribution, thereby yielding G^ and the decision rules dG^(Yi) for the independent problems. In particular, Robbins [1] considered Poisson Yi with mean θi, as in the case of the number of accidents by the ith driver in a sample of size n (in a given year) from a population of drivers, with distribution G for the accident-proneness parameter θ. In this case the Bayes estimate (with respect to squared error loss) of θi when Yi=y is observed is
(5.2)dg(y)=(y+1)g(y+1)g(y),y=0,1,…,
where g(y)=∫0∞‍θyeθdG(θ)/(y!). Using g^(k)=n-1∑i=1n‍I{Yi=k} to replace g(k) in (5.2) yields the (EB) estimate dg^(y). Stein [2] and subsequently James and Stein [60] considered the special case Yi~N(θi,σ2) with known σ but unknown θi.Unlike the pricing of credit derivatives, credibility theory for pricing insurance contracts uses the physical (real-world) measure rather than the risk neutral measure. This extends to insurance derivatives such as CAT (catastrophic) bonds used by the insurance industry to hedge exposures to catastrophic risks due to hurricanes, earthquakes, or other severe unexpected (rare) events. For example, an insurance company that wants protection for California earthquake losses between $30 million and $40 million can issue CAT bonds with $10 million total principal. In the event that the losses exceed $30 million, bondholders would lose some or all of their principal, or only their promised interest if the insurance company has made a bond issue large enough to cover this excess cost layer. Like insurance contracts, insurance derivatives are priced under the physical measure, using historical or actuarial data to estimate the parameters of future earthquake or other catastrophic losses. The value of the estimated expected payoff is discounted at the risk-free rate to price an insurance derivative; the uncertainties of their payoffs are not priced by the market.Robbins [61] gives a general description of linear EB models, which include Stein's normal model as a special case and which is also related to Bühlmann's linear estimation approach to credibility premium. Suppose the decision rule d(y) is restricted to linear functions of the form A+By and assume that the family of density functions f(y∣θ) satisfies
(5.3)E(Y∣θ)=θ,Var(Y∣θ)=a+bθ+cθ2
for some constants a, b, and c≠-1. Then the Bayes rule dG(y) that minimizes the Bayes risk, with respect to squared error loss, among linear estimators of θ is given by the linear regression formula
(5.4)dG(y)=Eθ+Cov(θ,Y)Var(Y)(y-EY).
Moreover, from (5.3), it follows that
(5.5)EY=Eθ,Cov(θ,Y)=Var(θ),(5.6)Var(Y)=E[Var(Y∣θ)]+Var[E(Y∣θ)]=a+bEθ+c(Eθ)2+(c+1)Var(θ),
and therefore
(5.7)Var(θ)=Var(Y)-{a+bEY+c(EY)2}c+1.
Since (θ1,Y1),…,(θn,Yn) are i.i.d. vectors from the distribution under which θi has prior distribution G and Yi has conditional density f(·∣θi) given Yi, the parameters Eθ, Cov(θ,Y), Var(Y), and EY in the Bayes rule (5.4) can be estimated consistently by using (5.5), (5.7), and the method of moments. This yields the linear EB estimate
(5.8)d^(y)=Y-+[1-cs2+a+bY-+cY-2(c+1)s2]+(y-Y-),
where x+=max(x,0), Y-=n-1∑i=1n‍Yi estimates the mean EY, and s2=∑i=1n‍(Yi-Y-)2/(n-1) estimates the variance Var(Y). In the case Yi∣θi~N(θi,σ2) with known σ2, (5.3) holds with a=σ2 and b=c=0, for which (5.8) corresponds to a variant of the James-Stein [60] estimator of (θ1,…,θn)T which has a smaller mean squared error than the maximum likelihood estimator (Y1,…,Yn)T. Moreover, if G is normal, then (5.4) minimizes the Bayes risk among all (not necessarily linear) estimators of θ.The linear EB estimate (5.8) was used by Bühlmann [58] to determine the credibility factorα in (5.1), writing α=αn to indicate its dependence on the training sample Y1,…,Yn of size n. The monograph by Bühlmann and Gisler [62] describes a variety of extensions of (5.8) to more general settings. For example, suppose there are I risk classes and Yij denotes the jth claim of the ith class, with E(Yij∣θi)=θi and Var(Yij∣θi)=σi2, 1≤j≤ni,1≤i≤I. Assuming a normal prior N(μ,τ2) for θi, the Bayes estimate of θi is
(5.9)E(θi∣Yi1,…,Yi,ni)=αiY-i+(1-αi)μ,
where αi=τ2/(τ2+σi2/ni) and Y-i=ni-1∑j=1ni‍Yij. Since E(Yij)=E(θi)=μ and Var(Yij)=Var(θi)+E(Var(Yij∣θi))=τ2+σi2, we can estimate μ,σi2, and τ2 by the method of moments when ni>1: μ^=(∑i=1I‍∑j=1ni‍Yij)/∑i=1I‍ni, σ^i2=∑j=1ni‍(Yij-Y-i)2/(ni-1), and τ^2=∑i=1I‍ni(Y-i-μ^)2/∑i=1I‍ni. Replacing μ, σi2, and τ2 in (5.9) by μ^, σ^i2, and τ^2 yields the EB estimate
(5.10)θ^i=α^iY-i+(1-α^i)μ^,
where α^i=τ^2/(τ^2+σ^i2/ni) is the credibility factor for the ith class. Another important extension, introduced by Hachemeister [63], is the credibility regression model that relates claim sizes to certain covariates; the credibility factor is a matrix in this case.Frees et al. [64, 65] unified various credibility models into the framework of linear mixed models (LMM) of the form
(5.11)Yij=βTxij+biTzij+ϵij,
with fixed effects β1,…,βp forming the vector β, subject-specific random effects bi1,…,biq forming the vector bi such that E(bi)=0, and zero-mean random disturbances ϵij that have variance σ2 and are uncorrelated with the random effects and the covariates xij and zij. In particular, the credibility regression model Yi=αiTxi+ϵi, in which the subject-specific regression parameters αi have distribution G, is a special case with j=1=ni, zi=xi, β being the mean of G, and bi=αi-β.


## 5.1. Linear Empirical Bayes and Credibility Factors

The empirical Bayes (EB) methodology introduced by Robbins [1] for n independent and structurally similar problems of statistical inference on unknown parameters θi based on observed data Yi (i=1,…,n), where Yi has probability density f(y∣θi) and θi and Yi may be vectors. The θi are assumed to have a common prior distribution G that has unspecified hyperparameters. Letting dG(y) be the Bayes decision rule (with respect to some loss function and assuming known hyperparameters) when Yi=y is observed, the basic principle underlying empirical Bayes is that dG can often be consistently estimated from Y1,…,Yn, leading to the empirical Bayes decision rule dG^. Thus, the n structurally similar problems can be pooled to provide information about unspecified hyperparameters in the prior distribution, thereby yielding G^ and the decision rules dG^(Yi) for the independent problems. In particular, Robbins [1] considered Poisson Yi with mean θi, as in the case of the number of accidents by the ith driver in a sample of size n (in a given year) from a population of drivers, with distribution G for the accident-proneness parameter θ. In this case the Bayes estimate (with respect to squared error loss) of θi when Yi=y is observed is
(5.2)dg(y)=(y+1)g(y+1)g(y),y=0,1,…,
where g(y)=∫0∞‍θyeθdG(θ)/(y!). Using g^(k)=n-1∑i=1n‍I{Yi=k} to replace g(k) in (5.2) yields the (EB) estimate dg^(y). Stein [2] and subsequently James and Stein [60] considered the special case Yi~N(θi,σ2) with known σ but unknown θi.Unlike the pricing of credit derivatives, credibility theory for pricing insurance contracts uses the physical (real-world) measure rather than the risk neutral measure. This extends to insurance derivatives such as CAT (catastrophic) bonds used by the insurance industry to hedge exposures to catastrophic risks due to hurricanes, earthquakes, or other severe unexpected (rare) events. For example, an insurance company that wants protection for California earthquake losses between $30 million and $40 million can issue CAT bonds with $10 million total principal. In the event that the losses exceed $30 million, bondholders would lose some or all of their principal, or only their promised interest if the insurance company has made a bond issue large enough to cover this excess cost layer. Like insurance contracts, insurance derivatives are priced under the physical measure, using historical or actuarial data to estimate the parameters of future earthquake or other catastrophic losses. The value of the estimated expected payoff is discounted at the risk-free rate to price an insurance derivative; the uncertainties of their payoffs are not priced by the market.Robbins [61] gives a general description of linear EB models, which include Stein's normal model as a special case and which is also related to Bühlmann's linear estimation approach to credibility premium. Suppose the decision rule d(y) is restricted to linear functions of the form A+By and assume that the family of density functions f(y∣θ) satisfies
(5.3)E(Y∣θ)=θ,Var(Y∣θ)=a+bθ+cθ2
for some constants a, b, and c≠-1. Then the Bayes rule dG(y) that minimizes the Bayes risk, with respect to squared error loss, among linear estimators of θ is given by the linear regression formula
(5.4)dG(y)=Eθ+Cov(θ,Y)Var(Y)(y-EY).
Moreover, from (5.3), it follows that
(5.5)EY=Eθ,Cov(θ,Y)=Var(θ),(5.6)Var(Y)=E[Var(Y∣θ)]+Var[E(Y∣θ)]=a+bEθ+c(Eθ)2+(c+1)Var(θ),
and therefore
(5.7)Var(θ)=Var(Y)-{a+bEY+c(EY)2}c+1.
Since (θ1,Y1),…,(θn,Yn) are i.i.d. vectors from the distribution under which θi has prior distribution G and Yi has conditional density f(·∣θi) given Yi, the parameters Eθ, Cov(θ,Y), Var(Y), and EY in the Bayes rule (5.4) can be estimated consistently by using (5.5), (5.7), and the method of moments. This yields the linear EB estimate
(5.8)d^(y)=Y-+[1-cs2+a+bY-+cY-2(c+1)s2]+(y-Y-),
where x+=max(x,0), Y-=n-1∑i=1n‍Yi estimates the mean EY, and s2=∑i=1n‍(Yi-Y-)2/(n-1) estimates the variance Var(Y). In the case Yi∣θi~N(θi,σ2) with known σ2, (5.3) holds with a=σ2 and b=c=0, for which (5.8) corresponds to a variant of the James-Stein [60] estimator of (θ1,…,θn)T which has a smaller mean squared error than the maximum likelihood estimator (Y1,…,Yn)T. Moreover, if G is normal, then (5.4) minimizes the Bayes risk among all (not necessarily linear) estimators of θ.The linear EB estimate (5.8) was used by Bühlmann [58] to determine the credibility factorα in (5.1), writing α=αn to indicate its dependence on the training sample Y1,…,Yn of size n. The monograph by Bühlmann and Gisler [62] describes a variety of extensions of (5.8) to more general settings. For example, suppose there are I risk classes and Yij denotes the jth claim of the ith class, with E(Yij∣θi)=θi and Var(Yij∣θi)=σi2, 1≤j≤ni,1≤i≤I. Assuming a normal prior N(μ,τ2) for θi, the Bayes estimate of θi is
(5.9)E(θi∣Yi1,…,Yi,ni)=αiY-i+(1-αi)μ,
where αi=τ2/(τ2+σi2/ni) and Y-i=ni-1∑j=1ni‍Yij. Since E(Yij)=E(θi)=μ and Var(Yij)=Var(θi)+E(Var(Yij∣θi))=τ2+σi2, we can estimate μ,σi2, and τ2 by the method of moments when ni>1: μ^=(∑i=1I‍∑j=1ni‍Yij)/∑i=1I‍ni, σ^i2=∑j=1ni‍(Yij-Y-i)2/(ni-1), and τ^2=∑i=1I‍ni(Y-i-μ^)2/∑i=1I‍ni. Replacing μ, σi2, and τ2 in (5.9) by μ^, σ^i2, and τ^2 yields the EB estimate
(5.10)θ^i=α^iY-i+(1-α^i)μ^,
where α^i=τ^2/(τ^2+σ^i2/ni) is the credibility factor for the ith class. Another important extension, introduced by Hachemeister [63], is the credibility regression model that relates claim sizes to certain covariates; the credibility factor is a matrix in this case.Frees et al. [64, 65] unified various credibility models into the framework of linear mixed models (LMM) of the form
(5.11)Yij=βTxij+biTzij+ϵij,
with fixed effects β1,…,βp forming the vector β, subject-specific random effects bi1,…,biq forming the vector bi such that E(bi)=0, and zero-mean random disturbances ϵij that have variance σ2 and are uncorrelated with the random effects and the covariates xij and zij. In particular, the credibility regression model Yi=αiTxi+ϵi, in which the subject-specific regression parameters αi have distribution G, is a special case with j=1=ni, zi=xi, β being the mean of G, and bi=αi-β.


## 6. Dynamic Empirical Bayes and Generalized Linear Mixed Models

Linear credibility theory, with a credibility premium of the form (5.1), has been extended in two directions. The first is linear evolutionary credibility that assumes a dynamic linear Bayesian model for the prior means over time; see [62]. The second is exact credibility theory for which the claims are assumed to belong to an exponential family of distributions so that generalized linear (or generalized linear mixed) models are used in lieu of linear regression (or linear mixed) models; see [66–72] that assume Yij to have a density function of the form
(6.1)f(y∣bi,zit,xit)=exp{[yθit-ψ(θit)]σ+c(y,σ)},
in which σ is a dispersion parameter and μit=dψ/dθ|θ=θit satisfies μit=g-1(βTxit) in the case of generalized linear models and
(6.2)μit=g-1(βTxit+biTzit)
in the case of generalized linear mixed models (GLMM), where g-1 is the inverse of a monotone link function g. The case g=dψ/dθ is called the “canonical link.” The random effects bi can contain an intercept term ai by augmenting the covariate vector to (1,zitT) in case ai is not included in bi; β is a vector of fixed effects and can likewise contain an intercept term. The density function (6.1) with σ=1 is that of an exponential family, which includes the Bernoulli and normal distributions as special cases.


### 6.1. Linear Evolutionary Credibility, Kalman Filtering, and Dynamic EB

Noting that insurers observe claims of risk classes over successive periods, Frees et al. [64, 65] in fact consider the setting of longitudinal data in their LMM approach to credibility theory; the index j in (5.11) is replaced by t that denotes time so that the Yit represent longitudinal data of claims. Bühlmann and Gisler [62] go further in assuming a dynamic Bayesian model in which the mean θit of Yit has a prior distribution with mean μt that changes with time according to the AR(1) model
(6.3)μt=ρμt-1+(1-ρ)μ+ηt,
with i.i.d. zero-mean ηt having variance V. Thus, the μt are latent states undergoing AR(1) dynamics, yielding a linear state-space model for which the unobserved states μt can be estimated from the observations Yis,s≤t, by the Kalman filter μ^t|t given by the recursion
(6.4)μ^t∣t=μ^t∣t-1+ρ-1Kt(Yt-μ^t∣t-11),μ^t+1∣t=ρμ^t∣t+(1-ρ)μ,
in which Yt=(Y1t,…,Ynt)T, 1=(1,…,1)T, and Kt is the 1×n Kalman gain matrix defined recursively in terms of the hyperparameters V=Var(ηt), vt=Var(Yit∣μt), and ρ by
(6.5)Kt=ρΣt∣t-11T(Σt∣t-111T+vtI)-1,(6.6)Σt+1∣t=(ρ-Kt1)2Σt∣t-1+V+vtKtKtT.
Under the assumption that θit∣μt has variance τ2 and Yit∣θit has variance a+bθit+cθit2 as in (5.3), vt=τ2+a+bθit+cθit2. The Kalman filter is the minimum-variance linear estimator of μt; it is also the Bayes estimator if Yit (conditional on μt) and ηt are normal.As in the linear EB approach described in Section5, we can use the method of moments to estimate μt and vt for every 1≤t≤T. Specifically, Y-t=n-1∑i=1n‍Yit is a consistent estimate of μt and st2=∑i=1n‍(Yit-Y-t)2/(n-1) is a consistent estimate of vt. Moreover, assuming E(μ0)=μ, it follows from (6.3) that E(μt)=μ for all t. Therefore, μ can be consistently estimated from the observations up to time t by μ^(t)=(∑s=1t‍Y-s)/t, and ρ and V can be consistently estimated by
(6.7)ρ^(t)=∑s=2t(Y-s-μ^(t))(Y-s-1-μ^(t))∑s=2t(Y-s-μ^(t))2,V^(t)=∑s=2t{Y-s-ρY-s-1-(1-ρ)μ^(t)}2(t-2),
for t≥2. Replacing the hyperparameters μ, ρ, V, and vt in the Bayes estimate μ^t|t of μt by their estimates μ^(t), ρ^(t), V^(t), and st2 yields the EB estimate of μt. Of particular interest in applications of longitudinal data is the prediction of Yi,t+1 based on the set of observations 𝒴t up to time t. The Bayes predictor, assuming normal priors and normal observations, is
(6.8)E(Yi,t+1∣𝒴t)=E(θi,t+1∣𝒴t)=E(μt+1∣𝒴t)=μ^t+1∣t,
which is the Kalman predictor given by (6.4). This is also the best linear predictor without assuming normality. The corresponding EB predictor replaces the hyperparameters in the Kalman predictor by their method-of-moment estimates μ^(t), ρ^(t), V^(t), and st2.Note that to estimate the hyperparameters in this EB approach via (6.3) and Kalman filtering, we have used the cross-sectional mean Y-t-1 of n independent observations that have mean μt-1. Lai and Sun [73] propose an alternative and more direct approach that replaces μt-1 in (6.3) by Y-t-1, leading to
(6.9)μt=ρY-t-1+ω+ηt
with ω=(1-ρ)μ, which can be written as an LMM for the longitudinal data Yij:
(6.10)Yit=ρY-t-1+ω+bi+ϵit,
by absorbing ηt into ϵit. Since (6.10) is in the form of a regression model, Lai and Sun [73] also include additional covariates to increase the predictive power of the model in the LMM
(6.11)Yit=ρY-t-1+ai+βTxit+biTzit+ϵit,
where ai and bi are subject-specific random effects, xit represents a vector of subject-specific covariates that are available prior to time t (for predicting Yit prior to observing it at time t), and zit denotes a vector of additional covariates that are associated with the random effects bi. A simulation study in [71] to compare the performance of parameter estimation and one-step-ahead prediction using (6.11) with that of Kalman filtering shows that the dynamic LMM (6.11) yields results comparable to those of Kalman filtering when the data are generated by a linear state-space model of the type (6.5), the parameters of which are not assumed to be known and are estimated by (6.7).


### 6.2. Dynamic EB Modeling via GLMM

Lai and Sun [73] assume the prior distribution for which the μit are i.i.d. with mean μt and
(6.12)g(μt)=∑j=1pθjg(Y-t-j)for1≤t≤T,whereY-s=n-1∑i=1nYis.
The dynamic model (6.12) for μt is an EB version of the Markov model introduced by Zeger and Qaqish [74], who consider a single time series Yt that has a density function of the form (6.1) and models μt by g(μt)=∑j=1p‍θjg(Yt-j), where g is a link function and n=1 in this case. For general n, note that μs is the mean of μis and can be consistently estimated by Y-s, which is the basic idea underlying the linear EB approach. This suggests using g(Y-t-j) in lieu of g(Yt-j) for n≥1 that leads to (6.12). As in the LMM (6.11), Lai and Sun [73] propose to increase the predictive power of the model by including fixed and random effects and other time-varying covariates of each subject i, thereby removing the dependence of g(μit)-g(μt) on t in the GLMM
(6.13)g(μit)=∑j=1pθjg(Y-t-j)+ai+βTxit+biTzit,
in which θ1,…,θp and β are the fixed effects and ai and bi are subject-specific random effects. Following [4], it is assumed that ai and bi are independent normal with zero means. Denote the unknown vector specifying the covariances of the random effects by α, and for notational simplicity, we can augment bi to include ai so that (6.13) can be written as g(μit)=∑j=1p‍θjg(Y-t-j)+xitTβ+(1,zitT)bi. Note that (6.13) is a natural extension of the dynamic linear EB model (6.11) to the exponential family. By including g(Y-t-j),j=1,…,p, in xit, we can implement the dynamic EB model (6.13) by using the methods described for general GLMMs.


### 6.3. Implementation and Model Selection

Breslow and Clayton assume thebi in (6.2) to have a common normal distribution with mean 0 and covariance matrix Σ that depends on an unknown parameter. Lai and Shih [75] have developed a nonparametric method to estimate the random effects distribution in the more general setting of nonlinear mixed effects models and have shown that there is very low resolution in estimating the distribution. This provides a heuristic explanation for why the choice of a normal distribution, with unspecified parameters, for the random effects distribution G does not result in a worse estimate of β and σ than the semiparametric approach that estimates G nonparametrically. An advantage of assuming a normal distribution for the random effects bi is that it only involves the covariance matrix of the random effects. The likelihood function of the GLMM defined by (6.1) and (6.2) is of the form ∏i=1n‍Li(σ,α,β), where
(6.14)Li(σ,α,β)=∫{∏t=1Tf(yit;θit,σ)}ϕα(b)db,
in which ϕα denotes the normal density function with mean 0 and covariance matrix depending on an unknown parameter α. We review three methods to compute the likelihood function, the maximizer of which gives the MLE of σ, α, and β.


#### 6.3.1. Laplace's Approximation

Lettingeli(b∣σ,α,β) be the integrand in the right-hand side of (6.14), Laplaces's asymptotic formula for integrals yields the approximation
(6.15)∫eli(b∣σ,α,β)db≈(2π)q/2{det[-l¨i(b^i∣σ,α,β)]}-1/2exp{li(b^i∣σ,α,β)},
where q is the dimension of bi, b^i=b^i(σ,α,β) is the maximizer of li(bσ,α,β), and l¨i denotes the Hessian matrix consisting of second partial derivatives of li with respect to the components of b. The R package lme4 computes the MLE by using the Laplace approximation to (6.14) or the restricted pseudolikelihood approach proposed by [76], as the user-specified option.


#### 6.3.2. Gauss-Hermite Quadrature

Laplaces's asymptotic formula (6.15) is derived from the asymptotic approximation of li by a quadratic function of b in a small neighborhood of b^i as λmin(-l¨i(b^i∣σ,α,β))→∞, where λmin(·) denotes the minimum eigenvalue of a symmetric matrix. Therefore, such formula may produce significant approximation error for Li if the corresponding λmin(-l¨i(b^i∣σ,α,β)) is not sufficiently large. One way to reduce the possible approximation error is to compute Li by using an adaptive Gauss-Hermite quadrature rule, as in [77]. The software package SAS uses adaptive Gauss-Hermite quadrature in the NLMIXED procedure to compute (6.14); the R function lmer(·) also uses Gaussian quadrature to compute (6.14) but only for certain special cases of the exponential family (6.1).


#### 6.3.3. Hybrid Method

The numerical integration procedure demands a much higher computational effort and becomes computationally infeasible ifn or q is large. To circumvent the issue of high-dimensional numerical integration, Antonio and Beirlant [72] propose to put prior distributions on the unknown parameters and to estimate them by using the Markov chain Monte Carlo (MCMC) method in a Bayesian way. The performance of the MCMC method, however, depends on how the prior parameters are set as well as the convergence rate of the Markov chain to its stationary distribution, which may not even exist. Yafune et al. [78] use direct Monte Carlo integration but point out that the computational time may be too long to be of practical interest. Instead of relying solely on Monte Carlo, Lai and Shih [79] proposed the following hybrid approach. Let Vi=-l¨i(b^i|σ,α,β). Since Laplaces's asymptotic formula (6.15) may be a poor approximation to (6.14) when λmin(Vi) is not sufficiently large, Monte Carlo integration, whose error is independent of q, can be used as an alternative method to evaluate (6.14).Lai and Shih [79] suggest using Monte Carlo integration instead of Laplace's asymptotic formula for i whose λmin(Vi)<c, where c is a positive threshold. Lai et al. [80, 81] further improve the hybrid method by using importance sampling in the Monte Carlo component of the hybrid method. Specifically, instead of sampling b(h) directly from ϕα as in [79], sample it from a mixture of the prior normal distribution with density ϕα and the posterior normal distribution N(b^i,[-l¨i(b^i∣σ,α,β)+εI]-1), where ε is a small positive number to ensure that the covariance matrix is invertible. We denote such choice of trial density as λ. It has the advantage of further incorporating the essence of Laplace's method in the Monte Carlo step such that the method is less dependent (than direct Monte Carlo) on the choice of the threshold c.(i)
Ifλmin(Vi)<c, evaluate Li(σ,α,β) in (6.14) by the Monte Carlo approximation
(6.16){∑h=1H[wh∏t=1Tf(yit;θit(b(h),β),σ)]}{∑h=1Hwh}
in which b(1),…,b(H) are independently sampled from the importance density λ and wh=ϕα(b(h))/λ(b(h)) are the importance weights.(ii)
Ifλmin(Vi)≥c, use Laplace approximation (6.15). Lai et al. [80] propose to use c=10,H=100, and a mixture distribution that assigns 0.2 probability to ϕα. Note that choosing λ=ϕα yields wh=1, for which (6.16) reduces to direct Monte Carlo evaluation of Li(σ,α,β) in [79]. The likelihoods Li(σ,α,β) thus computed then yield the likelihood function ∏i=1n‍Li(σ,α,β), whose logarithm can be maximized by using iterative gradient-type schemes.


#### 6.3.4. Model Selection

As shown by Breslow and Clayton [4], a consistent estimator of the asymptotic covariance matrix of (σ^,α^,β^) is V^-1={-∑i=1n‍R¨i(σ^,α^,β^)}-1, where Ri=logLi and R¨i denotes the Hessian matrix of second derivatives of Ri with respect to σ, α, and the components of β. As in [79], we compute the observed information matrix V^=-∑i=1n‍R¨i(σ^,α^,β^) by taking numerical derivatives and using the hybrid method to evaluate Li. Lai and Shih [79] have shown that V^-1/2(σ^-σ,α-α^,β-β^) has a limiting standard multivariate normal distribution as n→∞ under certain regularity conditions, and that the sandwich estimator
(6.17){-∑i=1nR¨i(σ^,α^,β^)}-1{∑i=1nR˙i(σ^,α^,β^)R˙iT(σ^,α^,β^)}{-∑i=1nR¨i(σ^,α^,β^)}-1
is more robust than V^-1, where R˙i denotes the gradient vector of partial derivatives.To choose the covariates in GLMM defined by (6.1) and (6.2), Lai et al. [80, 81] propose to use a model selection procedure that consists of stepwise forward addition followed by stepwise backward elimination. A forward addition step aims at adding to the current model the most significant covariate among all possible candidates, while a backward elimination step aims at excluding from the current model the least significant candidate covariate. Specifically, for forward addition, they propose to choose one explanatory variable, amongst K possible candidates, to add to a fitted model with k-1 covariates, with coefficients β~1,…,β~k-1, σ~, and α~. For the jth candidate, consider the kth component sk(j) of the score vector ∑i=1n‍R˙i(σ~,α~,(β~1,…,β~k-1,0)T), called the “Rao statistic,” that includes the k-1 chosen covariates and the jth candidate. Let v~k(j) be the kth diagonal element of -∑i=1n‍R¨i(σ~,α~,(β~1,…,β~k-1,0)T). The candidate that maximizes |sk(j)|/v~k(j) is chosen for forward addition. For backward elimination, they propose to eliminate one covariate from a fitted model with k covariates, with coefficients β~1,…,β~k, σ~, and α~. They consider the “Wald statistic” |β^j|/v^j for testing βj=0, where v^j is the square root of the jth diagonal element of V^-1={-∑i=1n‍R¨i(σ^,α^,β^)}-1. If the elimination of the covariate with the smallest Wald statistic leads to a new model with a smaller information criterion (e.g., BIC), then the new model is preferred. The stepwise forward addition part of the model selection procedure terminates when the chosen information criterion for model selection does not improve with the addition of a covariate. It also terminates at the current model if there is no more candidate covariate to be included. At that point, backward elimination starts and proceeds until the selected information criterion does not improve the current model by deleting an explanatory variable. Since the estimation algorithm is likelihood based, Schwarz's Bayesian Information Criterion (BIC) can still be used for model selection:
(6.18)BIC=-2∑i=1nRi(σ^,α^,β^)+(logn)(numberofparameters).
If the sandwich estimator (6.17) is to be used instead of V^-1, then in the forward addition and backward elimination procedure above, v^j should be replaced by the square root of the jth diagonal entry of (6.17).


### 6.4. Prediction in GLMM

The main goal of default modeling in many other applications in econometrics is to predict certain characteristics of a specific subject given its covariate vector, for example, to forecastμi,t+1=g-1(βTxi,t+1+biTzi,t+1) at time t with observed covariate vector (xi,t+1T,zi,t+1T)T. More generally, we would like to estimate a function of the unobservable bi and we denote it by h(bi). Now, if σ,α, and β are known, then h(bi) should be estimated by the posterior mean Eσ,α,β[h(bi)|ithsubject'sdata]. Without assuming σ,α,β being known, the empirical Bayes approach replaces them in the posterior mean by their estimates (σ^,α^,β^) so that h is estimated by
(6.19)h^=Eσ^,α^,β^[h(bi)∣ithsubject'sdata].
The expectation in (6.19) can be evaluated by the hybrid method, which we have used for likelihood calculation; see [80, 81]. First note that Laplace's approximation in (6.15) is based on the Taylor expansion
(6.20)li(b)≈li(b^i)+(b-b^i)l¨i(b^i)(b-b^i)T2.
Since the density function of bi given the ith risk class data is proportional to eli(b), it follows from (6.20) that the conditional distribution of bi given the data is approximately normal with mean b^i and covariance matrix Vi-1, where Vi=-l¨i(b^i). Hence for a subject with informative data (i.e., satisfying λmin(V)≥c), for whom we drop the subscript i in li, Vi, and b^i, we can use the preceding normal density to evaluate h^ in (6.19). Consider
(6.21)h(b)≈h(b^+V-1/2Z)≈h(b^)+h˙(b^)TV-1/2Z+12ZTV-1/2h¨(b^)TV-1/2Z,
where Z is a q-dimensional standard normal distribution. Using the preceding approximation gives
(6.22)h^=h(b^)+12V-1/2h¨(b^)TV-1/2
for the informative subject.When the eigenvalue criterion in the hybrid method fails (i.e., whenλmin(V)<c), (6.19) can be evaluated by using importance sampling for the Monte Carlo approximation
(6.23)h^≈∑h=1H[h(b(h))wh]∑h=1Hwh,
where {b(1),…,b(H)} are independent samples from the trial density λ(·) and the importance weights are wh=el(b(h))/λ(b(h)). Similar to the hybrid scheme presented in Section 6.3, λ(·) is chosen to be a mixture of the prior normal distribution with density ϕα and the posterior normal distribution N(b^i,[-l¨i(b^i∣σ,α,β)+εI]-1), where ε is a small positive number to ensure that the covariance matrix is invertible.


### 6.5. Information Sets and Selection of Variables for Prediction

We have assumed so far that the observations(Yi,t,xi,t,zi,t) are available at every 1≤t≤T, for all 1≤i≤n. In longitudinal data in biostatistics, however, there is often between-subject variations in the observation times. Lai et al. [82] recently addressed this difficulty by using a prediction approach that customizes the predictive model for an individual by choosing predictors from the individual's information set and using the data from all subjects whose information sets contain these predictors to estimate the parameters of the predictive model. Let Ti be the set of times when the ith subject is observed, and let ui,t be the subject's covariate vector at time t∈Ti. The marginal regression approach in longitudinal data analysis models E(Yi,t∣ui,t), t∈Ti, by using generalized linear models, under the independence assumption over time or other working dependence assumptions (such as exchangeability, or autoregression when Ti is regularly spaced) as in the GEE (generalized estimating equation) methodology introduced by Liang and Zeger [83]. It has been shown by Pepe and Anderson [84] that in ignorance of the true dependence structure over time and in the presence of time-varying covariates, the GEE parameter estimates may be biased unless the working independence assumption is used. In contrast to the marginal regression approach, the transitional modeling approach incorporates outcome and covariate histories to model E(Yi,t∣Yi,s,ui,τ;s∈Tiands<t,τ∈Tiandτ≤t) via generalized linear models, by estimating the model parameters from all subjects who are observed at all times t∈Ti. A compromise between marginal regression and transitional modeling is the partly conditional mean approach proposed by Pepe and Couper [85]. Instead of working with the subject-specific set Ti, this approach considers the set T of all possible times and introduces an indicator variable Ri(t) and a stratifying variable Ki(t) such that Ri(t)=1 (or 0) if the ith subject is observed (or not observed) at time t and Ki(t) takes the values 1,…,K to represent the stratum to which subject i belongs. Lai et al. [82, Section 2.1] have pointed out that partly conditional mean modeling can be viewed as a special case of a more general approach that chooses predictors from a subject's information set. This more general approach allows one to use ideas similar to those used in information criteria such as BIC or AIC to choose predictors from a subject's information set. More importantly, it provides a framework for Lai et al. [3] to extend and refine the dynamic EB approach as follows.To begin with, instead ofK strata in the partly conditional mean approach, we divide the subjects into K structurally similar subgroups. In many applications, subjects belonging to the same subgroup have similar observation times because of their structural similarity. For example, patients who have more serious ailments are monitored more frequently than others in a study cohort, causing the irregularity of observation times over different subgroups. Lai et al. [3] assume the cross-sectional dynamics (6.12) separately for each subgroup, that is, with μt and Y-t-j replaced by μt(k) and Y-t-j(k) for the kth subgroup, in which Y-s(k) is the sample average from all subjects (from the subgroup) who are observed at time s. Concerning (6.13) for μi,t with i belonging to the kth subgroup, we choose only predictors from the information set of subject i that are common to all subjects in the kth subgroup. Lai et al. [3] propose to use the BIC, details of which are given in [80, page 58-59], for the GLMM (6.13) with i belonging to the kth subgroup and Y-t-j replaced by Y-t-j(k), and give examples from the problem of predicting the batting performance of baseball players studied previously by Efron and Morris [86, 87] and Brown [88] to illustrate this point.


### 6.6. Evolutionary Credibility via GLMM (6.13)

In Section6.1 we have used the dynamic EB approach to elucidate the close connection between Kalman filtering in the linear evolutionary credibility theory in [62] and the linear mixed model (6.11). Section 6.2 extends (6.11) to the GLMM (6.13) to handle generalized linear models of the form (6.1). The usefulness of (6.1) for insurance claims data has been pointed out by Hsiao et al. [89] who note that these claims data often exhibit “excess zeros,” meaning that a large fraction of policies during a given year do not have claims. They propose to handle excess zeros by using a two-part model of the form
(6.24)Yit=IitZit,i=1,…,n;t=1,…,T,
where Iit=I{Yit≠0} and Zit has the conditional distribution of Yit given that Yit≠0. In this two-part model, if Yit is the claim size of the ith policy at time t, then Zit is the claim size given that the ith policy files a claim at time t.


#### 6.6.1. Two-Part Nonlinear State-Space Models

The state-space model for linear evolutionary credibility in Section6.1 can be extended to the two-part model (6.24) as follows:
(6.25)Iit~Bernoulli(πit),Zit~Gamma(κ,μitκ)ifIit=1,
where κ is the shape parameter and link functions of πit=E(Iit) and μit=E(Zit) undergo AR(1) dynamics of the form
(6.26)logit(πit)=ρlogit(πi,t-1)+(1-ρ)logit(πi)+ηit,log(μit)=ρ′log(μi,t-1)+(1-ρ′)log(μi)+ϵit.
The function logit(π)=log(π/(1-π)) is the canonical link function for the exponential family of Bernoulli random variables. The ηit and ϵit in (6.26) are assumed to be independent normal random variables with mean 0, Var(ηit)=V, and Var(ϵit)=V'.Whereas Kalman filtering can be used to estimate the unobserved statesμt in the linear state-space model for evolutionary credibility in Section 6.1, the state-space model (6.25)-(6.26) is nonlinear, for which the Kalman filter is replaced by nonlinear filters to estimate the unobservable states πit and μit. These nonlinear filters are highly complex and require Markov chain Monte Carlo methods for their implementation when the model hyperparameters ρ,ρ',V,V',π1,…,πn,μ1,…,μn are known. Since they are unknown in practice, there is the additional complexity of estimating them, for which a stochastic EM approach implemented by Monte Carlo methods can be used to address the computational difficulties, as we have discussed in Section 4.1 for another application.


#### 6.6.2. Two-Part Dynamic Modeling of Claims via GLMM

Since the Bernoulli and gamma distributions are special cases of the exponential family, we can apply the GLMM introduced in Section6.2 to obtain a more tractable two-part model of evolutionary credibility, replacing (6.26) by
(6.27)logit(πit)=ρ1logit(I-t-1)+μ'+ai,log(μit)=ρ2log(Z-t-1)+μ+bi,
in which ai~N(0,v') and bi~N(0,v) are random effects, I-s=n-1∑i=1n‍Iis, and
(6.28)Z-s=∑i:Iis≠0Zs∑i:Iis≠01
is the average claim size over all positive claims at time s. The parameters of the GLMM, which uses the cross-sectional means I-t-1 and Z-t-1 to substitute for the unobservable states πi,t-1 and μi,t-1 in (6.26), are μ,v,μ', and v'. As pointed out in Section 6.2, the GLMM can also incorporate additional covariates xit and zit as in (6.13) to improve the predictive power of the model.In the state-space model (6.25)-(6.26), the conditional distribution of the unobservable state vector (πit,μit) given the observed data Dt is difficult to represent analytically or algorithmically, unlike that in the linear Gaussian state-space model for which the conditional distribution is normal and its mean and covariance matrix have a closed-form recursive formula, that is, the Kalman filter. Instead of specifying the stochastic dynamics of the state vector precisely, West et al. [90] assume a conjugate family of prior distributions of the exponential family in the generalized linear model (6.1) so that the parameters of the conjugate family can be updated similarly to the Kalman filter. Lai and Sun [73] have compared the prediction performance of this approximation to the nonlinear HMM filters for the two-part state-space model (6.25)-(6.26) with that of the GLMM (6.27) in a simulation study, in which the data are generated by the two-part state-space model, and have found that the two methods perform similarly. An advantage of using GLMM to model evolutionary credibility over nonlinear state-space models or their approximations proposed in [90] is that GLMM has a relatively complete and easy-to-use methodology for likelihood inference and model selection; see [3, 80]. In particular, we can use model selection criteria such as BIC to choose the order p of an AR(p) model for μit, or to determine which of the variables in xit and zit should be included in the model to avoid overfitting.State-space models have traditionally been considered in the context of a singled-dimensional series. A distinctive feature of evolutionary credibility is that it involves many structurally similar time series, called “dynamic panel data” in econometrics, and we can exploit this feature to derive a dynamic empirical Bayes approach to evolutionary credibility. As pointed out in [91, page 335], traditional credibility theory can be called “empirical Bayes credibility” since it basically amounts to replacing the hyperparameters in the Bayes credibility (assuming Gaussian claims and a Gaussian prior) by their method-of-moment estimates. Extending this empirical Bayes approach to the Bayesian model of linear evolutionary credibility, which is a linear state-space model, yields the LMM in Section 6.1. Similarly, the GLMM in this section can be interpreted as a dynamic empirical Bayes approach to evolutionary credibility in generalized linear models.


### 6.7. Competing Risks and Retail Loans

We now return to the problem of predicting the default probabilities and loss given default of mortgages and corporate bonds that we have mentioned in Section1. The prediction problem is important for monitoring and management of a bank's credit risk. In 2001, the Basel Committee on Banking Supervision [92] proposed a new accord for risk evaluation, which was finalized in 2005 and is known as the Basel II Accord. Using the internal ratings-based (IRB) approach, the expected credit loss E(Lc) for a loan is decomposed into
(6.29)E(Lc)=EAD×LGD×PD,
where PD is the probability of default, LGD is the expected loss given default, which is expressed as a rate between 0 and 1, and EAD is the exposure at default. We first consider default probabilities of retail loans (such as personal loans, mortgages, credit cards, and car loans). “Retail” in banking means consumer-related activities, while “wholesale” means business-related loan transactions which include commercial loans, commercial real estate loans, and leases to businesses.Retail loans are typically grouped into several classes calledrisk buckets, within which the obligors can be regarded as having the same prior probability of default. Risk bucketing is tantamount to stratifying the obligors into classes so that the baseline covariate (such as credit score at the time of loan application) of an individual obligor can be replaced by the class membership. For retail loans, default means overdue payment for longer than 90 days. Thus, unlike a corporation that files for bankruptcy and the LGD of the loan is determined later by the recovery rate, a retail loan borrower can default and incur a loss to the lender while still holding the loan until it is foreclosed and taken off the balance sheet. Moreover, besides the possibility of default, another possible event that the lender faces is prepayment.


#### 6.7.1. Competing Risks

Calhoun and Deng [93] regard default and prepayment as competing risks and introduce a multilogit regression model of the trinomial variable Y(t), taking the values 2, 1, and 0 that correspond to prepayment, foreclosure due to default, and still surviving at time t, respectively. In our setting that considers the trinomial variable Yij(t) of the ith obligor in the jth risk class, the multilogit regression model relates Yij(t) to a covariate vector uij(t) by
(6.30)P(Yij(t)=k∣Yij(s)=0fors<t,uij(t))=exp(αjk+βjkTuij(t))1+∑h=12exp(αjh+βjhTuij(t))
for k=1,2. Competing-risks models can be viewed as “marked” failure-time models, in which a mark is assigned to a failure event, indicating the type of the event that has occurred.Competing-risks data are also prevalent in cancer studies for which the response to treatment can be classified in terms of failure from disease processes and non-disease-related causes. For example, after treatments for leukemia and lymphoma, patients may either relapse or die in remission. In quality-of-life studies, the risks of the different failure types need to be considered for cost-effectivness analyses. Although Cox's proportional hazards model has routinely been used to analyze competing-risks data with covariates, what it models is thecause-specific hazard rate(6.31)λj(t;Z)=limh→0h-1P{t≤Y≤t+h,J=j∣Y≥t,Z},
where Z is a vector of covariates, Y is the failure time, and J is the cause of failure, with the causes labeled 1,2,…,j*. This approach treats each cause separately as the only cause for the observed failures that can be attributed to it and treats the observed failures due to other causes as censoring outcomes. Viewing competing risks as multiple endpoints each of which precludes the occurrence of the others, Fine and Gray [94] proposed to consider, instead of (6.31), the hazard function of the cumulative incidence due to cause j:
(6.32)λ~j(t;Z)=limh→0h-1P[t≤Y≤t+h,J=j∣{Y≥t}∪{Y<t,J≠j},Z]=-ddtlog(1-Fj(t;Z)),
where Fj(t;Z) is the cumulative incidence function(6.33)Fj(t;Z)=P{Y≤t,J=j∣Z}
introduced by Gray [95]. Note that 1-Fj(t-;Z) = 1-P{Y<t,J=j∣Z} = P{Y≥t∣Z} + P{Y<t,J≠j∣Z}. In the absence of covariates, the cumulative incidence function can be estimated by the Aalen-Johansen estimator [27, 28] based on a finite-state Markov chain model of competing risks; see the survey by Klein and Shu [96]. Although the proportional hazards model λj(t;Z)=λj0(t)eβTZ can readily be extended to the hazard functions λ~j of the competing risks, regression analysis of λ~j under the model λ~j(t;Z)=λ~j0(t)exp(βTZ) becomes much more difficult. Three approaches to this problem have been developed. First, Fine and Gray [94] used the following modification of risk sets to extend partial likelihood to competing-risks data. The risk set at time t includes the subjects that are naturally at risk and those who have failed prior to time t due to other competing risks and are not censored at time t. A subject who has failed prior to time t due to other competing risks is weighted by the probability that it is not censored at time t given that it has not been censored at the failure time. The second method, proposed by Scheike et al. [97], uses weighted counting processes as responses in a binomial regression model for cumulative incidence functions. Both methods assume a common censoring distribution. Their performance is highly sensitive to how well the censoring distribution can be estimated. The third method, proposed by Klein and Andersen [98], uses jackknife pseudovalues to model the cumulative incidence functions directly. These pseudo-values are used in a generalized estimating equation to obtain estimates on model parameters, and Graw et al. [99] have recently shown this method to be consistent and asymptotically normal.An important difference between the retail loans dataYij(t) in (6.30) and those in cancer biostatistics is that the latter is subject to right-censoring due to patient withdrawal whereas the former remain in the lender's balance sheet until prepayment or default leading to foreclosure. Right-censoring results in much greater complexity than multilogit regression in modeling competing risks.


#### 6.7.2. GLMM for Loan Default Probabilities

Calhoun and Deng [93] fitted (6.30) to about 1.3 million conforming fixed-rate and adjustable-rate residential mortgages originated between 1979 and 1993. Their covariate vector includes loan age, dummy variables for origination years, the loan-to-value ratio (LTV), seasonality, occupancy status of the property, the original loan size, the ratio of the 10-year Constant Maturity Treasury rate to the 1-year Constant Maturity Treasury rate, mortgage premium value, and the probability of negative equity based on estimated drift and volatility of housing prices. Instead of the multilogit model (6.30) that generalizes logistic regression, Lai et al. [100] use the GLMM in Section 6.2 to model the trinomial variable Yij(t) that represents continuing payment, or default, or prepayment at time t of the ith obligor in the jth risk class. Whereas [93] does not divide the borrowers into risk classes, the statistical analysis of subprime mortgages in [100] conforms to risk bucketing that is commonly used for risk management of retail loans and is based on loan-level mortgage data consisting of 10,000 first-lien, 2–28 adjustable-rate subprime loans originated between January 2004 and June 2006; these loans are randomly selected from a database maintained by LoanPerformance, Inc., a mortgage information and analytics provider. The 10,000 loans are grouped into I=5 risk classes according to the FICO score when a loan was originated; the class with the lowest FICO scores is labeled 1 and the highest labeled 5. These data involve two time scales, calendar time t, and loan age s, which also arise in biostatistics in the context of time-sequential clinical trials with failure-time endpoints and staggered entry, where s=(t-τ)+ and τ is the entry time of the patient into the study (or in the case of mortgages the initiation date of the loan); see [101]. The case s=0 denotes that the loan has not been originated at calendar time t.Lai et al. [100] use the link function
(6.34)ηijt(k)=log(P{Yij(t)=k∣Yij(t-1)=0,uij(t)}P{Yij(t)=0∣Yij(t-1)=0,uij(t)}),k=1,2,
to express (6.30) in the generalized linear form, which can then be extended to the GLMM in Section 6.2:
(6.35)ηijt(k)=ρ(k)log(Y-t-1(k)Y-t-1(0))+aj(k)+β(k)TXijt+bj(k)TZijt,
where Xijt and Zijt are covariates that are available prior to t and Y-t-1(k) is the cross-sectional mean of the indicator variables I{Yij(t-1)=h,Yij(t-2)=0} for k=0,1,2. If the loan has not been initiated at time t-1, then the conditional probabilities in (6.34) are taken to be P{Yij(t)=·∣uij(t)} and the cross-sectional term involving Y-t-1(k)/Y-t-1(0) is also removed from (6.35). The covariate vector Xijt includes the loan age, occupancy status, initial mortgage rate, purpose (for purchase or refinance), a dummy variable for full documentation, or otherwise (important for reflecting the lender's standards in issuing the subprime loan), and Zt includes mortgage premium value and equity position as defined in [100], which involve the macroeconomic variables 1-year Constant Maturity Treasury rate and zip-code-based house price index, together with the loan's size, unpaid rate, mortgage rate, and LTV. These data do not contain the loan's payment history; ideally Yij(t) should also indicate if the payment is overdue for longer than 90 days even when the loan is not foreclosed, but such information is not available in the data obtained from LoanPerformance, only showing the trinomial variable if the loan is prepaid, foreclosed, or active at time t. The advantages of (6.35) over the multilogit model (6.30) in predicting defaults are discussed in [102].


#### 6.7.3. Loss Given Default

Lai and Wong [103, Section 5.2] have described how GLMM can be used to model the LGD of retail loans. They first review how LGD is calculated for corporate loans and bonds in the software package LossCalc of Moody's KMV, which is built from a training sample of recovery data of the defaulted corporate loans, bonds, and preferred stock. The training sample is used to fit a linear regression model of y=Φ-1(Fα,β(RR)) on nine predictors chosen from a set of regressions using a year-by-year out-of-sample prediction error criterion, where RR is the recovery rate implied by the market value (bid-side market quote) one month after default, Φ is the standard normal distribution function, and Fα,β is the distribution function of the Beta(α,β) distribution. The transformation Φ-1∘Fα,β is used to convert a beta random variable into a normal one so that linear regression can be applied. As shown in [103, page 237], a better approach is to fit the beta regression model, which is a generalized linear model.As pointed out in [103], unlike corporate loans are still quite liquid after default, personal loans, auto loans, and mortgage loans lack such liquidity and the usual recovery process involves (a) in-house departments, (b) selling off the debt at some fraction of the face value, and (c) collection or foreclosure agencies which take a fixed percentage of the recovered loan as commission. To describe simply the recovery process, the bank first tries (a) and switches to (b) or (c) if (a) is unsuccessful after, say, three months. In the case of mortgage loans, (b) means selling the mortgage to another mortgage company (e.g., Fannie Mae or Freddie Mac in the US), while (c) means foreclosure of the mortgage and repossession and sale of the property. The shortfall is the residual loan value plus planned interest minus the sales price of the mortgage or the repossessed property. For occasional missed payments, the “defaulted” obligors who have forgotten to pay or have been away usually pay within one month after they are contacted by the in-house department. Such a scenario corresponds to complete recovery and leads to an atom at 0 in the marginal distribution of the economic LGD. Lai and Wong [103] point out the importance of including the atom 0 and demonstrate the advantages of using beta regression over linear regression of the transformed recovery rate on covariates. They also note that the Basel II Accord requires LGD estimates to be grounded in historical recoveries and to incorporate business cycles since at times of economic downturn both PD and LGD tend to increase; the database should ideally cover at least one business cycle and must not be shorter than five years. They propose the following GLMM to estimate LGD using these data.LetDs denote the set of obligors who default in year s or before and who have available credit and collateral information in year s-1, with s=t,t-1,…,t-4. For i∈Ds, let Yi(s) be the economic LGD in year s. As noted above, Yi(s) has an atom at 0. Let πi(s)=P{Yi(s)>0}, which can be modeled by
(6.36)logit(πi(s))=ξi+ηTxi(s-1)+γTf(s),
where xi(s-1) is a vector of attributes of the obligor's creditworthiness and collateral in year s-1 and f(s) is a vector of macroeconomic factors in year s. The conditional density of Yi(s)∣Yi(s)>0 is assumed to belong to the exponential family of beta density functions
(6.37)f(y∣μi(s),ϕ)=Γ(ϕ)Γ(μi(s)ϕ)Γ((1-μi(s))ϕ)yμi(s)ϕ-1(1-y)(1-μi(s))ϕ-1
for 0<y<1, where μi(s)=E(Yi(s)∣Yi(s)>0),
(6.38)logit(μi(s))=ζi+κTxi(s-1)+νTf(s),Var(Yi(s)∣Yi(s)>0)=μi(s)(1-μi(s))/(1+ϕ), and ξi are i.i.d. normal and are independent of ζi which are also i.i.d. normal. This GLMM includes loan-specific random effects ξi and ζi in (6.36) and (6.38) and is similar to the two-part model in Section 6.6. It enables one to predict directly LGDi(t+1) of the ith obligor for year t+1 using xi(t) and the forecast of f(t+1), noting that LGDi(t+1)=πi(t+1)E(Yi(t+1)∣Yi(t+1)>0).


### 6.8. GLMM for Modeling Frailty and Contagion of Credit Portfolios

Similarly to the applications to evolutionary credibility and mortgage loans in the preceding two sections, the dynamic EB approach via GLMM can be applied to provide a new model for correlated default intensities of multiple firms in a portfolio of corporate bonds or a multiname credit derivative. This model can incorporate both frailty and contagion via the covariates of the GLMM and is considerably easier to implement. The details are given in [3, 100]. Note that the use of a mixture distribution for firm-specific random effects also appears in Vasicek's mixed binomial model (3.11) for the number of defaults in a loan portfolio up to a horizon T. It is shown in [3] that the prediction performance of the nonlinear state-space model of Duffie et al. [46] is similar to that of the much simpler GLMM model of dynamic frailty even when the data are actually generated by the nonlinear state-space model (4.2).


## 6.1. Linear Evolutionary Credibility, Kalman Filtering, and Dynamic EB

Noting that insurers observe claims of risk classes over successive periods, Frees et al. [64, 65] in fact consider the setting of longitudinal data in their LMM approach to credibility theory; the index j in (5.11) is replaced by t that denotes time so that the Yit represent longitudinal data of claims. Bühlmann and Gisler [62] go further in assuming a dynamic Bayesian model in which the mean θit of Yit has a prior distribution with mean μt that changes with time according to the AR(1) model
(6.3)μt=ρμt-1+(1-ρ)μ+ηt,
with i.i.d. zero-mean ηt having variance V. Thus, the μt are latent states undergoing AR(1) dynamics, yielding a linear state-space model for which the unobserved states μt can be estimated from the observations Yis,s≤t, by the Kalman filter μ^t|t given by the recursion
(6.4)μ^t∣t=μ^t∣t-1+ρ-1Kt(Yt-μ^t∣t-11),μ^t+1∣t=ρμ^t∣t+(1-ρ)μ,
in which Yt=(Y1t,…,Ynt)T, 1=(1,…,1)T, and Kt is the 1×n Kalman gain matrix defined recursively in terms of the hyperparameters V=Var(ηt), vt=Var(Yit∣μt), and ρ by
(6.5)Kt=ρΣt∣t-11T(Σt∣t-111T+vtI)-1,(6.6)Σt+1∣t=(ρ-Kt1)2Σt∣t-1+V+vtKtKtT.
Under the assumption that θit∣μt has variance τ2 and Yit∣θit has variance a+bθit+cθit2 as in (5.3), vt=τ2+a+bθit+cθit2. The Kalman filter is the minimum-variance linear estimator of μt; it is also the Bayes estimator if Yit (conditional on μt) and ηt are normal.As in the linear EB approach described in Section5, we can use the method of moments to estimate μt and vt for every 1≤t≤T. Specifically, Y-t=n-1∑i=1n‍Yit is a consistent estimate of μt and st2=∑i=1n‍(Yit-Y-t)2/(n-1) is a consistent estimate of vt. Moreover, assuming E(μ0)=μ, it follows from (6.3) that E(μt)=μ for all t. Therefore, μ can be consistently estimated from the observations up to time t by μ^(t)=(∑s=1t‍Y-s)/t, and ρ and V can be consistently estimated by
(6.7)ρ^(t)=∑s=2t(Y-s-μ^(t))(Y-s-1-μ^(t))∑s=2t(Y-s-μ^(t))2,V^(t)=∑s=2t{Y-s-ρY-s-1-(1-ρ)μ^(t)}2(t-2),
for t≥2. Replacing the hyperparameters μ, ρ, V, and vt in the Bayes estimate μ^t|t of μt by their estimates μ^(t), ρ^(t), V^(t), and st2 yields the EB estimate of μt. Of particular interest in applications of longitudinal data is the prediction of Yi,t+1 based on the set of observations 𝒴t up to time t. The Bayes predictor, assuming normal priors and normal observations, is
(6.8)E(Yi,t+1∣𝒴t)=E(θi,t+1∣𝒴t)=E(μt+1∣𝒴t)=μ^t+1∣t,
which is the Kalman predictor given by (6.4). This is also the best linear predictor without assuming normality. The corresponding EB predictor replaces the hyperparameters in the Kalman predictor by their method-of-moment estimates μ^(t), ρ^(t), V^(t), and st2.Note that to estimate the hyperparameters in this EB approach via (6.3) and Kalman filtering, we have used the cross-sectional mean Y-t-1 of n independent observations that have mean μt-1. Lai and Sun [73] propose an alternative and more direct approach that replaces μt-1 in (6.3) by Y-t-1, leading to
(6.9)μt=ρY-t-1+ω+ηt
with ω=(1-ρ)μ, which can be written as an LMM for the longitudinal data Yij:
(6.10)Yit=ρY-t-1+ω+bi+ϵit,
by absorbing ηt into ϵit. Since (6.10) is in the form of a regression model, Lai and Sun [73] also include additional covariates to increase the predictive power of the model in the LMM
(6.11)Yit=ρY-t-1+ai+βTxit+biTzit+ϵit,
where ai and bi are subject-specific random effects, xit represents a vector of subject-specific covariates that are available prior to time t (for predicting Yit prior to observing it at time t), and zit denotes a vector of additional covariates that are associated with the random effects bi. A simulation study in [71] to compare the performance of parameter estimation and one-step-ahead prediction using (6.11) with that of Kalman filtering shows that the dynamic LMM (6.11) yields results comparable to those of Kalman filtering when the data are generated by a linear state-space model of the type (6.5), the parameters of which are not assumed to be known and are estimated by (6.7).


## 6.2. Dynamic EB Modeling via GLMM

Lai and Sun [73] assume the prior distribution for which the μit are i.i.d. with mean μt and
(6.12)g(μt)=∑j=1pθjg(Y-t-j)for1≤t≤T,whereY-s=n-1∑i=1nYis.
The dynamic model (6.12) for μt is an EB version of the Markov model introduced by Zeger and Qaqish [74], who consider a single time series Yt that has a density function of the form (6.1) and models μt by g(μt)=∑j=1p‍θjg(Yt-j), where g is a link function and n=1 in this case. For general n, note that μs is the mean of μis and can be consistently estimated by Y-s, which is the basic idea underlying the linear EB approach. This suggests using g(Y-t-j) in lieu of g(Yt-j) for n≥1 that leads to (6.12). As in the LMM (6.11), Lai and Sun [73] propose to increase the predictive power of the model by including fixed and random effects and other time-varying covariates of each subject i, thereby removing the dependence of g(μit)-g(μt) on t in the GLMM
(6.13)g(μit)=∑j=1pθjg(Y-t-j)+ai+βTxit+biTzit,
in which θ1,…,θp and β are the fixed effects and ai and bi are subject-specific random effects. Following [4], it is assumed that ai and bi are independent normal with zero means. Denote the unknown vector specifying the covariances of the random effects by α, and for notational simplicity, we can augment bi to include ai so that (6.13) can be written as g(μit)=∑j=1p‍θjg(Y-t-j)+xitTβ+(1,zitT)bi. Note that (6.13) is a natural extension of the dynamic linear EB model (6.11) to the exponential family. By including g(Y-t-j),j=1,…,p, in xit, we can implement the dynamic EB model (6.13) by using the methods described for general GLMMs.


## 6.3. Implementation and Model Selection

Breslow and Clayton assume thebi in (6.2) to have a common normal distribution with mean 0 and covariance matrix Σ that depends on an unknown parameter. Lai and Shih [75] have developed a nonparametric method to estimate the random effects distribution in the more general setting of nonlinear mixed effects models and have shown that there is very low resolution in estimating the distribution. This provides a heuristic explanation for why the choice of a normal distribution, with unspecified parameters, for the random effects distribution G does not result in a worse estimate of β and σ than the semiparametric approach that estimates G nonparametrically. An advantage of assuming a normal distribution for the random effects bi is that it only involves the covariance matrix of the random effects. The likelihood function of the GLMM defined by (6.1) and (6.2) is of the form ∏i=1n‍Li(σ,α,β), where
(6.14)Li(σ,α,β)=∫{∏t=1Tf(yit;θit,σ)}ϕα(b)db,
in which ϕα denotes the normal density function with mean 0 and covariance matrix depending on an unknown parameter α. We review three methods to compute the likelihood function, the maximizer of which gives the MLE of σ, α, and β.


### 6.3.1. Laplace's Approximation

Lettingeli(b∣σ,α,β) be the integrand in the right-hand side of (6.14), Laplaces's asymptotic formula for integrals yields the approximation
(6.15)∫eli(b∣σ,α,β)db≈(2π)q/2{det[-l¨i(b^i∣σ,α,β)]}-1/2exp{li(b^i∣σ,α,β)},
where q is the dimension of bi, b^i=b^i(σ,α,β) is the maximizer of li(bσ,α,β), and l¨i denotes the Hessian matrix consisting of second partial derivatives of li with respect to the components of b. The R package lme4 computes the MLE by using the Laplace approximation to (6.14) or the restricted pseudolikelihood approach proposed by [76], as the user-specified option.


### 6.3.2. Gauss-Hermite Quadrature

Laplaces's asymptotic formula (6.15) is derived from the asymptotic approximation of li by a quadratic function of b in a small neighborhood of b^i as λmin(-l¨i(b^i∣σ,α,β))→∞, where λmin(·) denotes the minimum eigenvalue of a symmetric matrix. Therefore, such formula may produce significant approximation error for Li if the corresponding λmin(-l¨i(b^i∣σ,α,β)) is not sufficiently large. One way to reduce the possible approximation error is to compute Li by using an adaptive Gauss-Hermite quadrature rule, as in [77]. The software package SAS uses adaptive Gauss-Hermite quadrature in the NLMIXED procedure to compute (6.14); the R function lmer(·) also uses Gaussian quadrature to compute (6.14) but only for certain special cases of the exponential family (6.1).


### 6.3.3. Hybrid Method

The numerical integration procedure demands a much higher computational effort and becomes computationally infeasible ifn or q is large. To circumvent the issue of high-dimensional numerical integration, Antonio and Beirlant [72] propose to put prior distributions on the unknown parameters and to estimate them by using the Markov chain Monte Carlo (MCMC) method in a Bayesian way. The performance of the MCMC method, however, depends on how the prior parameters are set as well as the convergence rate of the Markov chain to its stationary distribution, which may not even exist. Yafune et al. [78] use direct Monte Carlo integration but point out that the computational time may be too long to be of practical interest. Instead of relying solely on Monte Carlo, Lai and Shih [79] proposed the following hybrid approach. Let Vi=-l¨i(b^i|σ,α,β). Since Laplaces's asymptotic formula (6.15) may be a poor approximation to (6.14) when λmin(Vi) is not sufficiently large, Monte Carlo integration, whose error is independent of q, can be used as an alternative method to evaluate (6.14).Lai and Shih [79] suggest using Monte Carlo integration instead of Laplace's asymptotic formula for i whose λmin(Vi)<c, where c is a positive threshold. Lai et al. [80, 81] further improve the hybrid method by using importance sampling in the Monte Carlo component of the hybrid method. Specifically, instead of sampling b(h) directly from ϕα as in [79], sample it from a mixture of the prior normal distribution with density ϕα and the posterior normal distribution N(b^i,[-l¨i(b^i∣σ,α,β)+εI]-1), where ε is a small positive number to ensure that the covariance matrix is invertible. We denote such choice of trial density as λ. It has the advantage of further incorporating the essence of Laplace's method in the Monte Carlo step such that the method is less dependent (than direct Monte Carlo) on the choice of the threshold c.(i)
Ifλmin(Vi)<c, evaluate Li(σ,α,β) in (6.14) by the Monte Carlo approximation
(6.16){∑h=1H[wh∏t=1Tf(yit;θit(b(h),β),σ)]}{∑h=1Hwh}
in which b(1),…,b(H) are independently sampled from the importance density λ and wh=ϕα(b(h))/λ(b(h)) are the importance weights.(ii)
Ifλmin(Vi)≥c, use Laplace approximation (6.15). Lai et al. [80] propose to use c=10,H=100, and a mixture distribution that assigns 0.2 probability to ϕα. Note that choosing λ=ϕα yields wh=1, for which (6.16) reduces to direct Monte Carlo evaluation of Li(σ,α,β) in [79]. The likelihoods Li(σ,α,β) thus computed then yield the likelihood function ∏i=1n‍Li(σ,α,β), whose logarithm can be maximized by using iterative gradient-type schemes.


### 6.3.4. Model Selection

As shown by Breslow and Clayton [4], a consistent estimator of the asymptotic covariance matrix of (σ^,α^,β^) is V^-1={-∑i=1n‍R¨i(σ^,α^,β^)}-1, where Ri=logLi and R¨i denotes the Hessian matrix of second derivatives of Ri with respect to σ, α, and the components of β. As in [79], we compute the observed information matrix V^=-∑i=1n‍R¨i(σ^,α^,β^) by taking numerical derivatives and using the hybrid method to evaluate Li. Lai and Shih [79] have shown that V^-1/2(σ^-σ,α-α^,β-β^) has a limiting standard multivariate normal distribution as n→∞ under certain regularity conditions, and that the sandwich estimator
(6.17){-∑i=1nR¨i(σ^,α^,β^)}-1{∑i=1nR˙i(σ^,α^,β^)R˙iT(σ^,α^,β^)}{-∑i=1nR¨i(σ^,α^,β^)}-1
is more robust than V^-1, where R˙i denotes the gradient vector of partial derivatives.To choose the covariates in GLMM defined by (6.1) and (6.2), Lai et al. [80, 81] propose to use a model selection procedure that consists of stepwise forward addition followed by stepwise backward elimination. A forward addition step aims at adding to the current model the most significant covariate among all possible candidates, while a backward elimination step aims at excluding from the current model the least significant candidate covariate. Specifically, for forward addition, they propose to choose one explanatory variable, amongst K possible candidates, to add to a fitted model with k-1 covariates, with coefficients β~1,…,β~k-1, σ~, and α~. For the jth candidate, consider the kth component sk(j) of the score vector ∑i=1n‍R˙i(σ~,α~,(β~1,…,β~k-1,0)T), called the “Rao statistic,” that includes the k-1 chosen covariates and the jth candidate. Let v~k(j) be the kth diagonal element of -∑i=1n‍R¨i(σ~,α~,(β~1,…,β~k-1,0)T). The candidate that maximizes |sk(j)|/v~k(j) is chosen for forward addition. For backward elimination, they propose to eliminate one covariate from a fitted model with k covariates, with coefficients β~1,…,β~k, σ~, and α~. They consider the “Wald statistic” |β^j|/v^j for testing βj=0, where v^j is the square root of the jth diagonal element of V^-1={-∑i=1n‍R¨i(σ^,α^,β^)}-1. If the elimination of the covariate with the smallest Wald statistic leads to a new model with a smaller information criterion (e.g., BIC), then the new model is preferred. The stepwise forward addition part of the model selection procedure terminates when the chosen information criterion for model selection does not improve with the addition of a covariate. It also terminates at the current model if there is no more candidate covariate to be included. At that point, backward elimination starts and proceeds until the selected information criterion does not improve the current model by deleting an explanatory variable. Since the estimation algorithm is likelihood based, Schwarz's Bayesian Information Criterion (BIC) can still be used for model selection:
(6.18)BIC=-2∑i=1nRi(σ^,α^,β^)+(logn)(numberofparameters).
If the sandwich estimator (6.17) is to be used instead of V^-1, then in the forward addition and backward elimination procedure above, v^j should be replaced by the square root of the jth diagonal entry of (6.17).


## 6.3.1. Laplace's Approximation

Lettingeli(b∣σ,α,β) be the integrand in the right-hand side of (6.14), Laplaces's asymptotic formula for integrals yields the approximation
(6.15)∫eli(b∣σ,α,β)db≈(2π)q/2{det[-l¨i(b^i∣σ,α,β)]}-1/2exp{li(b^i∣σ,α,β)},
where q is the dimension of bi, b^i=b^i(σ,α,β) is the maximizer of li(bσ,α,β), and l¨i denotes the Hessian matrix consisting of second partial derivatives of li with respect to the components of b. The R package lme4 computes the MLE by using the Laplace approximation to (6.14) or the restricted pseudolikelihood approach proposed by [76], as the user-specified option.


## 6.3.2. Gauss-Hermite Quadrature

Laplaces's asymptotic formula (6.15) is derived from the asymptotic approximation of li by a quadratic function of b in a small neighborhood of b^i as λmin(-l¨i(b^i∣σ,α,β))→∞, where λmin(·) denotes the minimum eigenvalue of a symmetric matrix. Therefore, such formula may produce significant approximation error for Li if the corresponding λmin(-l¨i(b^i∣σ,α,β)) is not sufficiently large. One way to reduce the possible approximation error is to compute Li by using an adaptive Gauss-Hermite quadrature rule, as in [77]. The software package SAS uses adaptive Gauss-Hermite quadrature in the NLMIXED procedure to compute (6.14); the R function lmer(·) also uses Gaussian quadrature to compute (6.14) but only for certain special cases of the exponential family (6.1).


## 6.3.3. Hybrid Method

The numerical integration procedure demands a much higher computational effort and becomes computationally infeasible ifn or q is large. To circumvent the issue of high-dimensional numerical integration, Antonio and Beirlant [72] propose to put prior distributions on the unknown parameters and to estimate them by using the Markov chain Monte Carlo (MCMC) method in a Bayesian way. The performance of the MCMC method, however, depends on how the prior parameters are set as well as the convergence rate of the Markov chain to its stationary distribution, which may not even exist. Yafune et al. [78] use direct Monte Carlo integration but point out that the computational time may be too long to be of practical interest. Instead of relying solely on Monte Carlo, Lai and Shih [79] proposed the following hybrid approach. Let Vi=-l¨i(b^i|σ,α,β). Since Laplaces's asymptotic formula (6.15) may be a poor approximation to (6.14) when λmin(Vi) is not sufficiently large, Monte Carlo integration, whose error is independent of q, can be used as an alternative method to evaluate (6.14).Lai and Shih [79] suggest using Monte Carlo integration instead of Laplace's asymptotic formula for i whose λmin(Vi)<c, where c is a positive threshold. Lai et al. [80, 81] further improve the hybrid method by using importance sampling in the Monte Carlo component of the hybrid method. Specifically, instead of sampling b(h) directly from ϕα as in [79], sample it from a mixture of the prior normal distribution with density ϕα and the posterior normal distribution N(b^i,[-l¨i(b^i∣σ,α,β)+εI]-1), where ε is a small positive number to ensure that the covariance matrix is invertible. We denote such choice of trial density as λ. It has the advantage of further incorporating the essence of Laplace's method in the Monte Carlo step such that the method is less dependent (than direct Monte Carlo) on the choice of the threshold c.(i)
Ifλmin(Vi)<c, evaluate Li(σ,α,β) in (6.14) by the Monte Carlo approximation
(6.16){∑h=1H[wh∏t=1Tf(yit;θit(b(h),β),σ)]}{∑h=1Hwh}
in which b(1),…,b(H) are independently sampled from the importance density λ and wh=ϕα(b(h))/λ(b(h)) are the importance weights.(ii)
Ifλmin(Vi)≥c, use Laplace approximation (6.15). Lai et al. [80] propose to use c=10,H=100, and a mixture distribution that assigns 0.2 probability to ϕα. Note that choosing λ=ϕα yields wh=1, for which (6.16) reduces to direct Monte Carlo evaluation of Li(σ,α,β) in [79]. The likelihoods Li(σ,α,β) thus computed then yield the likelihood function ∏i=1n‍Li(σ,α,β), whose logarithm can be maximized by using iterative gradient-type schemes.


## 6.3.4. Model Selection

As shown by Breslow and Clayton [4], a consistent estimator of the asymptotic covariance matrix of (σ^,α^,β^) is V^-1={-∑i=1n‍R¨i(σ^,α^,β^)}-1, where Ri=logLi and R¨i denotes the Hessian matrix of second derivatives of Ri with respect to σ, α, and the components of β. As in [79], we compute the observed information matrix V^=-∑i=1n‍R¨i(σ^,α^,β^) by taking numerical derivatives and using the hybrid method to evaluate Li. Lai and Shih [79] have shown that V^-1/2(σ^-σ,α-α^,β-β^) has a limiting standard multivariate normal distribution as n→∞ under certain regularity conditions, and that the sandwich estimator
(6.17){-∑i=1nR¨i(σ^,α^,β^)}-1{∑i=1nR˙i(σ^,α^,β^)R˙iT(σ^,α^,β^)}{-∑i=1nR¨i(σ^,α^,β^)}-1
is more robust than V^-1, where R˙i denotes the gradient vector of partial derivatives.To choose the covariates in GLMM defined by (6.1) and (6.2), Lai et al. [80, 81] propose to use a model selection procedure that consists of stepwise forward addition followed by stepwise backward elimination. A forward addition step aims at adding to the current model the most significant covariate among all possible candidates, while a backward elimination step aims at excluding from the current model the least significant candidate covariate. Specifically, for forward addition, they propose to choose one explanatory variable, amongst K possible candidates, to add to a fitted model with k-1 covariates, with coefficients β~1,…,β~k-1, σ~, and α~. For the jth candidate, consider the kth component sk(j) of the score vector ∑i=1n‍R˙i(σ~,α~,(β~1,…,β~k-1,0)T), called the “Rao statistic,” that includes the k-1 chosen covariates and the jth candidate. Let v~k(j) be the kth diagonal element of -∑i=1n‍R¨i(σ~,α~,(β~1,…,β~k-1,0)T). The candidate that maximizes |sk(j)|/v~k(j) is chosen for forward addition. For backward elimination, they propose to eliminate one covariate from a fitted model with k covariates, with coefficients β~1,…,β~k, σ~, and α~. They consider the “Wald statistic” |β^j|/v^j for testing βj=0, where v^j is the square root of the jth diagonal element of V^-1={-∑i=1n‍R¨i(σ^,α^,β^)}-1. If the elimination of the covariate with the smallest Wald statistic leads to a new model with a smaller information criterion (e.g., BIC), then the new model is preferred. The stepwise forward addition part of the model selection procedure terminates when the chosen information criterion for model selection does not improve with the addition of a covariate. It also terminates at the current model if there is no more candidate covariate to be included. At that point, backward elimination starts and proceeds until the selected information criterion does not improve the current model by deleting an explanatory variable. Since the estimation algorithm is likelihood based, Schwarz's Bayesian Information Criterion (BIC) can still be used for model selection:
(6.18)BIC=-2∑i=1nRi(σ^,α^,β^)+(logn)(numberofparameters).
If the sandwich estimator (6.17) is to be used instead of V^-1, then in the forward addition and backward elimination procedure above, v^j should be replaced by the square root of the jth diagonal entry of (6.17).


## 6.4. Prediction in GLMM

The main goal of default modeling in many other applications in econometrics is to predict certain characteristics of a specific subject given its covariate vector, for example, to forecastμi,t+1=g-1(βTxi,t+1+biTzi,t+1) at time t with observed covariate vector (xi,t+1T,zi,t+1T)T. More generally, we would like to estimate a function of the unobservable bi and we denote it by h(bi). Now, if σ,α, and β are known, then h(bi) should be estimated by the posterior mean Eσ,α,β[h(bi)|ithsubject'sdata]. Without assuming σ,α,β being known, the empirical Bayes approach replaces them in the posterior mean by their estimates (σ^,α^,β^) so that h is estimated by
(6.19)h^=Eσ^,α^,β^[h(bi)∣ithsubject'sdata].
The expectation in (6.19) can be evaluated by the hybrid method, which we have used for likelihood calculation; see [80, 81]. First note that Laplace's approximation in (6.15) is based on the Taylor expansion
(6.20)li(b)≈li(b^i)+(b-b^i)l¨i(b^i)(b-b^i)T2.
Since the density function of bi given the ith risk class data is proportional to eli(b), it follows from (6.20) that the conditional distribution of bi given the data is approximately normal with mean b^i and covariance matrix Vi-1, where Vi=-l¨i(b^i). Hence for a subject with informative data (i.e., satisfying λmin(V)≥c), for whom we drop the subscript i in li, Vi, and b^i, we can use the preceding normal density to evaluate h^ in (6.19). Consider
(6.21)h(b)≈h(b^+V-1/2Z)≈h(b^)+h˙(b^)TV-1/2Z+12ZTV-1/2h¨(b^)TV-1/2Z,
where Z is a q-dimensional standard normal distribution. Using the preceding approximation gives
(6.22)h^=h(b^)+12V-1/2h¨(b^)TV-1/2
for the informative subject.When the eigenvalue criterion in the hybrid method fails (i.e., whenλmin(V)<c), (6.19) can be evaluated by using importance sampling for the Monte Carlo approximation
(6.23)h^≈∑h=1H[h(b(h))wh]∑h=1Hwh,
where {b(1),…,b(H)} are independent samples from the trial density λ(·) and the importance weights are wh=el(b(h))/λ(b(h)). Similar to the hybrid scheme presented in Section 6.3, λ(·) is chosen to be a mixture of the prior normal distribution with density ϕα and the posterior normal distribution N(b^i,[-l¨i(b^i∣σ,α,β)+εI]-1), where ε is a small positive number to ensure that the covariance matrix is invertible.


## 6.5. Information Sets and Selection of Variables for Prediction

We have assumed so far that the observations(Yi,t,xi,t,zi,t) are available at every 1≤t≤T, for all 1≤i≤n. In longitudinal data in biostatistics, however, there is often between-subject variations in the observation times. Lai et al. [82] recently addressed this difficulty by using a prediction approach that customizes the predictive model for an individual by choosing predictors from the individual's information set and using the data from all subjects whose information sets contain these predictors to estimate the parameters of the predictive model. Let Ti be the set of times when the ith subject is observed, and let ui,t be the subject's covariate vector at time t∈Ti. The marginal regression approach in longitudinal data analysis models E(Yi,t∣ui,t), t∈Ti, by using generalized linear models, under the independence assumption over time or other working dependence assumptions (such as exchangeability, or autoregression when Ti is regularly spaced) as in the GEE (generalized estimating equation) methodology introduced by Liang and Zeger [83]. It has been shown by Pepe and Anderson [84] that in ignorance of the true dependence structure over time and in the presence of time-varying covariates, the GEE parameter estimates may be biased unless the working independence assumption is used. In contrast to the marginal regression approach, the transitional modeling approach incorporates outcome and covariate histories to model E(Yi,t∣Yi,s,ui,τ;s∈Tiands<t,τ∈Tiandτ≤t) via generalized linear models, by estimating the model parameters from all subjects who are observed at all times t∈Ti. A compromise between marginal regression and transitional modeling is the partly conditional mean approach proposed by Pepe and Couper [85]. Instead of working with the subject-specific set Ti, this approach considers the set T of all possible times and introduces an indicator variable Ri(t) and a stratifying variable Ki(t) such that Ri(t)=1 (or 0) if the ith subject is observed (or not observed) at time t and Ki(t) takes the values 1,…,K to represent the stratum to which subject i belongs. Lai et al. [82, Section 2.1] have pointed out that partly conditional mean modeling can be viewed as a special case of a more general approach that chooses predictors from a subject's information set. This more general approach allows one to use ideas similar to those used in information criteria such as BIC or AIC to choose predictors from a subject's information set. More importantly, it provides a framework for Lai et al. [3] to extend and refine the dynamic EB approach as follows.To begin with, instead ofK strata in the partly conditional mean approach, we divide the subjects into K structurally similar subgroups. In many applications, subjects belonging to the same subgroup have similar observation times because of their structural similarity. For example, patients who have more serious ailments are monitored more frequently than others in a study cohort, causing the irregularity of observation times over different subgroups. Lai et al. [3] assume the cross-sectional dynamics (6.12) separately for each subgroup, that is, with μt and Y-t-j replaced by μt(k) and Y-t-j(k) for the kth subgroup, in which Y-s(k) is the sample average from all subjects (from the subgroup) who are observed at time s. Concerning (6.13) for μi,t with i belonging to the kth subgroup, we choose only predictors from the information set of subject i that are common to all subjects in the kth subgroup. Lai et al. [3] propose to use the BIC, details of which are given in [80, page 58-59], for the GLMM (6.13) with i belonging to the kth subgroup and Y-t-j replaced by Y-t-j(k), and give examples from the problem of predicting the batting performance of baseball players studied previously by Efron and Morris [86, 87] and Brown [88] to illustrate this point.


## 6.6. Evolutionary Credibility via GLMM (6.13)

In Section6.1 we have used the dynamic EB approach to elucidate the close connection between Kalman filtering in the linear evolutionary credibility theory in [62] and the linear mixed model (6.11). Section 6.2 extends (6.11) to the GLMM (6.13) to handle generalized linear models of the form (6.1). The usefulness of (6.1) for insurance claims data has been pointed out by Hsiao et al. [89] who note that these claims data often exhibit “excess zeros,” meaning that a large fraction of policies during a given year do not have claims. They propose to handle excess zeros by using a two-part model of the form
(6.24)Yit=IitZit,i=1,…,n;t=1,…,T,
where Iit=I{Yit≠0} and Zit has the conditional distribution of Yit given that Yit≠0. In this two-part model, if Yit is the claim size of the ith policy at time t, then Zit is the claim size given that the ith policy files a claim at time t.


### 6.6.1. Two-Part Nonlinear State-Space Models

The state-space model for linear evolutionary credibility in Section6.1 can be extended to the two-part model (6.24) as follows:
(6.25)Iit~Bernoulli(πit),Zit~Gamma(κ,μitκ)ifIit=1,
where κ is the shape parameter and link functions of πit=E(Iit) and μit=E(Zit) undergo AR(1) dynamics of the form
(6.26)logit(πit)=ρlogit(πi,t-1)+(1-ρ)logit(πi)+ηit,log(μit)=ρ′log(μi,t-1)+(1-ρ′)log(μi)+ϵit.
The function logit(π)=log(π/(1-π)) is the canonical link function for the exponential family of Bernoulli random variables. The ηit and ϵit in (6.26) are assumed to be independent normal random variables with mean 0, Var(ηit)=V, and Var(ϵit)=V'.Whereas Kalman filtering can be used to estimate the unobserved statesμt in the linear state-space model for evolutionary credibility in Section 6.1, the state-space model (6.25)-(6.26) is nonlinear, for which the Kalman filter is replaced by nonlinear filters to estimate the unobservable states πit and μit. These nonlinear filters are highly complex and require Markov chain Monte Carlo methods for their implementation when the model hyperparameters ρ,ρ',V,V',π1,…,πn,μ1,…,μn are known. Since they are unknown in practice, there is the additional complexity of estimating them, for which a stochastic EM approach implemented by Monte Carlo methods can be used to address the computational difficulties, as we have discussed in Section 4.1 for another application.


### 6.6.2. Two-Part Dynamic Modeling of Claims via GLMM

Since the Bernoulli and gamma distributions are special cases of the exponential family, we can apply the GLMM introduced in Section6.2 to obtain a more tractable two-part model of evolutionary credibility, replacing (6.26) by
(6.27)logit(πit)=ρ1logit(I-t-1)+μ'+ai,log(μit)=ρ2log(Z-t-1)+μ+bi,
in which ai~N(0,v') and bi~N(0,v) are random effects, I-s=n-1∑i=1n‍Iis, and
(6.28)Z-s=∑i:Iis≠0Zs∑i:Iis≠01
is the average claim size over all positive claims at time s. The parameters of the GLMM, which uses the cross-sectional means I-t-1 and Z-t-1 to substitute for the unobservable states πi,t-1 and μi,t-1 in (6.26), are μ,v,μ', and v'. As pointed out in Section 6.2, the GLMM can also incorporate additional covariates xit and zit as in (6.13) to improve the predictive power of the model.In the state-space model (6.25)-(6.26), the conditional distribution of the unobservable state vector (πit,μit) given the observed data Dt is difficult to represent analytically or algorithmically, unlike that in the linear Gaussian state-space model for which the conditional distribution is normal and its mean and covariance matrix have a closed-form recursive formula, that is, the Kalman filter. Instead of specifying the stochastic dynamics of the state vector precisely, West et al. [90] assume a conjugate family of prior distributions of the exponential family in the generalized linear model (6.1) so that the parameters of the conjugate family can be updated similarly to the Kalman filter. Lai and Sun [73] have compared the prediction performance of this approximation to the nonlinear HMM filters for the two-part state-space model (6.25)-(6.26) with that of the GLMM (6.27) in a simulation study, in which the data are generated by the two-part state-space model, and have found that the two methods perform similarly. An advantage of using GLMM to model evolutionary credibility over nonlinear state-space models or their approximations proposed in [90] is that GLMM has a relatively complete and easy-to-use methodology for likelihood inference and model selection; see [3, 80]. In particular, we can use model selection criteria such as BIC to choose the order p of an AR(p) model for μit, or to determine which of the variables in xit and zit should be included in the model to avoid overfitting.State-space models have traditionally been considered in the context of a singled-dimensional series. A distinctive feature of evolutionary credibility is that it involves many structurally similar time series, called “dynamic panel data” in econometrics, and we can exploit this feature to derive a dynamic empirical Bayes approach to evolutionary credibility. As pointed out in [91, page 335], traditional credibility theory can be called “empirical Bayes credibility” since it basically amounts to replacing the hyperparameters in the Bayes credibility (assuming Gaussian claims and a Gaussian prior) by their method-of-moment estimates. Extending this empirical Bayes approach to the Bayesian model of linear evolutionary credibility, which is a linear state-space model, yields the LMM in Section 6.1. Similarly, the GLMM in this section can be interpreted as a dynamic empirical Bayes approach to evolutionary credibility in generalized linear models.


## 6.6.1. Two-Part Nonlinear State-Space Models

The state-space model for linear evolutionary credibility in Section6.1 can be extended to the two-part model (6.24) as follows:
(6.25)Iit~Bernoulli(πit),Zit~Gamma(κ,μitκ)ifIit=1,
where κ is the shape parameter and link functions of πit=E(Iit) and μit=E(Zit) undergo AR(1) dynamics of the form
(6.26)logit(πit)=ρlogit(πi,t-1)+(1-ρ)logit(πi)+ηit,log(μit)=ρ′log(μi,t-1)+(1-ρ′)log(μi)+ϵit.
The function logit(π)=log(π/(1-π)) is the canonical link function for the exponential family of Bernoulli random variables. The ηit and ϵit in (6.26) are assumed to be independent normal random variables with mean 0, Var(ηit)=V, and Var(ϵit)=V'.Whereas Kalman filtering can be used to estimate the unobserved statesμt in the linear state-space model for evolutionary credibility in Section 6.1, the state-space model (6.25)-(6.26) is nonlinear, for which the Kalman filter is replaced by nonlinear filters to estimate the unobservable states πit and μit. These nonlinear filters are highly complex and require Markov chain Monte Carlo methods for their implementation when the model hyperparameters ρ,ρ',V,V',π1,…,πn,μ1,…,μn are known. Since they are unknown in practice, there is the additional complexity of estimating them, for which a stochastic EM approach implemented by Monte Carlo methods can be used to address the computational difficulties, as we have discussed in Section 4.1 for another application.


## 6.6.2. Two-Part Dynamic Modeling of Claims via GLMM

Since the Bernoulli and gamma distributions are special cases of the exponential family, we can apply the GLMM introduced in Section6.2 to obtain a more tractable two-part model of evolutionary credibility, replacing (6.26) by
(6.27)logit(πit)=ρ1logit(I-t-1)+μ'+ai,log(μit)=ρ2log(Z-t-1)+μ+bi,
in which ai~N(0,v') and bi~N(0,v) are random effects, I-s=n-1∑i=1n‍Iis, and
(6.28)Z-s=∑i:Iis≠0Zs∑i:Iis≠01
is the average claim size over all positive claims at time s. The parameters of the GLMM, which uses the cross-sectional means I-t-1 and Z-t-1 to substitute for the unobservable states πi,t-1 and μi,t-1 in (6.26), are μ,v,μ', and v'. As pointed out in Section 6.2, the GLMM can also incorporate additional covariates xit and zit as in (6.13) to improve the predictive power of the model.In the state-space model (6.25)-(6.26), the conditional distribution of the unobservable state vector (πit,μit) given the observed data Dt is difficult to represent analytically or algorithmically, unlike that in the linear Gaussian state-space model for which the conditional distribution is normal and its mean and covariance matrix have a closed-form recursive formula, that is, the Kalman filter. Instead of specifying the stochastic dynamics of the state vector precisely, West et al. [90] assume a conjugate family of prior distributions of the exponential family in the generalized linear model (6.1) so that the parameters of the conjugate family can be updated similarly to the Kalman filter. Lai and Sun [73] have compared the prediction performance of this approximation to the nonlinear HMM filters for the two-part state-space model (6.25)-(6.26) with that of the GLMM (6.27) in a simulation study, in which the data are generated by the two-part state-space model, and have found that the two methods perform similarly. An advantage of using GLMM to model evolutionary credibility over nonlinear state-space models or their approximations proposed in [90] is that GLMM has a relatively complete and easy-to-use methodology for likelihood inference and model selection; see [3, 80]. In particular, we can use model selection criteria such as BIC to choose the order p of an AR(p) model for μit, or to determine which of the variables in xit and zit should be included in the model to avoid overfitting.State-space models have traditionally been considered in the context of a singled-dimensional series. A distinctive feature of evolutionary credibility is that it involves many structurally similar time series, called “dynamic panel data” in econometrics, and we can exploit this feature to derive a dynamic empirical Bayes approach to evolutionary credibility. As pointed out in [91, page 335], traditional credibility theory can be called “empirical Bayes credibility” since it basically amounts to replacing the hyperparameters in the Bayes credibility (assuming Gaussian claims and a Gaussian prior) by their method-of-moment estimates. Extending this empirical Bayes approach to the Bayesian model of linear evolutionary credibility, which is a linear state-space model, yields the LMM in Section 6.1. Similarly, the GLMM in this section can be interpreted as a dynamic empirical Bayes approach to evolutionary credibility in generalized linear models.


## 6.7. Competing Risks and Retail Loans

We now return to the problem of predicting the default probabilities and loss given default of mortgages and corporate bonds that we have mentioned in Section1. The prediction problem is important for monitoring and management of a bank's credit risk. In 2001, the Basel Committee on Banking Supervision [92] proposed a new accord for risk evaluation, which was finalized in 2005 and is known as the Basel II Accord. Using the internal ratings-based (IRB) approach, the expected credit loss E(Lc) for a loan is decomposed into
(6.29)E(Lc)=EAD×LGD×PD,
where PD is the probability of default, LGD is the expected loss given default, which is expressed as a rate between 0 and 1, and EAD is the exposure at default. We first consider default probabilities of retail loans (such as personal loans, mortgages, credit cards, and car loans). “Retail” in banking means consumer-related activities, while “wholesale” means business-related loan transactions which include commercial loans, commercial real estate loans, and leases to businesses.Retail loans are typically grouped into several classes calledrisk buckets, within which the obligors can be regarded as having the same prior probability of default. Risk bucketing is tantamount to stratifying the obligors into classes so that the baseline covariate (such as credit score at the time of loan application) of an individual obligor can be replaced by the class membership. For retail loans, default means overdue payment for longer than 90 days. Thus, unlike a corporation that files for bankruptcy and the LGD of the loan is determined later by the recovery rate, a retail loan borrower can default and incur a loss to the lender while still holding the loan until it is foreclosed and taken off the balance sheet. Moreover, besides the possibility of default, another possible event that the lender faces is prepayment.


### 6.7.1. Competing Risks

Calhoun and Deng [93] regard default and prepayment as competing risks and introduce a multilogit regression model of the trinomial variable Y(t), taking the values 2, 1, and 0 that correspond to prepayment, foreclosure due to default, and still surviving at time t, respectively. In our setting that considers the trinomial variable Yij(t) of the ith obligor in the jth risk class, the multilogit regression model relates Yij(t) to a covariate vector uij(t) by
(6.30)P(Yij(t)=k∣Yij(s)=0fors<t,uij(t))=exp(αjk+βjkTuij(t))1+∑h=12exp(αjh+βjhTuij(t))
for k=1,2. Competing-risks models can be viewed as “marked” failure-time models, in which a mark is assigned to a failure event, indicating the type of the event that has occurred.Competing-risks data are also prevalent in cancer studies for which the response to treatment can be classified in terms of failure from disease processes and non-disease-related causes. For example, after treatments for leukemia and lymphoma, patients may either relapse or die in remission. In quality-of-life studies, the risks of the different failure types need to be considered for cost-effectivness analyses. Although Cox's proportional hazards model has routinely been used to analyze competing-risks data with covariates, what it models is thecause-specific hazard rate(6.31)λj(t;Z)=limh→0h-1P{t≤Y≤t+h,J=j∣Y≥t,Z},
where Z is a vector of covariates, Y is the failure time, and J is the cause of failure, with the causes labeled 1,2,…,j*. This approach treats each cause separately as the only cause for the observed failures that can be attributed to it and treats the observed failures due to other causes as censoring outcomes. Viewing competing risks as multiple endpoints each of which precludes the occurrence of the others, Fine and Gray [94] proposed to consider, instead of (6.31), the hazard function of the cumulative incidence due to cause j:
(6.32)λ~j(t;Z)=limh→0h-1P[t≤Y≤t+h,J=j∣{Y≥t}∪{Y<t,J≠j},Z]=-ddtlog(1-Fj(t;Z)),
where Fj(t;Z) is the cumulative incidence function(6.33)Fj(t;Z)=P{Y≤t,J=j∣Z}
introduced by Gray [95]. Note that 1-Fj(t-;Z) = 1-P{Y<t,J=j∣Z} = P{Y≥t∣Z} + P{Y<t,J≠j∣Z}. In the absence of covariates, the cumulative incidence function can be estimated by the Aalen-Johansen estimator [27, 28] based on a finite-state Markov chain model of competing risks; see the survey by Klein and Shu [96]. Although the proportional hazards model λj(t;Z)=λj0(t)eβTZ can readily be extended to the hazard functions λ~j of the competing risks, regression analysis of λ~j under the model λ~j(t;Z)=λ~j0(t)exp(βTZ) becomes much more difficult. Three approaches to this problem have been developed. First, Fine and Gray [94] used the following modification of risk sets to extend partial likelihood to competing-risks data. The risk set at time t includes the subjects that are naturally at risk and those who have failed prior to time t due to other competing risks and are not censored at time t. A subject who has failed prior to time t due to other competing risks is weighted by the probability that it is not censored at time t given that it has not been censored at the failure time. The second method, proposed by Scheike et al. [97], uses weighted counting processes as responses in a binomial regression model for cumulative incidence functions. Both methods assume a common censoring distribution. Their performance is highly sensitive to how well the censoring distribution can be estimated. The third method, proposed by Klein and Andersen [98], uses jackknife pseudovalues to model the cumulative incidence functions directly. These pseudo-values are used in a generalized estimating equation to obtain estimates on model parameters, and Graw et al. [99] have recently shown this method to be consistent and asymptotically normal.An important difference between the retail loans dataYij(t) in (6.30) and those in cancer biostatistics is that the latter is subject to right-censoring due to patient withdrawal whereas the former remain in the lender's balance sheet until prepayment or default leading to foreclosure. Right-censoring results in much greater complexity than multilogit regression in modeling competing risks.


### 6.7.2. GLMM for Loan Default Probabilities

Calhoun and Deng [93] fitted (6.30) to about 1.3 million conforming fixed-rate and adjustable-rate residential mortgages originated between 1979 and 1993. Their covariate vector includes loan age, dummy variables for origination years, the loan-to-value ratio (LTV), seasonality, occupancy status of the property, the original loan size, the ratio of the 10-year Constant Maturity Treasury rate to the 1-year Constant Maturity Treasury rate, mortgage premium value, and the probability of negative equity based on estimated drift and volatility of housing prices. Instead of the multilogit model (6.30) that generalizes logistic regression, Lai et al. [100] use the GLMM in Section 6.2 to model the trinomial variable Yij(t) that represents continuing payment, or default, or prepayment at time t of the ith obligor in the jth risk class. Whereas [93] does not divide the borrowers into risk classes, the statistical analysis of subprime mortgages in [100] conforms to risk bucketing that is commonly used for risk management of retail loans and is based on loan-level mortgage data consisting of 10,000 first-lien, 2–28 adjustable-rate subprime loans originated between January 2004 and June 2006; these loans are randomly selected from a database maintained by LoanPerformance, Inc., a mortgage information and analytics provider. The 10,000 loans are grouped into I=5 risk classes according to the FICO score when a loan was originated; the class with the lowest FICO scores is labeled 1 and the highest labeled 5. These data involve two time scales, calendar time t, and loan age s, which also arise in biostatistics in the context of time-sequential clinical trials with failure-time endpoints and staggered entry, where s=(t-τ)+ and τ is the entry time of the patient into the study (or in the case of mortgages the initiation date of the loan); see [101]. The case s=0 denotes that the loan has not been originated at calendar time t.Lai et al. [100] use the link function
(6.34)ηijt(k)=log(P{Yij(t)=k∣Yij(t-1)=0,uij(t)}P{Yij(t)=0∣Yij(t-1)=0,uij(t)}),k=1,2,
to express (6.30) in the generalized linear form, which can then be extended to the GLMM in Section 6.2:
(6.35)ηijt(k)=ρ(k)log(Y-t-1(k)Y-t-1(0))+aj(k)+β(k)TXijt+bj(k)TZijt,
where Xijt and Zijt are covariates that are available prior to t and Y-t-1(k) is the cross-sectional mean of the indicator variables I{Yij(t-1)=h,Yij(t-2)=0} for k=0,1,2. If the loan has not been initiated at time t-1, then the conditional probabilities in (6.34) are taken to be P{Yij(t)=·∣uij(t)} and the cross-sectional term involving Y-t-1(k)/Y-t-1(0) is also removed from (6.35). The covariate vector Xijt includes the loan age, occupancy status, initial mortgage rate, purpose (for purchase or refinance), a dummy variable for full documentation, or otherwise (important for reflecting the lender's standards in issuing the subprime loan), and Zt includes mortgage premium value and equity position as defined in [100], which involve the macroeconomic variables 1-year Constant Maturity Treasury rate and zip-code-based house price index, together with the loan's size, unpaid rate, mortgage rate, and LTV. These data do not contain the loan's payment history; ideally Yij(t) should also indicate if the payment is overdue for longer than 90 days even when the loan is not foreclosed, but such information is not available in the data obtained from LoanPerformance, only showing the trinomial variable if the loan is prepaid, foreclosed, or active at time t. The advantages of (6.35) over the multilogit model (6.30) in predicting defaults are discussed in [102].


### 6.7.3. Loss Given Default

Lai and Wong [103, Section 5.2] have described how GLMM can be used to model the LGD of retail loans. They first review how LGD is calculated for corporate loans and bonds in the software package LossCalc of Moody's KMV, which is built from a training sample of recovery data of the defaulted corporate loans, bonds, and preferred stock. The training sample is used to fit a linear regression model of y=Φ-1(Fα,β(RR)) on nine predictors chosen from a set of regressions using a year-by-year out-of-sample prediction error criterion, where RR is the recovery rate implied by the market value (bid-side market quote) one month after default, Φ is the standard normal distribution function, and Fα,β is the distribution function of the Beta(α,β) distribution. The transformation Φ-1∘Fα,β is used to convert a beta random variable into a normal one so that linear regression can be applied. As shown in [103, page 237], a better approach is to fit the beta regression model, which is a generalized linear model.As pointed out in [103], unlike corporate loans are still quite liquid after default, personal loans, auto loans, and mortgage loans lack such liquidity and the usual recovery process involves (a) in-house departments, (b) selling off the debt at some fraction of the face value, and (c) collection or foreclosure agencies which take a fixed percentage of the recovered loan as commission. To describe simply the recovery process, the bank first tries (a) and switches to (b) or (c) if (a) is unsuccessful after, say, three months. In the case of mortgage loans, (b) means selling the mortgage to another mortgage company (e.g., Fannie Mae or Freddie Mac in the US), while (c) means foreclosure of the mortgage and repossession and sale of the property. The shortfall is the residual loan value plus planned interest minus the sales price of the mortgage or the repossessed property. For occasional missed payments, the “defaulted” obligors who have forgotten to pay or have been away usually pay within one month after they are contacted by the in-house department. Such a scenario corresponds to complete recovery and leads to an atom at 0 in the marginal distribution of the economic LGD. Lai and Wong [103] point out the importance of including the atom 0 and demonstrate the advantages of using beta regression over linear regression of the transformed recovery rate on covariates. They also note that the Basel II Accord requires LGD estimates to be grounded in historical recoveries and to incorporate business cycles since at times of economic downturn both PD and LGD tend to increase; the database should ideally cover at least one business cycle and must not be shorter than five years. They propose the following GLMM to estimate LGD using these data.LetDs denote the set of obligors who default in year s or before and who have available credit and collateral information in year s-1, with s=t,t-1,…,t-4. For i∈Ds, let Yi(s) be the economic LGD in year s. As noted above, Yi(s) has an atom at 0. Let πi(s)=P{Yi(s)>0}, which can be modeled by
(6.36)logit(πi(s))=ξi+ηTxi(s-1)+γTf(s),
where xi(s-1) is a vector of attributes of the obligor's creditworthiness and collateral in year s-1 and f(s) is a vector of macroeconomic factors in year s. The conditional density of Yi(s)∣Yi(s)>0 is assumed to belong to the exponential family of beta density functions
(6.37)f(y∣μi(s),ϕ)=Γ(ϕ)Γ(μi(s)ϕ)Γ((1-μi(s))ϕ)yμi(s)ϕ-1(1-y)(1-μi(s))ϕ-1
for 0<y<1, where μi(s)=E(Yi(s)∣Yi(s)>0),
(6.38)logit(μi(s))=ζi+κTxi(s-1)+νTf(s),Var(Yi(s)∣Yi(s)>0)=μi(s)(1-μi(s))/(1+ϕ), and ξi are i.i.d. normal and are independent of ζi which are also i.i.d. normal. This GLMM includes loan-specific random effects ξi and ζi in (6.36) and (6.38) and is similar to the two-part model in Section 6.6. It enables one to predict directly LGDi(t+1) of the ith obligor for year t+1 using xi(t) and the forecast of f(t+1), noting that LGDi(t+1)=πi(t+1)E(Yi(t+1)∣Yi(t+1)>0).


## 6.7.1. Competing Risks

Calhoun and Deng [93] regard default and prepayment as competing risks and introduce a multilogit regression model of the trinomial variable Y(t), taking the values 2, 1, and 0 that correspond to prepayment, foreclosure due to default, and still surviving at time t, respectively. In our setting that considers the trinomial variable Yij(t) of the ith obligor in the jth risk class, the multilogit regression model relates Yij(t) to a covariate vector uij(t) by
(6.30)P(Yij(t)=k∣Yij(s)=0fors<t,uij(t))=exp(αjk+βjkTuij(t))1+∑h=12exp(αjh+βjhTuij(t))
for k=1,2. Competing-risks models can be viewed as “marked” failure-time models, in which a mark is assigned to a failure event, indicating the type of the event that has occurred.Competing-risks data are also prevalent in cancer studies for which the response to treatment can be classified in terms of failure from disease processes and non-disease-related causes. For example, after treatments for leukemia and lymphoma, patients may either relapse or die in remission. In quality-of-life studies, the risks of the different failure types need to be considered for cost-effectivness analyses. Although Cox's proportional hazards model has routinely been used to analyze competing-risks data with covariates, what it models is thecause-specific hazard rate(6.31)λj(t;Z)=limh→0h-1P{t≤Y≤t+h,J=j∣Y≥t,Z},
where Z is a vector of covariates, Y is the failure time, and J is the cause of failure, with the causes labeled 1,2,…,j*. This approach treats each cause separately as the only cause for the observed failures that can be attributed to it and treats the observed failures due to other causes as censoring outcomes. Viewing competing risks as multiple endpoints each of which precludes the occurrence of the others, Fine and Gray [94] proposed to consider, instead of (6.31), the hazard function of the cumulative incidence due to cause j:
(6.32)λ~j(t;Z)=limh→0h-1P[t≤Y≤t+h,J=j∣{Y≥t}∪{Y<t,J≠j},Z]=-ddtlog(1-Fj(t;Z)),
where Fj(t;Z) is the cumulative incidence function(6.33)Fj(t;Z)=P{Y≤t,J=j∣Z}
introduced by Gray [95]. Note that 1-Fj(t-;Z) = 1-P{Y<t,J=j∣Z} = P{Y≥t∣Z} + P{Y<t,J≠j∣Z}. In the absence of covariates, the cumulative incidence function can be estimated by the Aalen-Johansen estimator [27, 28] based on a finite-state Markov chain model of competing risks; see the survey by Klein and Shu [96]. Although the proportional hazards model λj(t;Z)=λj0(t)eβTZ can readily be extended to the hazard functions λ~j of the competing risks, regression analysis of λ~j under the model λ~j(t;Z)=λ~j0(t)exp(βTZ) becomes much more difficult. Three approaches to this problem have been developed. First, Fine and Gray [94] used the following modification of risk sets to extend partial likelihood to competing-risks data. The risk set at time t includes the subjects that are naturally at risk and those who have failed prior to time t due to other competing risks and are not censored at time t. A subject who has failed prior to time t due to other competing risks is weighted by the probability that it is not censored at time t given that it has not been censored at the failure time. The second method, proposed by Scheike et al. [97], uses weighted counting processes as responses in a binomial regression model for cumulative incidence functions. Both methods assume a common censoring distribution. Their performance is highly sensitive to how well the censoring distribution can be estimated. The third method, proposed by Klein and Andersen [98], uses jackknife pseudovalues to model the cumulative incidence functions directly. These pseudo-values are used in a generalized estimating equation to obtain estimates on model parameters, and Graw et al. [99] have recently shown this method to be consistent and asymptotically normal.An important difference between the retail loans dataYij(t) in (6.30) and those in cancer biostatistics is that the latter is subject to right-censoring due to patient withdrawal whereas the former remain in the lender's balance sheet until prepayment or default leading to foreclosure. Right-censoring results in much greater complexity than multilogit regression in modeling competing risks.


## 6.7.2. GLMM for Loan Default Probabilities

Calhoun and Deng [93] fitted (6.30) to about 1.3 million conforming fixed-rate and adjustable-rate residential mortgages originated between 1979 and 1993. Their covariate vector includes loan age, dummy variables for origination years, the loan-to-value ratio (LTV), seasonality, occupancy status of the property, the original loan size, the ratio of the 10-year Constant Maturity Treasury rate to the 1-year Constant Maturity Treasury rate, mortgage premium value, and the probability of negative equity based on estimated drift and volatility of housing prices. Instead of the multilogit model (6.30) that generalizes logistic regression, Lai et al. [100] use the GLMM in Section 6.2 to model the trinomial variable Yij(t) that represents continuing payment, or default, or prepayment at time t of the ith obligor in the jth risk class. Whereas [93] does not divide the borrowers into risk classes, the statistical analysis of subprime mortgages in [100] conforms to risk bucketing that is commonly used for risk management of retail loans and is based on loan-level mortgage data consisting of 10,000 first-lien, 2–28 adjustable-rate subprime loans originated between January 2004 and June 2006; these loans are randomly selected from a database maintained by LoanPerformance, Inc., a mortgage information and analytics provider. The 10,000 loans are grouped into I=5 risk classes according to the FICO score when a loan was originated; the class with the lowest FICO scores is labeled 1 and the highest labeled 5. These data involve two time scales, calendar time t, and loan age s, which also arise in biostatistics in the context of time-sequential clinical trials with failure-time endpoints and staggered entry, where s=(t-τ)+ and τ is the entry time of the patient into the study (or in the case of mortgages the initiation date of the loan); see [101]. The case s=0 denotes that the loan has not been originated at calendar time t.Lai et al. [100] use the link function
(6.34)ηijt(k)=log(P{Yij(t)=k∣Yij(t-1)=0,uij(t)}P{Yij(t)=0∣Yij(t-1)=0,uij(t)}),k=1,2,
to express (6.30) in the generalized linear form, which can then be extended to the GLMM in Section 6.2:
(6.35)ηijt(k)=ρ(k)log(Y-t-1(k)Y-t-1(0))+aj(k)+β(k)TXijt+bj(k)TZijt,
where Xijt and Zijt are covariates that are available prior to t and Y-t-1(k) is the cross-sectional mean of the indicator variables I{Yij(t-1)=h,Yij(t-2)=0} for k=0,1,2. If the loan has not been initiated at time t-1, then the conditional probabilities in (6.34) are taken to be P{Yij(t)=·∣uij(t)} and the cross-sectional term involving Y-t-1(k)/Y-t-1(0) is also removed from (6.35). The covariate vector Xijt includes the loan age, occupancy status, initial mortgage rate, purpose (for purchase or refinance), a dummy variable for full documentation, or otherwise (important for reflecting the lender's standards in issuing the subprime loan), and Zt includes mortgage premium value and equity position as defined in [100], which involve the macroeconomic variables 1-year Constant Maturity Treasury rate and zip-code-based house price index, together with the loan's size, unpaid rate, mortgage rate, and LTV. These data do not contain the loan's payment history; ideally Yij(t) should also indicate if the payment is overdue for longer than 90 days even when the loan is not foreclosed, but such information is not available in the data obtained from LoanPerformance, only showing the trinomial variable if the loan is prepaid, foreclosed, or active at time t. The advantages of (6.35) over the multilogit model (6.30) in predicting defaults are discussed in [102].


## 6.7.3. Loss Given Default

Lai and Wong [103, Section 5.2] have described how GLMM can be used to model the LGD of retail loans. They first review how LGD is calculated for corporate loans and bonds in the software package LossCalc of Moody's KMV, which is built from a training sample of recovery data of the defaulted corporate loans, bonds, and preferred stock. The training sample is used to fit a linear regression model of y=Φ-1(Fα,β(RR)) on nine predictors chosen from a set of regressions using a year-by-year out-of-sample prediction error criterion, where RR is the recovery rate implied by the market value (bid-side market quote) one month after default, Φ is the standard normal distribution function, and Fα,β is the distribution function of the Beta(α,β) distribution. The transformation Φ-1∘Fα,β is used to convert a beta random variable into a normal one so that linear regression can be applied. As shown in [103, page 237], a better approach is to fit the beta regression model, which is a generalized linear model.As pointed out in [103], unlike corporate loans are still quite liquid after default, personal loans, auto loans, and mortgage loans lack such liquidity and the usual recovery process involves (a) in-house departments, (b) selling off the debt at some fraction of the face value, and (c) collection or foreclosure agencies which take a fixed percentage of the recovered loan as commission. To describe simply the recovery process, the bank first tries (a) and switches to (b) or (c) if (a) is unsuccessful after, say, three months. In the case of mortgage loans, (b) means selling the mortgage to another mortgage company (e.g., Fannie Mae or Freddie Mac in the US), while (c) means foreclosure of the mortgage and repossession and sale of the property. The shortfall is the residual loan value plus planned interest minus the sales price of the mortgage or the repossessed property. For occasional missed payments, the “defaulted” obligors who have forgotten to pay or have been away usually pay within one month after they are contacted by the in-house department. Such a scenario corresponds to complete recovery and leads to an atom at 0 in the marginal distribution of the economic LGD. Lai and Wong [103] point out the importance of including the atom 0 and demonstrate the advantages of using beta regression over linear regression of the transformed recovery rate on covariates. They also note that the Basel II Accord requires LGD estimates to be grounded in historical recoveries and to incorporate business cycles since at times of economic downturn both PD and LGD tend to increase; the database should ideally cover at least one business cycle and must not be shorter than five years. They propose the following GLMM to estimate LGD using these data.LetDs denote the set of obligors who default in year s or before and who have available credit and collateral information in year s-1, with s=t,t-1,…,t-4. For i∈Ds, let Yi(s) be the economic LGD in year s. As noted above, Yi(s) has an atom at 0. Let πi(s)=P{Yi(s)>0}, which can be modeled by
(6.36)logit(πi(s))=ξi+ηTxi(s-1)+γTf(s),
where xi(s-1) is a vector of attributes of the obligor's creditworthiness and collateral in year s-1 and f(s) is a vector of macroeconomic factors in year s. The conditional density of Yi(s)∣Yi(s)>0 is assumed to belong to the exponential family of beta density functions
(6.37)f(y∣μi(s),ϕ)=Γ(ϕ)Γ(μi(s)ϕ)Γ((1-μi(s))ϕ)yμi(s)ϕ-1(1-y)(1-μi(s))ϕ-1
for 0<y<1, where μi(s)=E(Yi(s)∣Yi(s)>0),
(6.38)logit(μi(s))=ζi+κTxi(s-1)+νTf(s),Var(Yi(s)∣Yi(s)>0)=μi(s)(1-μi(s))/(1+ϕ), and ξi are i.i.d. normal and are independent of ζi which are also i.i.d. normal. This GLMM includes loan-specific random effects ξi and ζi in (6.36) and (6.38) and is similar to the two-part model in Section 6.6. It enables one to predict directly LGDi(t+1) of the ith obligor for year t+1 using xi(t) and the forecast of f(t+1), noting that LGDi(t+1)=πi(t+1)E(Yi(t+1)∣Yi(t+1)>0).


## 6.8. GLMM for Modeling Frailty and Contagion of Credit Portfolios

Similarly to the applications to evolutionary credibility and mortgage loans in the preceding two sections, the dynamic EB approach via GLMM can be applied to provide a new model for correlated default intensities of multiple firms in a portfolio of corporate bonds or a multiname credit derivative. This model can incorporate both frailty and contagion via the covariates of the GLMM and is considerably easier to implement. The details are given in [3, 100]. Note that the use of a mixture distribution for firm-specific random effects also appears in Vasicek's mixed binomial model (3.11) for the number of defaults in a loan portfolio up to a horizon T. It is shown in [3] that the prediction performance of the nonlinear state-space model of Duffie et al. [46] is similar to that of the much simpler GLMM model of dynamic frailty even when the data are actually generated by the nonlinear state-space model (4.2).


## 7. Conclusion

This paper shows the methodological connections between (a) default modeling and prediction of mortgage loans and portfolios of corporate bonds, a topic of timely interest in the wake of the financial crisis following the events in 2007 and 2008 described in Section1, (b) insurance rate-making in evolutionary credibility theory, and (c) analysis of longitudinal data in biostatistics via GLMM. The underlying statistical principle linking these disparate fields is dynamic empirical Bayes developed in [3] on the foundational works on empirical Bayes by Robbins [1] and Stein [2]. Although survival analysis may seem to be a more natural candidate from biostatistics to connect to bond or loan default—and indeed we have provided such connections in Section 2—a major difference between survival data in biostatistics and default data in finance lies in the event rate, for which the latter typically has very low rate. In fact, importance sampling is often needed to simulate the “rare events” of loan default in managing portfolio credit risk; see [104]. Another difference lies in the objectives of the statistical analysis. In biostatistics the goal is to compare survival patterns of different groups (e.g., one receiving an experimental treatment and the other receiving standard of care in cancer patients), whereas in finance the goal is to predict default and to use the prediction for pricing, hedging, and risk management. It should be noted that the complications in the analysis of survival data in biostatistics are mostly caused by censoring. On the other hand, with loans, censoring is not an issue if one introduces competing risks, as shown in Section 6.7.



---

*Source: 832175-2012-12-23.xml*