Modelling the Time to Write-Off of Non-Performing Loans Using a Promotion Time Cure Model with Parametric Frailty

Abstract

Modelling the outcome after loan default is receiving increasing attention, and survival analysis is particularly suitable for this purpose due to the likely presence of censoring in the data. In this study, we suggest that the time to loan write-off may be influenced by latent competing risks, as well as by common, unobservable drivers, such as the state of the economy. We therefore expand on the promotion time cure model and include a parametric frailty parameter to account for common, unobservable factors and for possible observable covariates not included in the model. We opt for a parametric model due to its interpretability and analytical tractability, which are desirable properties in bank risk management. Both a gamma and inverse Gaussian frailty parameter are considered for the univariate case, and we also consider a shared frailty model. A Monte Carlo study demonstrates that the parameter estimation of the models is reliable, after which they are fitted to a real-world dataset in respect of large corporate loans in the US. The results show that a more flexible hazard function is possible by including a frailty parameter. Furthermore, the shared frailty model shows potential to capture dependence in write-off times within industry groups.

1. Introduction

Risk management in the credit risk environment has advanced considerably over the last few decades, and during the same time, statistical models developed for this purpose have grown in variety and sophistication. Initially, statistical modelling focused on the assessment of a counterparty’s creditworthiness (i.e., default risk), but the modelling of the recovery process and outcome after default are receiving increasing attention. Two main outcomes or resolutions are possible after default: either the outstanding amount is fully recovered or the loan is written off. Modelling the time to write-off of a defaulted loan is valuable to a bank, as the time on book of such loans not only affects its operational decisions but also its profitability [1]. Survival analysis has enjoyed increasing popularity in the credit risk environment for modelling of the time to default and, more recently, to model the time to resolution of a defaulted loan (see, e.g., [2,3,4,5]). Since write-offs are typically associated with longer resolution times than full recoveries, the number of write-offs is usually under-represented in samples of defaulted loans [6]. This makes survival analysis particularly suitable for modelling the time to write-off because the status at the end of the observation period is often recorded as “not resolved”, i.e., censoring is present.

With the fast-growing body of literature on the use of survival analysis (and specifically the Cox proportional hazard (CPH) model) to model default risk, an exploration of the use of the methodology to describe the post-default outcome of a loan is a natural outflow. Several authors have made use of Cox regression to demonstrate that loan and obligor characteristics, as well as macroeconomic factors, drive the resolution time of defaulted loans (see e.g., [3,4,7,8]).

Recently, both mixture cure and promotion time cure (PTC) models have found applications in the credit risk environment (see e.g., [9,10,11,12]). Recognising that there is a proportion of defaulted loans that is not “susceptible” to a full recovery, de Oliveira and Louzada [5] modelled the recovery process of defaulted loans using a PTC model. The model considers the full repayment of a defaulted loan as the event of interest and allows for latent factors that may impact the post-default resolution outcome.

The time to write-off or full recovery may, in addition to latent competing risks, also be influenced by common, unobservable drivers, such as the state of the economy [3]. In this study, we propose an expansion of the PTC model by including a frailty parameter to control for unobserved heterogeneity within our population of defaulted loans. Underlying the post-default recovery process is the macroeconomic environment, which has a systematic effect on all borrowers’ ability to resolve their defaulted obligations. Although some of these are observable and can be included as covariates, we propose a frailty parameter that can account for common, unobservable factors, as well as for possible observable covariates not included in the model. We consider both the univariate and the shared frailty parameter. Global regulatory standards, specifically those set by the Basel Accords, require banks to have a comprehensive framework for model risk management, which includes ongoing monitoring and assessment of model performance. We demonstrate how inclusion of a frailty parameter in the PTC model results in improved model performance when estimating the time to write-off of a non-performing loan.

Our paper offers the following contributions. First, we formulate the PTC model with gamma frailty and derive its likelihood function. Second, we compare the performance of this model to that of the standard PTC model, as well as to a PTC model that includes an inverse Gaussian frailty parameter, as proposed by Barriga et al. [13]. For the comparison, we first perform a simulation study, then fit the models to both a simulated and a real-world dataset relating to the resolution outcomes of defaulted corporate loans. Third, we extend the PTC model to include a shared frailty parameter. We derive the likelihood function of this model and evaluate its performance on the same real-world dataset.

The probability of write-off (or recovery) is a key component within two-step loss given default (LGD) models, where the post-default outcome probability and loss severity are characterised separately (see, e.g., [14]). Improving on the accuracy of LGD models should therefore enable a bank to attain more accurate estimates of this risk parameter for the determination of credit loss provisions, which should, in turn, result in more stable earnings. It will also enable the bank to more accurately assess its regulatory and/or economic capital requirements, which are crucial for a bank’s financial stability. Obtaining a more accurate estimate of the time to write-off will, furthermore, enable the bank to better manage its cash flow requirements and reduce the liquidity risk to which it is exposed.

The remainder of this article is structured as follows. In the next section, an overview of survival modelling is provided, including cure models and their application in the credit risk environment. The formulation of the proposed model, namely the PTC model with parametric frailty, follows in Section 3. In Section 4, the likelihood function of this model is derived, and the finite sample performance of maximum likelihood estimators (MLEs) of the frailty models is analysed through a Monte Carlo simulation study. We apply our model to both a simulated dataset and a real-world loan loss dataset in Section 5. In Section 6, we describe the extension of the PTC model with a shared frailty parameter, including an application to a real-world dataset. We provide some concluding remarks and avenues for future research in Section 7.

2. Background

Survival modelling is a statistical method used to analyse and predict the time until an event of interest occurs. The event of interest is often referred to as the “failure” or “death” event, and the time until this event is referred to as the “survival time”. The event of interest in our context is the write-off of a defaulted loan. In general, the survival function is defined as $S\left(t\right)=P(T>t)$, where T is the random variable representing the time to the event, and t is a specific point in time.

Within the resolution process of defaulted loans, random right censoring is normally present due to the fixed duration of the observation period. Not all subjects experience the event of interest, i.e., write-off, within that time period. In right censoring, there are two latent random variables: T and C, where, T is the time of failure, and C is the censoring time. The random observed variable is $(Y,\delta )$, where

$$Y=\min (T,C)$$

and

$$\begin{array}{cc}\hfill \delta =& \left\{\begin{array}{cc}1\hfill & \mathrm{if}T\le C\hfill \\ 0\hfill & \mathrm{if}T>C\hfill \end{array}\right..\hfill \end{array}$$

The simplest way of modelling the time to write-off is to assign a distribution with positive support to T, for example, the gamma or Weibull distributions.

However, for credit risk modelling, specifying $S\left(t\right)$ unconditionally without taking covariates into account is too simplistic, and the CPH model has therefore dominated.

The CPH model, which was proposed in the seminal paper of Cox [15], assumes that the hazard function is a function of a set of covariates and that the effect of the covariates on the hazard function is proportional over time.

The hazard function $h\left(t\right)=f\left(t\right)/S\left(t\right)$ of a survival process represents the instantaneous risk of the event of interest at a given time, where f is the density of T. For the CPH model, the hazard function, conditional on a vector of covariates $\mathbf{x}$ = ${\left(x_1,x_2,\dots ,x_m\right)}^T,$ has the form

$$h\left(t\right|\mathbf{x})=h_0\left(t\right)e^{{\mathbf{\beta }}^T\mathbf{x}},$$

(1)

where $h_0\left(t\right)$ is the baseline hazard function, and $\mathbf{\beta }$ = ${\left({\beta }_1,{\beta }_2,\dots ,{\beta }_m\right)}^T$ is the vector of unknown regression parameters associated with $\mathbf{x}$. The conditional survival function is then given by

$$S\left(t\right|\mathbf{x})=P\left(T>t|\mathbf{x}\right)=S_0{\left(t\right)}^{\exp \left({\mathbf{\beta }}^T\mathbf{x}\right)},$$

where $S_0\left(t\right)=S\left(t\right|\mathbf{x}=\mathbf{0})$ is the baseline survival function.

The earliest studies to make use of the CPH model to model default risk include those of Banasik et al. [16], Stepanova and Thomas [17] and Cao et al. [18]; subsequently, its application to credit risk problems has increased markedly. Several authors, including Betz et al. [3], Malwandla et al. [4] and Joubert et al. [8], have since made use of Cox regression to model the post-default outcome of non-performing loans.

Notwithstanding the popularity of the CPH model, cure models have recently enjoyed increasing attention within the credit risk domain. This is particularly true for the modelling of the probability of default, which, in many high-quality loan portfolios, represents a rare event (see, e.g., [9,10,11]). Cure models, as all survival models, had their origin in medicine, where Boag [19] first suggested that there exists a certain proportion of cancer patients who will never experience the event of interest, i.e., they can be considered cured. By not accounting for the cure proportion, the estimates of the survival of susceptible subjects may be biased [20].

In early work by Berkson and Gage [21] the following simple model was proposed:

$$S\left(t\right)={\pi }_0{S}^c\left(t\right)+(1-{\pi }_0)S^s\left(t\right),$$

(2)

where $S^c$ and $S^s$ are the survival functions for non-susceptible and susceptible individuals, respectively, and ${\pi }_0$ is the proportion of individuals considered non-susceptible, also known as the cure proportion. Note that $S^c\left(t\right)=1$ for all $t\ge 0$, as the individuals being modelled are not susceptible to the event of interest and therefore have a probability of survival equal to one.

The model in (2) was extended by Farewell [22] to allow for covariates. Suppose another set of covariates $\mathbf{y}$ that may be equal to $\mathbf{x}$ is observed. The conditional survival function for a mixture cure model is defined as:

$$S\left(t\right|\mathbf{x},\mathbf{y})={\pi }_0\left(\mathbf{x}\right)+\left(1-{\pi }_0\left(\mathbf{x}\right)\right)S^s\left(t\right|\mathbf{y}),$$

(3)

where ${\pi }_0\left(\mathbf{x}\right)$ is commonly modelled using a logistic regression model, and $S^s\left(t\right|\mathbf{y})$ is modelled by a CPH or accelerated failure time model. For a detailed overview of mixture cure models, the reader is referred to the review paper by Amico and van Keilegom [23].

The second type of cure model, known as the promotion time cure (PTC) model or, alternatively, the bounded cumulative hazard model, was first introduced by Yakovlev et al. [24]. The PTC model has an intuitive biological interpretation: an individual cancer patient has a number of carcinogenic cells (denoted by N) left active after initial treatment, which have the potential to activate over time and produce a significant relapse of cancer. Chen et al. [25] further proposed that an individual is at risk of relapse if he or she is exposed to at least one of the so-called latent factors, and relapse (or failure) is observed when one (or more) of these latent factors becomes activated. If the individual is not considered to be at risk, he or she is considered cured.

Formally, let N be the unobservable number of causes of the event of interest (write-off in our case) for a given subject in the population. For a subject for which $N\ge 1$, the event time is defined as

$$T=min\left(Z_1,Z_2,\dots ,Z_N\right),$$

where $Z_1,Z_2,\dots ,Z_N$ represents the time required for the jth unobservable cause to be realised. A subject is cured if $N=0$, which implies that $T=\infty$, i.e., write-off will not happen. The set of observed covariates is denoted by $\mathbf{x}={\left(1,x_1,x_2,\dots ,x_m\right)}^T$. Rodrigues et al. [26] demonstrate that the long-term survival function satisfies the proportional hazards property if and only if the number of risks follows a Poisson distribution. Assume that N follows a Poisson distribution with a mean of $\theta \left(\mathbf{x}\right)>0$, where $\theta \left(\mathbf{x}\right)=g\left({\mathbf{\beta }}^T\mathbf{x}\right)$ for a strictly increasing function g (see e.g., Zhang and Peng [27] and Portier et al. [28]). Furthermore, assume that $Z_1,Z_2,\dots ,Z_N$ are i.i.d. from a distribution $F_{\mathbf{\eta }}$ with parameter $\mathbf{\eta }$, which could be a vector and independent of N. We then arrive at the well-known PTC model,

$$S_{\mathbf{\eta }}\left(t\right|\mathbf{x})=P\left(T>t|\mathbf{x}\right)=e^{-\theta \left(\mathbf{x}\right)F_{\mathbf{\eta }}\left(t\right)}.$$

(4)

For a short derivation of (4), the interested reader is referred to the work of Chen et al. [25]. It is clear from (4) that the cure fraction (i.e., for those not susceptible to write-off) is given by

$$\underset{t\to \infty }{lim}{S}_{\mathbf{\eta }}\left(t\right|\mathbf{x})=e^{-\theta \left(\mathbf{x}\right)}={\pi }_0\left(\mathbf{x}\right).$$

To obtain a similar expression as in (3), the survival function in (4) can be decomposed to obtain

$$S_{\mathbf{\eta }}\left(t\right|\mathbf{x})={\pi }_0\left(\mathbf{x}\right)+\left(1-{\pi }_0\left(\mathbf{x}\right)\right)S_{\mathbf{\eta }}^s\left(t\right|\mathbf{x}),$$

where

$$S_{\mathbf{\eta }}^s\left(t\right|\mathbf{x})=\frac{e^{-\theta \left(\mathbf{x}\right)F_{\mathbf{\eta }}\left(t\right)}-{\pi }_0\left(\mathbf{x}\right)}{1-{\pi }_0\left(\mathbf{x}\right)}$$

(5)

represents the survival function of those who are susceptible to write-off. Note that $S_{\mathbf{\eta }}^s\left(t\right|\mathbf{x}$) is a proper survival function in that $\underset{t\to \infty }{lim}{S}_{\mathbf{\eta }}^s\left(t\right|\mathbf{x})=0$.

The distribution function $F_{\mathbf{\eta }}$ can be specified parametrically, e.g., by considering an exponential, Weibull or gamma distribution or, alternatively, left completely unspecified, resulting in a semiparametric PTC model. For an in-depth discussion relating to the estimation and inference of (4), the reader is referred to the review paper by Amico and van Keilegom [23].

We argue that the PTC model is also interpretable in the post-default loan resolution context. There exist several obligor- and loan-specific factors that may influence a loan’s probability of being written off. Some of these factors are observable and can be included in a model as covariates, but some factors, for example, an individual’s level of discretionary expenditure and undisclosed debt, are latent and therefore analogous to the carcinogenic cells modelled by Chen et al. [25]. There also exists a proportion of defaulted loans for which the outstanding amounts will be fully recovered and is therefore not exposed to write-off. In our context, these loans can be considered “cured”.

de Oliveira and Louzada [5] made use of a PTC model to model the recovery process of defaulted loans within the context of latent competing risks. The considered event of interest is full loan recovery, and the cure proportion represents the proportion of obligors not susceptible to fully repay the outstanding loan amount.

The time to write-off or full recovery may, in addition to latent competing risks, also be influenced by common, unobservable drivers, such as the state of the economy. Within the credit risk environment and specifically in modelling the time to default, several authors have accounted for these unobservable factors by including a frailty parameter within the CPH framework (see e.g., [29,30,31]). An advantageous feature of the CPH model is the ease with which a frailty parameter can be included by multiplication with the baseline hazard function. A frailty parameter is a random variable that can account for unobserved heterogeneity in a population [32]. With the univariate frailty model, the frailty term is assumed to be a random variable that affects the failure times of individuals independently. In a shared frailty model, on the other hand, the frailty term is assumed to be shared among the individuals within a cluster or group [33]. For the CPH model, the frailty parameter is incorporated as follows:

$$h\left(t\right|\mathbf{x},W)=W{h}_0\left(t\right)e^{{\mathbf{\beta }}^T\mathbf{x}},$$

(6)

where the frailty term W is a random variable with non-negative support, such as the gamma distribution, inverse Gaussian (IG) distribution or log-normal distribution. The conditional survival function of the CPH model, given the frailty $W=w$, is

$$S\left(t\right|\mathbf{x},w)=P\left(T>t|\mathbf{x},W=w\right)=S_0{\left(t\right)}^{w\exp \left({\mathbf{\beta }}^T\mathbf{x}\right)}.$$

In the context of the post-default loan resolution process, Betz et al. [3] accounted for unobservable systematic effects in the comovement of the resolution times of defaulted loans by including a random frailty parameter in the CPH model.

We propose including a frailty component in the PTC model of Chen et al. [25] to model the time to write-off. The frailty parameter, which is assumed to have a parametric form, can control for common unobservable factors, as well as for possible observable covariates not included in the model. The PTC model with an IG frailty was used by Barriga et al. [13] in a survival study of colorectal cancer patients, which is, to the best of our knowledge, the only study to include a frailty effect within the PTC framework.

3. Proposed Model: The Promotion Time Cure Model with Parametric Frailty

To account for the heterogeneity in the time to write-off of individual obligors, we propose an extension of the model in (4) by including a parametric frailty component. We opt for a parametric model due to its interpretability and analytical tractability. These characteristics are highly valued in the bank risk management environment, where most models are subjected to intense regulatory scrutiny. A further benefit of parametric models is the ease with which they can be subjected to stress and sensitivity testing.

We first consider the univariate case and include a gamma distributed frailty parameter as follows. Assume that N has a Poisson distribution with parameter $W\theta \left(\mathbf{x}\right),$ where W is an unobserved gamma distributed random variable with parameter ${\sigma }^2>0$ and density function

$$f_{{\sigma }^2}\left(w\right)=\frac{{\sigma }^{-2/{\sigma }^2}}{\Gamma (1/{\sigma }^2)}{w}^{1/{\sigma }^2-1}{e}^{-w/{\sigma }^2},w\ge 0.$$

(7)

Li et al. [34] showed that the model is identifiable if the distribution function of W is specified. In this case, $\mathrm{E}\left[W\right]=1$, which assures that the model is identifiable, and $\mathrm{Var}\left[W\right]={\sigma }^2$.

From (4), we derive the conditional survival function, given $W=w$, as

$${\tilde{S}}_{\mathbf{\eta }}\left(t\right|\mathbf{x},w)=P\left(T>t|\mathbf{x},W=w\right)=e^{-w\theta \left(\mathbf{x}\right)F_{\mathbf{\eta }}\left(t\right)}.$$

(8)

The unconditional survival function is then

$$\begin{array}{cc}\hfill {\tilde{S}}_{\mathbf{\eta }}\left(t\right|\mathbf{x})& ={\int }_0^{\infty }{\tilde{S}}_{\mathbf{\eta }}\left(t\right|\mathbf{x},w)f_W\left(w\right)\mathrm{d}w\hfill \\ \hfill & ={\int }_0^{\infty }{e}^{-w\theta \left(\mathbf{x}\right)F_{\mathbf{\eta }}\left(t\right)}{f}_W\left(w\right)\mathrm{d}w\hfill \\ \hfill & ={\phi }_w\left(\theta \left(\mathbf{x}\right)F_{\mathbf{\eta }}\left(t\right)\right)\hfill \\ \hfill & ={\left\{1+{\sigma }^2\theta \left(\mathbf{x}\right)F_{\mathbf{\eta }}\left(t\right)\right\}}^{-1/{\sigma }^2},\hfill \end{array}$$

(9)

where ${\phi }_w\left(s\right)=\mathrm{E}\left[e^{-sW}\right]$ is the Laplace transform of the frailty W.

It is clear from (9) that

$$\underset{t\to \infty }{lim}{\tilde{S}}_{\mathbf{\eta }}\left(t\right|\mathbf{x})={\left\{1+{\sigma }^2\theta \left(\mathbf{x}\right)\right\}}^{-1/{\sigma }^2}={\pi }_0\left(\mathbf{x}\right),$$

(10)

where ${\pi }_0$ is the cure fraction of the PTC model with gamma frailty, and ${\sigma }^2$ can be seen as a heterogeneity parameter.

From some algebraic manipulation, it follows that (9) can be expressed as

$${\tilde{S}}_{\mathbf{\eta }}\left(t\right|\mathbf{x})={\pi }_0\left(\mathbf{x}\right)+\left(1-{\pi }_0\left(\mathbf{x}\right)\right){\tilde{S}}_{\mathbf{\eta }}^s\left(t\right|\mathbf{x}),$$

(11)

where

$${\tilde{S}}_{\mathbf{\eta }}^s\left(t\right|\mathbf{x})=\frac{{\left\{1+{\sigma }^2\theta \left(\mathbf{x}\right)F_{\mathbf{\eta }}\left(t\right)\right\}}^{-1/{\sigma }^2}-{\pi }_0\left(\mathbf{x}\right)}{1-{\pi }_0\left(\mathbf{x}\right)}$$

(12)

represents the survival function of those who are susceptible to write-off. Note that ${\tilde{S}}_{\mathbf{\eta }}^s\left(t\right|\mathbf{x})$ is a proper survival function in the sense that $\underset{t\to \infty }{lim}{\tilde{S}}_{\mathbf{\eta }}^s\left(t\right|\mathbf{x})=0$ and $\underset{t\to 0}{lim}{\tilde{S}}_{\mathbf{\eta }}^s\left(t\right|\mathbf{x})=1$. From (11), it is clear that the survival function of the non-susceptible individuals is 1, whereas the proper survival function of those who are susceptible is a function of the frailty parameter, ${\sigma }^2$. However, from (10), it follows that the cure fraction is dependent on the frailty parameter.

4. Estimation and Simulation

In this section, we discuss estimation of the parametric PTC model with a gamma frailty parameter using maximum likelihood estimation (MLE). We illustrate the small sample performance of the MLEs through a simulation study.

4.1. Maximum Likelihood Estimation of the PTC Model with Gamma Frailty

For the jth individual, $j=1,2,\dots ,n$, we observe the triplet $(Y_j,{\delta }_j,{\mathbf{x}}_j)$, where $Y_j=min(T_j,C_j)$ is the jth observed time, ${\delta }_j=I\left(T_j\le C_j\right)$ is the corresponding indicator, $T_j$ is the time to write-off, $C_j$ the censoring time and ${\mathbf{x}}_j={\left(1,x_{1j},x_{2j},\dots ,x_{mj}\right)}^T$ is the vector of m covariates associated with the jth individual. Throughout, we assume that the model specified in (8) is a parametric model, i.e., we assume that F has a parametric form with an unknown parameter $\mathbf{\eta }$ and that $\theta$ is related to the covariates by the relationship ${\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_j\right)=e^{{\mathbf{\beta }}^T{\mathbf{x}}_j},$ where $\mathbf{\beta }={\left({\beta }_0,{\beta }_1,\dots ,{\beta }_m\right)}^T$ is a vector of unknown regression parameters.

The likelihood function is given by

$$\begin{array}{cc}\hfill L& =\prod _{j=1}^n{\tilde{S}}_{\mathbf{\eta }}\left(Y_j\right|{\mathbf{x}}_j)h{\left(Y_j\right|{\mathbf{x}}_j)}^{{\delta }_j}\hfill \\ \hfill & =\prod _{j=1}^n{\tilde{S}}_{\mathbf{\eta }}\left(Y_j\right|{\mathbf{x}}_j){\left(-\frac{\mathrm{d}}{\mathrm{d}t}log{\tilde{S}}_{\mathbf{\eta }}\left(Y_j\right|{\mathbf{x}}_j\right)}^{{\delta }_j},\hfill \end{array}$$

where ${\tilde{S}}_{\mathbf{\eta }}\left(t\right|{\mathbf{x}}_j)$ is given in (9). The log-likelihood of the PTC model with gamma frailty is therefore

$$\ell \left(\mathbf{\beta },\mathbf{\eta },{\sigma }^2\right)=-\sum _{j=1}^n\frac{1}{{\sigma }^2}log\left(1+{\sigma }^2{\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_j\right)F_{\mathbf{\eta }}\left(Y_j\right)\right)+\sum _{j=1}^n{\delta }_{j}log\left(\frac{{\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_j\right)f_{\mathbf{\eta }}\left(Y_j\right)}{1+{\sigma }^2{\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_j\right)F_{\mathbf{\eta }}\left(Y_j\right)}\right).$$

(13)

For mixture cure models, a popular approach to estimate model parameters is to apply the EM algorithm. To this end, the mixture cure representation of the model in (11) is used. Recall that in a mixture cure model,

$$S^s\left(t\right|\mathbf{x})=P\left(T>t|B=1,\mathbf{x}\right),$$

where B is the cure status (or susceptibility) indicator, $B=1$ indicates that an individual is susceptible to the event of interest (i.e., uncured) and $B=0$ corresponds to individuals who are non-susceptible (i.e., cured). For the proposed model, (10) indicates that

$$P\left(B=0|\mathbf{x}\right)={\pi }_0\left(\mathbf{x}\right)={\left\{1+{\sigma }^2\theta \left(\mathbf{x}\right)\right\}}^{-1/{\sigma }^2}.$$

Suppose that for the jth individual, $j=1,2,\dots ,n$, we observe not only the triplet $(Y_j,{\delta }_j,{\mathbf{x}}_j)$ but also $B_j$, the indicator showing whether an observation is susceptible or not. It is then possible to define the complete likelihood (see [35]) as

$$\begin{array}{cc}\hfill L& =\prod _{j=1}^n{\left[\left(1-{\pi }_0\left({\mathbf{x}}_j\right)\right){\tilde{h}}_{\mathbf{\eta }}^s\left(Y_j\right|\mathbf{x}){\tilde{S}}_{\mathbf{\eta }}^s\left(Y_j\right|\mathbf{x})\right]}^{{\delta }_j{B}_j}\hfill \\ \hfill & \times \prod _{j=1}^n{\left[{\pi }_0{\left({\mathbf{x}}_j\right)}^{1-B_j}{\left\{\left(1-{\pi }_0\left({\mathbf{x}}_j\right)\right){\tilde{S}}_{\mathbf{\eta }}^s\left(Y_j\right|{\mathbf{x}}_j)\right\}}^{B_j}\right]}^{1-{\delta }_j},\hfill \end{array}$$

where ${\pi }_0\left({\mathbf{x}}_j\right)$ and ${\tilde{S}}_{\mathbf{\eta }}^s$ are given in (10) and (12), respectively, while

$${\tilde{h}}_{\mathbf{\eta }}^s\left(t\right|{\mathbf{x}}_j)=-\frac{\mathrm{d}}{\mathrm{d}t}log{\tilde{S}}_{\mathbf{\eta }}^s\left(t\right|{\mathbf{x}}_j)$$

represents the hazard function of those who are susceptible to write-off.

This likelihood can be rewritten as

$$\begin{array}{cc}\hfill L& =\prod _{j=1}^n{\left(1-{\pi }_0\left({\mathbf{x}}_j\right)\right)}^{{\delta }_j{B}_j}{\pi }_0{\left({\mathbf{x}}_j\right)}^{(1-{\delta }_j)(1-B_j)}{\left(1-{\pi }_0\left({\mathbf{x}}_j\right)\right)}^{(1-{\delta }_j)B_j}\hfill \\ \hfill & \times \prod _{j=1}^n{\left[{\tilde{h}}_{\mathbf{\eta }}^s\left(Y_j\right|\mathbf{x}){\tilde{S}}_{\mathbf{\eta }}^s\left(Y_j\right|\mathbf{x})\right]}^{{\delta }_j{B}_j}{\left\{{\tilde{S}}_{\mathbf{\eta }}^s\left(Y_j\right|{\mathbf{x}}_j)\right\}}^{(1-{\delta }_j)B_j}.\hfill \end{array}$$

For the mixture cure model in (3), the survival function of the susceptible individuals $S^s\left(t\right|\mathbf{y})$ is not a direct function of the cure fraction ${\pi }_0$. Unfortunately, this is not the case for the survival function ${\tilde{S}}^s\left(t\right|\mathbf{x}),$ when the PTC model is written in its mixture cure form, as can be seen from (12). Consequently, unlike for the “traditional” mixture cure model, the likelihood presented above cannot be split into two separate components, namely a latency part and an incidence part, which could have then be maximised separately (see, e.g., [9,36]). We therefore settle for the conventional method to obtain maximum likelihood estimators for the parameters of the PTC model with gamma frailty, i.e., we find the set of parameter values ($\beta$, $\mathbf{\nu }$ and ${\sigma }^2$) that maximises the log-likelihood function in (13). The required optimisation was performed by applying the Nelder–Mead algorithm, which is included in the optim function of the R statistical software package [37].

Similar to the study by Barriga et al. [13], we also consider a PTC model with IG frailty. In this case, the log-likelihood function is given by

$$\begin{array}{cc}\hfill \ell \left(\mathbf{\beta },\lambda ,k,\zeta \right)& =\frac{n}{\zeta }+\sum _{j=1}^n{\delta }_{j}log{\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_j\right)+\sum _{j=1}^n{\delta }_{j}log{F}_{\mathbf{\eta }}\left(Y_j\right)\hfill \\ \hfill & -\frac{1}{2}\sum _{j=1}^n{\delta }_{j}log\left(1+2\zeta {\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_j\right)F_{\mathbf{\eta }}\left(Y_j\right)\right)-\frac{1}{\zeta }\sum _{j=1}^n\sqrt{1+2\zeta {\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_j\right)F_{\mathbf{\eta }}\left(Y_j\right)},\hfill \end{array}$$

(14)

where $\zeta$ is the frailty parameter.

In the simulation study, we choose $F_{\mathbf{\eta }}$ as the Weibull and linear failure rate (LFR) distributions. In the case of the Weibull distribution, $\mathbf{\eta }={(\lambda ,k)}^T$, where $\lambda >0$ is a scale parameter and $k>0$ is a shape parameter with distribution function

$$F_{\lambda ,k}\left(t\right)=1-e^{-{\left(\frac{t}{\lambda }\right)}^k},t\ge 0$$

(15)

and density function

$$f_{\lambda ,k}\left(t\right)=\frac{k}{\lambda }{\left(\frac{t}{\lambda }\right)}^{k-1}{e}^{-{\left(\frac{t}{\lambda }\right)}^k},t\ge 0.$$

(16)

In the case of the LFR distribution, $\mathbf{\eta }={(a,b)}^T$, where $a>0$ and $b>0$ are shape parameters, and the distribution function is given by

$$F_{a,b}\left(t\right)=1-e^{-at-\frac{b}{2}{t}^2},t\ge 0$$

with density function

$$f_{a,b}\left(t\right)=\left(a+bt\right)e^{-at-\frac{b}{2}{t}^2},t\ge 0.$$

4.2. Simulation Study

To evaluate the performance of the MLEs, based on the log-likelihood given in (13) (and also for comparison purposes, the log-likelihood given in (14)), we conducted a simulation study as follows. Two covariates are used; the first, ${\mathbf{x}}_{j1},j=1,\dots ,n$, is generated from a standard exponential distribution, whereas the second, ${\mathbf{x}}_{j2},j=1,\dots ,n$, is generated from a Bernoulli distribution with a probability of success of $0.6$. These covariates remain fixed throughout the Monte Carlo simulation study. The cure proportion is given by

$$p_j={\left\{1+{\sigma }^2{\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_j\right)\right\}}^{-1/{\sigma }^2},$$

where

$${\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_j\right)=e^{{\mathbf{\beta }}^T{\mathbf{x}}_j}=e^{{\beta }_0+{\beta }_1{x}_{j1}+{\beta }_2{x}_{j2}}.$$

As mentioned earlier, we take $F_{\mathbf{\eta }}$ to be either the Weibull or LFR distribution. In the former case, $\mathbf{\eta }={(\lambda ,k)}^T={(2,1)}^T$, and in the latter case, $\mathbf{\eta }={(a,b)}^T={(1.5,2)}^T$. Censoring times $C_j,j=1,\dots ,n$ are sampled from a standard uniform distribution, and the frailty parameter, ${\sigma }^2$, is set to $0.5$. In order to control the cured and censoring proportions, we consider two different values for the regression coefficient vector, namely $\mathbf{\beta }={(2,-1,2)}^T$ and $\mathbf{\beta }={(0.5,-1,2)}^T$. The first calibration (in conjunction with the parameter choices for $F_{\mathbf{\eta }}$) results in an average cure proportion of $10\%$ and a $20\%$ censoring proportion, while the second results in an average cure proportion of $30\%$, with a $45\%$ censoring proportion. The jth resolution time $T_j$ from a PTC model with gamma frailty is generated as follows:

• Generate $B_j$ from a Bernoulli distribution:

  $$B_j=\left\{\begin{array}{cc}1\hfill & \mathrm{with}\mathrm{probability}1-p_j\hfill \\ 0\hfill & \mathrm{with}\mathrm{probability}{p}_j\hfill \end{array}\right.$$
• If $B_j=0$, then set $T_j=\infty$;
• If $B_j=1$, generate $U_j\sim U[0,1]$, and calculate

  $$T_j=F_{\mathbf{\eta }}^{-1}\left\{\frac{{\left(U_j(1-p_j)+p_j\right)}^{-{\sigma }^2}-1}{{\sigma }^2{\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_j\right)}\right\},$$

  where $F_{\mathbf{\eta }}^{-1}$ is the inverse cumulative distribution function;
• Generate censoring times $C_j\sim U[0,1]$;
• The simulated data are then $(Y_j,{\delta }_j,{\mathbf{x}}_{j1},{\mathbf{x}}_{j2}),j=1,\dots ,n$, where

  $$Y_j=min(T_j,C_j),$$

  and

  $${\delta }_j=I(T_j\le C_j).$$

All calculated values are based on 1000 independent Monte Carlo replications for samples of size $n=200,300,500$ and 1000.

Table 1. Simulation results with $10\%$ cure fraction and $20\%$ censoring.

Parameter	$\mathit{n}=200$	$\mathit{n}=300$	$\mathit{n}=500$	$\mathit{n}=1000$
PTC with gamma frailty and $F_{\mathbf{\eta }}$ a Weibull distribution
${\beta }_0=2$	2.314 (0.988)	2.165 (0.706)	2.120 (0.512)	2.049 (0.330)
${\beta }_1=-1$	−1.035 (0.191)	−1.016 (0.154)	−1.016 (0.116)	−1.004 (0.085)
${\beta }_2=2$	2.047 (0.347)	2.025 (0.276)	2.017 (0.216)	2.011 (0.152)
$\lambda =2$	2.070 (1.217)	2.064 (0.988)	2.027 (0.782)	1.998 (0.560)
$\nu =1$	1.018 (0.113)	1.012 (0.090)	1.009 (0.069)	1.002 (0.048)
${\sigma }^2=0.5$	0.517 (0.302)	0.506 (0.247)	0.516 (0.179)	0.502 (0.127)
PTC with gamma frailty and $F_{\mathbf{\eta }}$ a LFR distribution
${\beta }_0=2$	2.318 (0.833)	2.199 (0.651)	2.170 (0.493)	2.110 (0.353)
${\beta }_1=-1$	−1.059 (0.195)	−1.028 (0.152)	−1.022 (0.108)	−1.009 (0.076)
${\beta }_2=2$	2.103 (0.361)	2.050 (0.262)	2.029 (0.197)	2.023 (0.135)
$a=1.5$	1.338 (0.804)	1.391 (0.682)	1.381 (0.551)	1.398 (0.415)
$b=2$	3.337 (4.018)	2.820 (2.966)	2.531 (2.510)	2.257 (1.882)
${\sigma }^2=0.5$	0.569 (0.308)	0.529 (0.238)	0.528 (0.161)	0.513 (0.107)
PTC with IG frailty and $F_{\mathbf{\eta }}$ a Weibull distribution
${\beta }_0=2$	2.556 (1.461)	2.385 (1.227)	2.312 (1.018)	2.107 (0.496)
${\beta }_1=-1$	−1.079 (0.234)	−1.057 (0.210)	−1.050 (0.167)	−1.018 (0.122)
${\beta }_2=2$	2.134 (0.447)	2.110 (0.407)	2.089 (0.330)	2.040 (0.231)
$\lambda =2$	2.228 (1.346)	2.181 (1.044)	2.096 (0.806)	2.046 (0.558)
$\nu =1$	1.060 (0.168)	1.051 (0.159)	1.042 (0.130)	1.016 (0.090)
$\zeta =0.5$	1.638 (3.506)	1.554 (3.553)	1.213 (2.284)	0.768 (1.261)
PTC with IG frailty and $F_{\mathbf{\eta }}$ a LFR distribution
${\beta }_0=2$	2.344 (0.973)	2.226 (0.703)	2.183 (0.497)	2.106 (0.332)
${\beta }_1=-1$	−1.058 (0.184)	−1.035 (0.150)	−1.023 (0.105)	−1.012 (0.077)
${\beta }_2=2$	2.095 (0.337)	2.069 (0.278)	2.037 (0.196)	2.028 (0.139)
$a=1.5$	1.362 (0.816)	1.382 (0.673)	1.374 (0.551)	1.404 (0.402)
$b=2$	3.459 (4.147)	3.052 (3.301)	2.585 (2.561)	2.350 (2.054)
$\zeta =0.5$	0.927 (1.549)	0.886 (2.742)	0.637 (0.503)	0.567 (0.293)

contains the approximate average values of the MLEs, as well as the corresponding standard errors given in parenthesis, for the PTC model with either the gamma or IG frailty and $F_{\mathbf{\eta }}$ as either a Weibull or LFR distribution. The average cure fraction is $10\%$, and the censoring proportion is $20\%$.

Table 2. Simulation results with $30\%$ cure fraction and $45\%$ censoring.

Parameter	$\mathit{n}=200$	$\mathit{n}=300$	$\mathit{n}=500$	$\mathit{n}=1000$
PTC with gamma frailty and $F_{\mathbf{\eta }}$ a Weibull distribution
${\beta }_0=0.5$	0.789 (0.887)	0.686 (0.755)	0.581 (0.409)	0.546 (0.254)
${\beta }_1=-1$	−1.067 (0.245)	−1.023 (0.189)	−1.022 (0.141)	−1.010 (0.102)
${\beta }_2=2$	2.097 (0.426)	2.038 (0.337)	2.020 (0.247)	2.021 (0.179)
$\lambda =2$	2.153 (1.188)	2.078 (0.982)	2.077 (0.751)	1.999 (0.533)
$\nu =1$	1.039 (0.135)	1.020 (0.109)	1.013 (0.079)	1.005 (0.055)
${\sigma }^2=0.5$	0.595 (0.480)	0.530 (0.375)	0.525 (0.273)	0.517 (0.182)
PTC with gamma frailty and $F_{\mathbf{\eta }}$ a LFR distribution
${\beta }_0=0.5$	0.812 (0.874)	0.695 (0.591)	0.618 (0.395)	0.579 (0.262)
${\beta }_1=-1$	−1.092 (0.264)	−1.039 (0.200)	−1.030 (0.142)	−1.014 (0.098)
${\beta }_2=2$	2.158 (0.494)	2.073 (0.361)	2.038 (0.252)	2.029 (0.175)
$a=1.5$	1.328 (0.830)	1.344 (0.669)	1.399 (0.518)	1.409 (0.370)
$b=2$	3.045 (3.104)	2.705 (2.683)	2.476 (2.226)	2.145 (1.645)
${\sigma }^2=0.5$	0.671 (0.578)	0.572 (0.401)	0.548 (0.280)	0.526 (0.174)
PTC with IG frailty and $F_{\mathbf{\eta }}$ a Weibull distribution
${\beta }_0=0.5$	0.958 (1.203)	0.802 (0.950)	0.746 (0.791)	0.623 (0.433)
${\beta }_1=-1$	−1.095 (0.271)	−1.062 (0.241)	−1.054 (0.205)	−1.033 (0.146)
${\beta }_2=2$	2.164 (0.51)	2.135 (0.466)	2.108 (0.383)	2.069 (0.286)
$\lambda =2$	2.202 (1.254)	2.206 (1.014)	2.131 (0.771)	2.018 (0.521)
$\nu =1$	1.070 (0.179)	1.060 (0.166)	1.050 (0.137)	1.023 (0.100)
$\zeta =0.5$	2.139 (5.047)	1.988 (4.589)	1.544 (3.243)	1.045 (2.389)
PTC with IG frailty and $F_{\mathbf{\eta }}$ a LFR distribution
${\beta }_0=0.5$	0.932 (0.970)	0.829 (0.970)	0.690 (0.488)	0.628 (0.341)
${\beta }_1=-1$	−1.098 (0.266)	−1.068 (0.235)	−1.038 (0.170)	−1.028 (0.120)
${\beta }_2=2$	2.168 (0.483)	2.142 (0.426)	2.080 (0.318)	2.058 (0.238)
$a=1.5$	1.256 (0.795)	1.298 (0.684)	1.350 (0.558)	1.374 (0.412)
$b=2$	3.138 (3.305)	2.888 (2.833)	2.584 (2.381)	2.248 (1.888)
$\zeta =0.5$	5.661 (70.567)	2.762 (20.101)	1.108 (2.504)	0.781 (1.191)

shows the results for a $30\%$ cure fraction and a $45\%$ censoring proportion, with all other aspects identical to the set-up of Table 1.

Table 1 and Table 2 show the properties that we expect from MLEs, namely that they are unbiased and consistent. This is shown by the average estimates approaching the true values of the parameters, as well as the decreasing standard errors, as the sample size increases. Table 1 shows that for the IG frailty parameter $\zeta$, the bias and the standard errors of the estimates are greater for smaller samples when $F_{\mathbf{\eta }}$ is the Weibull distribution compared to the LFR distribution. We also observe a slower rate of convergence in this case. The same holds for the rate of convergence in the case of a higher cure fraction and censoring percentage, as observed in Table 2. When $F_{\mathbf{\eta }}$ is the LFR distribution, larger standard errors are observed for smaller sample sizes ($n=200$ and $n=300$). Overall, the results demonstrate that the MLEs, based on the log-likelihood for the PTC model with gamma frailty, are accurate for small samples, especially when comparing their performance to those of the PTC model with an IG frailty. Although not shown in this paper, confidence intervals for the parameters can be estimated using the following bootstrap algorithm, which makes use of the ”cases bootstrap” approach discussed in [38]:

• Resample $(Y_1^{*},{\delta }_1^{*},{\mathbf{x}}_1^{*}),\dots ,(Y_n^{*},{\delta }_n^{*},{\mathbf{x}}_n^{*})$ i.i.d from the empirical distribution function of $(Y_1,{\delta }_1,{\mathbf{x}}_1),\dots ,(Y_n,{\delta }_n,{\mathbf{x}}_n)$;
• Fit the model in (9) on the resampled data to obtain estimates ${\widehat{\mathbf{\vartheta }}}^{*}={\left({\widehat{\mathbf{\beta }}}^{*},{\widehat{\mathbf{\eta }}}^{*},{\widehat{\sigma }}^{*}\right)}^T$. Denote the first of these by ${\widehat{\mathbf{\vartheta }}}^{*}\left(1\right)$;
• Repeat steps 1 and 2 B times to obtain ${\widehat{\mathbf{\vartheta }}}^{*}\left(1\right),{\widehat{\mathbf{\vartheta }}}^{*}\left(2\right),\dots ,{\widehat{\mathbf{\vartheta }}}^{*}\left(B\right)$.

Now, suppose that we want to construct a confidence interval for the parameter ${\beta }_1$. From ${\widehat{\mathbf{\vartheta }}}^{*}\left(1\right),{\widehat{\mathbf{\vartheta }}}^{*}\left(2\right),\dots ,{\widehat{\mathbf{\vartheta }}}^{*}\left(B\right)$, we have the bootstrap estimators ${\widehat{\beta }}_1^{*}\left(1\right),\dots ,{\widehat{\beta }}_1^{*}\left(B\right)$, from which a confidence interval can be obtained. Alternatively, one can make use of the large sample properties of MLEs to construct a confidence interval (see, e.g., [13]).

5. Application

In this section, we analyse whether the proposed models, as described in Section 3, are able to account for common unobservable factors by fitting them to two datasets. Before fitting the frailty models to a real-world dataset, we apply the models to a simulated dataset, where we simulate an “unobservable” variable not directly accounted for in our frailty models. We then proceed to apply the proposed models to a real-world loss dataset obtained from Global Credit Data (GCD).

5.1. Simulated Dataset

In order to further investigate the performance of the PTC model with gamma frailty, we simulate a single dataset $(Y_j,{\delta }_j,x_{j1},x_{j2}),j=1,\dots ,n,$ of size $n=1000$ from a PTC model without frailty, i.e., from the model in (4). The two covariates are generated as in Section 4.2, whereas the cure proportion is now given by $p_j=e^{-{\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_j\right)}$, where

$${\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_j\right)=e^{{\mathbf{\beta }}^T{\mathbf{x}}_j}=e^{{\beta }_0+{\beta }_1{x}_{j1}+{\beta }_2{x}_{j2}}$$

with $\mathbf{\beta }={(2,-1,2)}^T$. $F_{\mathbf{\eta }}$ is taken to be the Weibull distribution, with the same parameters as in Section 4.2 and censoring times are again assumed to follow a standard uniform distribution. With this configuration, our simulated dataset has a non-susceptible percentage of $6.4\%$, and $16.2\%$ of the resolution times are censored.

One of the covariates is now dropped from the simulated dataset so that we only “observe” the set $(Y_j,{\delta }_j,x_{j1}),j=1,\dots ,n$. On this observed dataset, we fit three different models:

• The PTC model with no frailty in (4);
• The PTC model in (8) with gamma frailty and frailty parameter ${\sigma }^2$;
• The PTC model in (8) with an IG frailty and frailty parameter $\zeta$.

For all three models, we have

$${\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_j\right)=e^{{\mathbf{\beta }}^T{\mathbf{x}}_j}=e^{{\beta }_0+{\beta }_1{x}_{j1}}.$$

In

Table 3. Parameter estimates and goodness-of-fit measures for all three models fitted to the simulated dataset.

Parameter	PTC	PTC with Gamma Frailty	PTC with IG Frailty
${\beta }_0$	$1.724$	$3.695$	$3.603$
${\beta }_1$	$-0.629$	$-0.954$	$-0.978$
$\lambda$	$3.513$	$1.589$	$2.778$
k	$0.790$	$0.967$	$1.025$
${\sigma }^2$		$0.825$ *
$\zeta$			$2.244$ *
log-likelihood	$1485.493$	$1493.001$	$1496.808$
AIC	$-2962.986$	$-2976.002$	$-2983.616$
BIC	$-2943.355$	$-2951.463$	$-2959.077$

* These parameters are significant at a 5% significance level according to the adjusted likelihood ratio test (see Claeskens et al. [39]).

, the estimates of all three models are shown, as well as their respective calculated log likelihoods, Akaike information criteria (AIC) and Bayesian information criteria (BIC).

Based on both the AIC and BIC, the two frailty models outperform the PTC model, with the IG frailty yielding the best performance. The estimates of the ${\beta }_0$ and ${\beta }_1$ parameters for the two frailty models are very similar but differ notably from the estimate for the PTC model. This highlights the idiosyncratic attributes of the frailty models compared to the PTC model and, in conjunction with the lower AIC and BIC values, suggests that the frailty models are better at capturing the effect of the unobserved covariate than the PTC model.

5.2. Real-World Dataset

Our real-world loss dataset was provided by GCD, formally known as the Pan-European Credit Data Consortium (PECDC). GCD is a non-profit organisation that collects and analyses historical loss data from member banks [40].

We use a subsample of the loss database consisting of large corporations (LCs) from the United States. All loans that defaulted in the period 1 January 2005 to 31 December 2019 are included. For materiality purposes, we further filter the data to only include loans for which the exposure amount is larger than USD 500. The dataset includes information on when the loan went into default, whether (and when) the loan was fully settled, whether (and when) the loan was written off and whether the case has been resolved.

We include the industry of the obligor as a categorical covariate in the model based on the study by Schuermann [41], which demonstrates that the industry of an obligor has an impact on the recovery rate. In industries in which tangible assets represent a substantial portion of a company’s assets, such as manufacturing or transportation, the bankruptcy process may be resolved more quickly, since the assets can be more easily sold to repay creditors.

Obligor industries are aggregated into five industry groups, as indicated in

Table 4. Summary of industry groupings, industry frequency and censored proportion of US large corporation dataset.

Industry Group		Industry Description	Frequency	Censored
1	Primary Industries	Agriculture, Hunting and Forestry	$14.5\%$	$28.8\%$
		Fishing and Fishing Products
		Mining
2	Manufacturing, Utilities and Construction	Manufacturing	$33.3\%$	$33.5\%$
		Utilities
		Construction
3	Trade Industries	Wholesale and Retail Trade	$24\%$	$28.1\%$
		Transportation and Storage
4	Service Industries	Hotels and Restaurants	$17\%$	$27.9\%$
		Finance and Insurance
		Real Estate Rental and Leasing
		Professional, Scientific and Technical Services
		Private Sector Services (Household)
5	Other	Communications	$4.7\%$	$35.6\%$
		Public Administration and Defence
		Education
		Health and Social Services
		Other Community, Social and Personal Services
−1	Unknown	Unknown	$6.5\%$	$30.7\%$

, with a sixth (denoted by $-1$) for loans for which the industry is unknown. Also shown in the rightmost column of Table 4 is the proportion of censored cases for each industry group. We consider defaulted loans that were fully repaid, as well as those for which the resolution process had not yet been completed at the end of the observation period, as censored.

In Figure 1, the Kaplan–Meier estimate of the survival probability of all data is shown on the left, and that of the data stratified by industry ($-1,1,\dots ,5$) is shown on the right. From the left pane graph, it is evident that a cure proportion is present in the aggregated dataset. The narrow confidence interval can be attributed to the large sample size ($n=1558$). From the right pane graph, notable differences in the cure proportions for the different industry groups can be observed, suggesting that the industry of a borrower has a definite impact on the probability of recovering the full loan amount after default.

Besides industry groups, which are included as dummy variables with corresponding coefficients (${\beta }_2,\dots ,{\beta }_6$), we include the annualised seasonally adjusted GDP growth in the month of default as a covariate, with the associated coefficient ${\beta }_1$. Our choice of GDP growth is motivated by the study of Khieu et al. [42], in which the authors report a positive relationship between annual GDP growth and recovery rates of defaulted corporate debt in the US. Furthermore, Betz et al. [43] found positive dependencies between default resolution times and final loan loss rates.

The MLEs of the model parameters are shown in the first three columns of

Table 5. Estimation results of different models on the US large corporation dataset.

Parameter	PTC	PTC with Gamma Frailty	PTC with IG Frailty	PTC with Shared Frailty
${\beta }_0$	3.920	6.968	2.141	0.843
${\beta }_1$	−0.025	−0.021	−0.029	−0.023
${\beta }_2$	0.187	0.094	0.134
${\beta }_3$	0.037	0.254	0.083
${\beta }_4$	0.175	0.501	0.297
${\beta }_5$	0.050	0.064	0.093
${\beta }_6$	0.076	0.168	0.117
$\lambda$	4.462 $\times {10}^{-4}$	4.644 $\times {10}^{-8}$	2.523 $\times {10}^{-5}$	3.343 $\times {10}^{-4}$
k	0.852	1.579	1.327	1.073
${\sigma }^2$		2.157 *		3.926
$\zeta$			2.360 *
log-likelihood	−8813.626	−8735.423	−8753.194	−8787.344
AIC	17,645.200	17,490.850	17,526.390	17,584.690
BIC	17,693.412	17,544.358	17,579.900	17,611.444

* These parameters are significant at the 5% level according to the adjusted likelihood ratio test.

. For both frailty models, the frailty parameters (${\sigma }^2$ and $\zeta$, respectively) are significant at the 5% level based on the adjusted log-likelihood test. In line with the findings of Khieu et al. [42] and Betz et al. [43], a negative relationship between GDP growth and time to write-off is observed. This means that higher GDP growth at the time of default is associated with a higher probability of full loan repayment. This makes intuitive sense, especially when considering the price at which collateral can be sold or the prospects for a firm to recover from financial distress and settle outstanding loan obligations. Also shown in Table 5 are the calculated log likelihoods, AIC and BIC. Based on all three goodness-of-fit measures, the PTC model with gamma frailty performs the best out of the three models.

In Figure 2, the survival curves (left pane) and hazard rates (right pane) for the PTC model, as well as the PTC with gamma and IG frailties, are plotted. From the left pane graph, it is evident that the fitted PTC model exhibits a much smaller cure fraction than the two frailty models, for which the frailty parameters were found to be significant. This further motivates the inclusion of a frailty parameter within the PTC model to capture residual heterogeneity. We observe that the hazard functions of the frailty models have very different shapes when compared to the standard PTC model. The hazard function has a significant impact on the survival function, especially on the cure rate, as can be seen from the graph on the left-hand side. The hazard functions of the two frailty models both exhibit a non-monotone shape. By including a frailty parameter in the model, the hazard function can become more flexible and may better capture non-standard patterns. Consequently, due to the flexibility of the frailty parameter, possibly less reliance may be placed on the choice of $F_{\mathbf{\eta }}$.

6. Extension to the Shared Gamma Frailty Model

In many survival studies, the failure times of individuals are not independent but dependent within groups or clusters. This can be due to shared environmental factors, such as exposure to air pollution in the case of cancer data or the state of the economy in the case of credit risk. By including a shared frailty term in a survival model, the dependence between the failure times of individual subjects within a cluster can be allowed for, leading to more accurate estimates of the hazard rates and the survival probabilities.

Suppose that we have n clusters of size $n_i,i=1,2,\dots ,n$, with each cluster sharing a common risk, assumed to be random, and modelled by an i.i.d. gamma distributed frailty $W_i,i=1,\dots ,n$, each with a density of $f_{{\sigma }^2}$ given in (7). The multivariate conditional survival function, given $W_i=w_i$, is then

$$\begin{array}{cc}\hfill {\tilde{S}}_{\mathbf{\eta }}(t_{i1},\dots ,t_{i{n}_i}|{\mathbf{x}}_i,w_i)& =P\left(T_{i1}>t_{i1},\dots ,T_{i{n}_i}>t_{i{n}_i}|{\mathbf{x}}_i,W_i=w_i\right)\hfill \\ \hfill & =P\left(T_{i1}>t_{i1}|{\mathbf{x}}_i,W_i=w_i\right)\dots P\left(T_{i{n}_i}>t_{i{n}_i}|{\mathbf{x}}_i,W_i=w_i\right)\hfill \\ \hfill & =e^{-w_i{\sum }_{j=1}^{n_i}{\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_{ij}\right)F_{\mathbf{\eta }}\left(t_{ij}\right)},\hfill \end{array}$$

where ${\mathbf{x}}_i=\left({\mathbf{x}}_{i1},\dots ,{\mathbf{x}}_{i{n}_i}\right)$ is the covariate matrix of the individual in the ith cluster, and $F_{\mathbf{\eta }}$ is, as in the univariate case, some specified distribution with parameter vector $\mathbf{\eta }$.

For the jth individual in the ith cluster, $j=1,2,\dots ,n_i$, we observe the triplet $(Y_{ij},{\delta }_{ij},{\mathbf{x}}_{ij})$, where $Y_{ij}=min(T_{ij},C_{ij})$, ${\delta }_{ij}=I\left(T_{ij}\le C_{ij}\right)$, and ${\mathbf{x}}_{ij}={(1,x_{1ij},x_{2ij},\dots ,x_{mij})}^T$ is the vector of m covariates associated with the jth individual in the ith cluster. Similar to the univariate case, ${\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_{ij}\right)=e^{{\mathbf{\beta }}^T{\mathbf{x}}_{ij}},$ where $\mathbf{\beta }={({\beta }_0,{\beta }_1,\dots ,{\beta }_m)}^T$ is a vector of unknown regression parameters. Furthermore, we denote the number of observed events in the ith cluster by $d_i={\sum }_{j=1}^{n_i}{\delta }_{ij}$. The likelihood is then given by

$$L=\prod _{i=1}^n{\int }_0^{\infty }\left(\prod _{j=1}^{n_i}{\left(w_i{\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_{ij}\right)f_{\mathbf{\eta }}\left(Y_{ij}\right)\right)}^{{\delta }_{ij}}{e}^{-w_i{\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_{ij}\right)F_{\mathbf{\eta }}\left(Y_{ij}\right)}\right)f_{{\sigma }^2}\left(w_i\right)\mathrm{d}{w}_i,$$

where $f_{\mathbf{\eta }}$ is the density associated with distribution $F_{\mathbf{\eta }}$.

By some algebraic manipulation, and by taking the natural logarithm, the log-likelihood function in closed form is given by

$$\begin{array}{cc}\hfill \ell (\mathbf{\beta },\sigma ,\mathbf{\eta })& =\sum _{i=1}^n\left(\sum _{j=1}^{n_i}{\delta }_{ij}log\left\{{\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_{ij}\right)f_{\mathbf{\eta }}\left(Y_{ij}\right)\right\}\right)+\sum _{i=1}^{n}log\Gamma \left(\frac{1}{{\sigma }^2}+d_i\right)\hfill \\ \hfill & -\sum _{i=1}^n\left(\frac{1}{{\sigma }^2}+d_i\right)log\left(\frac{1}{{\sigma }^2}+\sum _{j=1}^{n_i}{\theta }_{\mathbf{\beta }}\left({\mathbf{x}}_{ij}\right)F_{\mathbf{\eta }}\left(Y_{ij}\right))\right)\hfill \\ \hfill & -\frac{2n}{{\sigma }^2}log\sigma -nlog\Gamma \left(\frac{1}{{\sigma }^2}\right).\hfill \end{array}$$

For the application of the shared frailty model to our real-world dataset described in Section 5.2, we cluster our data with reference to industry group based on the intuitive reasoning that obligors within the same industry group are subject to similar forces in the economy, which is a source of dependency within these groups. Instead of treating the industry group as a covariate, a shared frailty parameter is fitted as a means to account for the dependence of write-off times within industry groups. Again, we assume $F_{\mathbf{\eta }}$ to be the Weibull distribution function with $\mathbf{\eta }=(\lambda ,k)$ as defined in (15).

The final column of Table 5 lists the MLEs of the PTC model with shared frailty. It is important to note that the significant frailty variance parameter ${\sigma }^2$ of the univariate gamma frailty model confirms that the model succeeds in capturing unobserved heterogeneity. However, for the shared frailty model, a significant frailty variance parameter ${\sigma }^2$ confirms that a relationship exists within each industry cluster and that its value can be interpreted to measure the dependence between loan write-off times within each industry cluster (see discussion in [32] (p. 144)).

According to the values of the calculated log likelihood, AIC and BIC, the shared frailty model performs favourably compared to the standard PTC model, despite the shared frailty model being a more parsimonious model. This also suggests that the shared frailty component does succeed in capturing some of the dependence between loan write-off times attributable to the industry of the obligor.

7. Concluding Remarks

In this paper, a PTC model with gamma frailty was developed and applied to a real-world loss dataset. The applicability of a cure proportion in the context of the post default loan resolution process was demonstrated, and we also showed how the inclusion of a frailty parameter in the PTC model can allow for the modelling of a more flexible shape of the hazard function. This may improve the accuracy of the survival model and, equally importantly, provide a better understanding of the underlying processes that drive the occurrence of write-off over time. Furthermore, the flexibility offered by the inclusion of a frailty parameter affords for possibly less reliance to be placed on the selection of a “baseline” distribution $F_{\mathbf{\eta }}$ in the PTC model.

The proportional hazard assumption of the popular CPH model implies that the relative hazard rate remains constant over time with different covariate levels. By including a parametric frailty parameter in the PTC model with fewer covariates, a more flexible relationship between potential (latent) predictors and the time to write-off can be accommodated. Through a Monte Carlo study, we demonstrated that the parameter estimates of the PTC model with gamma frailty are reliable.

We also showed that, by including a shared gamma frailty term in a PTC model, the dependence between the failure times of individual subjects within a cluster can be accounted for, which, in turn, may lead to more accurate estimates of the hazard rate and the write-off times. Our initial exploration of the use of a shared frailty model warrants future investigation, with specific focus on the small sample performance of shared frailty models. Other PTC models with shared frailties can also be considered, such as a log-normal shared frailty model or a compound Poisson shared frailty model that includes two important special cases, namely gamma and IG distributions.

A shared frailty model takes the dependency of the event times within clusters into account, but other forms of dependency may also exist that were not explored in this study. The models we considered assume independence between the censoring time and write-off time of individual loans. A model that allows for possible dependence between these times may be considered in a future study. This is referred to as dependent censoring and is rapidly gaining traction in contemporary research (see, e.g., [44,45,46]).

In a highly regulated industry such as banking, parametric models are often preferred to non-parametric models due to their analytical tractability and interpretability. Parametric models also easily allow for sensitivity analysis and stress testing by varying parameter values. As demonstrated by this study, parametric frailty models offer the potential to characterise unexplained heterogeneity in the write-off times of individual loans. Besides the time to write-off, which can be predicted more accurately, these models also offer the prospect of modelling the probability of write-off and full repayment in the context of censored data. As such, the model offers a viable alternative for modelling the resolution outcome in two-step loss, given default modelling, while making use of the most recent (and therefore most relevant) loan resolution data. The reliance of these models on accurate and trustworthy data requires banks to make substantial investments in systems, expertise and validation efforts. Modellers should take care to avoid spurious accuracy and to understand the economic and operational drivers of credit risk before and after loan default. Inclusion of a frailty parameter to capture residual heterogeneity should therefore not be regarded as a substitute for including known risk drivers as covariates within a prediction model.