Bootstrap Bandwidth Selection and Confidence Regions for Double Smoothed Default Probability Estimation

Abstract

For a fixed time, t, and a horizon time, b, the probability of default (PD) measures the probability that an obligor, that has paid his/her credit until time t, runs into arrears not later that time $t+b$. This probability is one of the most crucial elements that influences the risk in credits. Previous works have proposed nonparametric estimators for the probability of default derived from Beran’s estimator and a doubly smoothed Beran’s estimator of the conditional survival function for censored data. They have also found asymptotic expressions for the bias and variance of the estimators, but they do not provide any practical way to choose the smoothing parameters involved. In this paper, resampling methods based on bootstrap techniques are proposed to approximate the bandwidths on which Beran and smoothed Beran’s estimators of the PD depend. Bootstrap algorithms for the calculation of confidence regions of the probability of default are also proposed. Extensive simulation studies show the good behavior of the presented algorithms. The bandwidth selector and the confidence region algorithm are applied to a German credit dataset to analyze the probability of default conditional on the credit scoring.

1. Introduction

The debts coming from clients with unpaid credits have an important impact in the solvency of banks and other credit institutions. According to the Basel Committee on Banking Supervision of the Bank for International Settlements, one of the crucial elements for the risk measurement of capital assignments is the probability of default. For a fixed time, t, and a horizon time, b, the probability of default (PD) can be defined as the probability that a credit that has been paid until time t becomes unpaid not later than time $t+b$. The PD is allowed to depend on the credit scoring x, which is, usually, some linear combination of informative covariates of the credit and the clients. Standard methods to estimate the PD include logistic models and other binary response parametric regression models. See the studies of [1,2,3,4,5,6], among others.

In recent years, survival analysis has started to be considered an interesting tool in credit risk problems. Since the work by [7], some literature has been developed following this line. See the works of [8,9,10,11,12,13,14]. Nonparametric estimators of the probability of default based on conditional survival function estimators are presented in [12,13] since the probability of default can be written in terms of the conditional survival function in a right censored context. All these nonparametric estimators are based on covariate smoothing. In the recent work [14], a general nonparametric estimator of the probability of default with double smoothing, both in the covariate and in the time variable, was proposed and studied.

In particular, Beran’s estimator and a doubly smoothed Beran’s estimator of the PD were presented in [13,14], respectively, and their asymptotic properties were analyzed. Simulation studies carried out in these papers show a good performance of the estimators, especially of the doubly smoothed Beran’s estimator. However, these previous studies were carried out using theoretical smoothing parameters. Since the integrated mean square error expressions are complex and depend on several population parameters, they are not useful in practice to obtain plug-in estimations of these theoretical bandwidths. The goal of this work is to propose resampling techniques to approximate them.

Bootstrap has become a strong tool in many statistical applications since it was first introduced by [15]. Bootstrap for right censored data was first proposed by [16] and the bootstrap method and its applications were studied in [17]. Asymptotic theory of bootstrap for right censored data was stablished by [18,19]. In [20], bootstrap for nonparametric regression with right censored observations at fixed covariate values was studied. A bootstrap approach for the nonparametric censored regression setup was studied in [21]. In [22] a local cross-validation bandwidth selector was proposed.

Our approach follows the ideas of [21], and it is based on the obvious bootstrap. Both Beran and smoothed Beran’s estimators are bootstrapped in order to approximate their corresponding optimal bandwidths. The bootstrap is also useful to compute confidence regions, as the existing theoretical result only allows to compute pointwise and theoretical confidence intervals, which are not computable in practice since the variance of the estimator again depends on unknown population quantities. A bootstrap algorithm to compute confidence regions for the probability of default is also proposed.

The remainder of this paper is organized as follows. In Section 2, bootstrap selectors for the bandwidths of Beran and smoothed Beran’s estimators are proposed. In Section 3, a simulation study shows the behavior of the PD estimators with bootstrap bandwidths. The issue of obtaining confidence regions for the probability of default, $PD\left(t\right|x)$, for a fixed value of $x\in I\subseteq \mathbb{R}$ and t covering the interval $I_T\subseteq {\mathbb{R}}^{+}$ is addressed in Section 4 using Beran and smoothed Beran’s estimators. A simulation study on the proposed methods for computing confidence regions is shown in Section 5. In Section 6, Beran’s estimator and the smoothed Beran’s estimator with bootstrap bandwidths are used to estimate the probability of default function conditional on the credit scoring for a German credit dataset. Section 7 contains some concluding remarks.

2. Bandwidth Selection for Beran and Smoothed Beran’s PD Estimators

Consider the right censored simple random sample ${\left\{(X_i,Z_i,{\delta }_i)\right\}}_{i=1}^n$ of $(X,Z,\delta )$ where $X_i$ represents the covariate, $Z_i=min\{T_i,C_i\}$ the observed lifetime and ${\delta }_i=I_{\{T_i\le C_i\}}$ the censoring indicator, where $T_i\ge 0$ and $C_i\ge 0$ are the time to occurrence of the event and the censoring time for the i-th individual of the sample with $i=1,\dots ,n$. In credit risk, usually, X is the credit scoring, Z is the observed maturity, T is the time to default, and C is the time until the end of the study or the anticipated cancellation of the credit. The distribution function of T is denoted by $F\left(t\right)$ and the survival function by $S\left(t\right)=1-F\left(t\right)$. The functions $F\left(t\right|x)$ and $S\left(t\right|x)$ are the conditional distribution and survival functions of T evaluated at t given $X=x$. The conditional distribution function of Z is denoted by $H\left(t\right|x)$, and the conditional distribution function of C is denoted by $G\left(t\right|x)$. It is assumed that an unknown relationship between T and X exists. Given $X=x$, the distributions of T and C conditional to $X=x$ are supposed to be independent.

Let $x\in I$ be a fixed value of the covariate X and $b>0$ a fixed positive time. Then, the probability of default in a time horizon $t+b$ from a maturity time, t, is defined as follows:

$$\begin{array}{ccc}PD\left(t\right|x)& =& P(T\le t+b|T>t,X=x)\hfill \\ & =& {\displaystyle \frac{F(t+b|x)-F(t\left|x\right)}{1-F\left(t\right|x)}}=1-{\displaystyle \frac{S(t+b|x)}{S\left(t\right|x)}}.\hfill \end{array}$$

(1)

In this section, methods for the automatic selection of the bandwidths for Beran’s estimator and the smoothed Beran’s estimator of the probability of default are proposed.

2.1. Beran’s Estimator

Beran’s estimator of the conditional survival function is given by

$$\begin{array}{ccc}{\widehat{S}}_h\left(t\right|x)& =& {\displaystyle \prod _{i=1}^n\left(1-{\displaystyle \frac{I_{\{Z_i\le t,{\delta }_i=1\}}{w}_{h,i}\left(x\right)}{1-{\sum }_{j=1}^n{I}_{\{Z_j<Z_i\}}{w}_{h,j}\left(x\right)}}\right)}\hfill \end{array}$$

(2)

where $w_{h,i}\left(x\right)={\displaystyle \frac{K\left((x-X_i)/h\right)}{{\sum }_{j=1}^{n}K\left((x-X_j)/h\right)}}$ with $i=1,...,n$, $K\left(u\right)$ is a kernel function and $h=h_n$ is the bandwidth that determines the smoothness introduced in the estimator through the covariate X.

Plugging (2) into (1), the probability of default estimator based on Beran’s PD estimator is obtained:

$$\begin{array}{ccc}{\widehat{PD}}_h\left(t\right|x)& =& 1-{\displaystyle \frac{{\widehat{S}}_h(t+b|x)}{{\widehat{S}}_h\left(t\right|x)}}.\hfill \end{array}$$

(3)

In this section, bootstrap methods are proposed for the automatic selection of h. There are two classic methods for bootstrap resampling in a censoring context: the obvious bootstrap and the simple bootstrap. The equivalence between both methods in an unconditional setup is proved in [16]. In [21], this result is extended to the case where a covariate is involved, assuming there is no ties in the sample values of the covariate. This was done by proving the equivalence of the two resampling methods, the obvious bootstrap and the simple weighted bootstrap. In this paper, the following obvious bootstrap method combined with a smoothed bootstrap for the covariate is proposed.

2.1.1. Algorithm for Bootstrap Resampling

Let $I_1\subseteq \mathbb{R}$ be an interval containing appropriate bandwidth values and let $r\in I_1$ be the pilot bandwidth for the bootstrap resampling:

• Obtain $U_1,\dots ,U_n$ iid with $U_i\sim U(0,1)$ and $V_1,\dots ,V_n$ iid with common density K for all $i=1,\dots ,n$.
• For each $i=1,\dots ,n$, define

  $$X_i^{*}=X_{\left[n{U}_i\right]+1}+r{V}_i,$$

  where $\left[u\right]$ is the integer part of u. Generate $T_i^{*}$ from Beran’s estimator of the conditional distribution of T using the sample ${\left\{(X_i,Z_i,{\delta }_i)\right\}}_{i=1}^n$ and bandwidth r, denoted by ${\widehat{F}}_r\left(t\right|X_i^{*})$, and $C_i^{*}$ from the Beran’s estimator of the conditional distribution of C using the sample ${\left\{(X_i,Z_i,1-{\delta }_i)\right\}}_{i=1}^n$ and bandwidth r, denoted by ${\widehat{G}}_r\left(t\right|X_i^{*})$.
  The estimators ${\widehat{F}}_r\left(t\right|X_i^{*})$ and ${\widehat{G}}_r\left(t\right|X_i^{*})$ are forced to be equal to one from the last observed lifetime ($max\{Z_i:i=1,\dots ,n\}$) onwards.
• For each $i=1,\dots ,n$, obtain

  $$Z_i^{*}=min\{T_i^{*},C_i^{*}\},$$

  $${\delta }_i^{*}=I\left(T_i^{*}\le C_i^{*}\right).$$
• Consider the bootstrap resample ${\left\{(X_i^{*},Z_i^{*},{\delta }_i^{*})\right\}}_{i=1}^n$.

In this paper, we want to estimate the probability of default function, $PD\left(t\right|x)$, for a fixed $x\in I$ and t covering the interval $I_T\subset \mathbb{R}$. Therefore, our goal is to get the bandwidth $h_{MISE}\in I_1$ that minimizes the mean integrated squared error given by

$$MIS{E}_x\left(h\right)=E\left({\int }_{I_T}{\left({\widehat{PD}}_h\left(t\right|x)-PD\left(t\right|x)\right)}^{2}dt\right)$$

(4)

whose bootstrap approximation is

$$MIS{E}_x^{*}\left(h\right)=E\left({\int }_{I_T}{\left({\widehat{PD}}_h^{*}\left(t\right|x)-{\widehat{PD}}_r\left(t\right|x)\right)}^{2}dt\right)$$

where ${\widehat{PD}}_r\left(t\right|x)$ is the estimation of the theoretical PD with pilot bandwidth r, using the sample ${\left\{(X_i,Z_i,{\delta }_i)\right\}}_{i=1}^n$ and ${\widehat{PD}}_h^{*}\left(t\right|x)$ is the bootstrap estimation of $PD$ with bandwidth h, using the bootstrap resample ${\left\{(X_i^{*},Z_i^{*},{\delta }_i^{*})\right\}}_{i=1}^n$.

The resampling distribution of ${\widehat{PD}}_h^{*}\left(t\right|x)$ cannot be computed in a close form, so the Monte Carlo method is used. It is based on obtaining B bootstrap resamples and estimating ${\widehat{PD}}_h^{*}\left(t\right|x)$ for each of them. Thus, the distribution of ${\widehat{PD}}_h^{*}\left(t\right|x)$ is approximated by the empirical one of ${\widehat{PD}}_h^{*,1}\left(t\right|x),\dots ,{\widehat{PD}}_h^{*,B}\left(t\right|x)$, obtained from B bootstrap resamples and the bootstrap version of the estimation error committed by Beran’s estimator for any smoothing parameter h is given by

$${MISE}_x^{*}\left(h\right)\simeq {\displaystyle \frac{1}{B}}\sum _{k=1}^B\left({\int }_{I_T}{\left({\widehat{PD}}_h^{*,k}\left(t\right|x)-{\widehat{PD}}_r\left(t\right|x)\right)}^{2}dt\right).$$

(5)

Likewise, the integral is approximated by a Riemann sum.

2.1.2. Algorithm for Bootstrap Bandwidth Selector Based on Beran’s Estimator

Let $x\in I$ be a fixed value of the covariate, $t\in I_T$ and $r\in I_1$:

• Compute ${\widehat{PD}}_r\left(t\right|x)$ from the original sample ${\left\{(X_i,Z_i,{\delta }_i)\right\}}_{i=1}^n$.
• Obtain B bootstrap resamples of the form ${\left\{(X_i^{*,k},Z_i^{*,k},{\delta }_i^{*,k})\right\}}_{i=1}^n$ with $k=1,...,B$ using the smoothed bootstrap with pilot bandwidth $r\in I_1$ and calculate ${\widehat{PD}}_h^{*,k}\left(t\right|x)$ for each of them.
• Approximate $MIS{E}_x^{*}\left(h\right)$ according to (5).
• Repeat Steps 1–3 for values of h in a grid of $I_1$.
• Select the value of h that provides the smallest $MIS{E}_x^{*}\left(h\right)$ as the bootstrap bandwidth $h^{*}$.

Concerning the auxiliary bandwidth $r\in I_1$, a preliminary analysis not shown here suggests

$$r={\displaystyle \frac{3}{4}}\left(Q_X\left(0.975\right)-Q_X\left(0.025\right)\right){\left(\sum _{i=1}^n{\delta }_i\right)}^{-1/3}$$

(6)

where $Q_X\left(u\right)$ is the u quantile of the sample ${\left\{X_i\right\}}_{i=1}^n$, as a suitable pilot bandwidth in this context.

Note that the proposed algorithm is also valid to obtain a bootstrap approximation of the optimal bandwidth for the estimation of $PD\left(t\right|x)$ for fixed values of $t\in I_T$ and $x\in I$ by replacing $MIS{E}_x^{*}\left(h\right)$ by $MS{E}_{t,x}^{*}\left(h\right)$, which is the bootstrap analogue of

$$MS{E}_{t,x}\left(h\right)=E\left({\left({\widehat{PD}}_h\left(t\right|x)-PD\left(t\right|x)\right)}^2\right).$$

2.2. Smoothed Beran’s Estimator

Nonparametric estimators for the probability of default such as Beran’s estimator, ${\widehat{PD}}_h\left(t\right|x)$, are smoothed in the covariate X. It is interesting to consider estimators with double smoothing both in the covariate, X, and in the time variable, T. This idea was previously used in [23,24,25] to obtain doubly smoothed estimators of the conditional survival function. In [14], a doubly smoothed version of PD Beran’s estimator was proposed based on a further smoothing, in time, of Beran’s estimator of the conditional survival function ([26]). Asymptotic properties and simulation studies carried out in [14] show that the doubly smoothed Beran’s estimator performs better than the classical Beran’s estimator when estimating the probability of the default curve. This section presents a method for automatic selection of the bandwidths on which it depends.

Let ${\widehat{S}}_h\left(t\right|x)$ be Beran’s estimator of the conditional survival function given in (2) with $h=h_n$ being the smoothing parameter for the covariate. Then, the expression for the doubly smoothed Beran’s estimator of the conditional survival function defined in [25] is as follows:

$$\begin{array}{ccc}{\tilde{S}}_{h,g}\left(t\right|x)& =& {\displaystyle 1-\sum _{i=1}^n{s}_{\left(i\right)}\mathbb{K}\left({\displaystyle \frac{t-Z_{\left(i\right)}}{g}}\right)}\hfill \end{array}$$

(7)

with $s_{\left(i\right)}={\widehat{S}}_h\left(Z_{(i-1)}\right|x)-{\widehat{S}}_h\left(Z_{\left(i\right)}\right|x)$ where $Z_{\left(i\right)}$ is the i-th element of the sorted sample of Z, $\mathbb{K}\left(t\right)$ is the distribution function of a kernel K, $\mathbb{K}\left(t\right)={\int }_{-\infty }^{t}K\left(u\right)du$, and $g=g_n$ is the smoothing parameter for the time variable. Then, plugging (7) into (1), the probability of the default estimator based on the smoothed Beran’s survival estimator is obtained:

$$\begin{array}{ccc}{\widetilde{PD}}_{h,g}\left(t\right|x)& =& 1-{\displaystyle \frac{{\tilde{S}}_{h,g}(t+b|x)}{{\tilde{S}}_{h,g}\left(t\right|x)}}.\hfill \end{array}$$

(8)

A bootstrap method is proposed for the automatic selection of the bivariate bandwidth $(h,g)$.

2.2.1. Algorithm for Bootstrap Resampling

Let $I_1\subseteq \mathbb{R}$ and $I_2\subseteq \mathbb{R}$ be intervals containing appropriate bandwidth values and let $r\in I_1$ and $s\in I_2$ be pilot bandwidths for the smoothed resample of X, T and C:

• Obtain $U_1,\dots ,U_n$ iid with $U_i\sim U(0,1)$ and $V_1,\dots ,V_n$ iid with common density K, $W_1^1,\dots ,W_n^1$ iid with common density K and $W_1^2,\dots ,W_n^2$ iid with common density K for all $i=1,\dots ,n$.
• For each $i=1,\dots ,n$, obtain

  $$X_i^{*}=X_{\left[n{U}_i\right]+1}+r{V}_i,$$

  $$T_i^{*}=T_{0,i}^{*}+s{W}_i^1$$

  $$C_i^{*}=C_{0,i}^{*}+s{W}_i^2$$

  where $T_{0,i}^{*}$ is resampled from ${\widehat{F}}_r\left(t\right|X_i^{*})$ constructed using Beran’s estimator with the sample ${\left\{(X_i,Z_i,{\delta }_i)\right\}}_{i=1}^n$ and $C_{0,i}^{*}$ is resampled from ${\widehat{G}}_r\left(t\right|X_i^{*})$ constructed using Beran’s estimator with the sample ${\left\{(X_i,Z_i,1-{\delta }_i)\right\}}_{i=1}^n$.
• For each $i=1,\dots ,n$, obtain

  $$Z_i^{*}=min\{T_i^{*},C_i^{*}\},$$

  $${\delta }_i^{*}=I\left(T_i^{*}\le C_i^{*}\right).$$
• Consider the bootstrap resample ${\left\{(X_i^{*},Z_i^{*},{\delta }_i^{*})\right\}}_{i=1}^n$.

The conditional distribution functions of $T^{*}|X^{*}$ and $C^{*}|X^{*}$ are, respectively, the smoothed Beran’s estimators ${\tilde{F}}_{r,s}\left(t\right|X_i^{*})$ and ${\tilde{G}}_{r,s}\left(t\right|X_i^{*})$.

The optimal bivariate bandwidth, $(h_{MISE},g_{MISE})\in I_1\times I_2$ is defined as the pair of bandwidths that minimizes the mean integrated squared error given by

$$MIS{E}_x(h,g)=E\left({\displaystyle {\int }_{I_T}{\left({\widetilde{PD}}_{h,g}\left(t\right|x)-PD\left(t\right|x)\right)}^{2}dt}\right).$$

(9)

The bootstrap version of $MIS{E}_x(h,g)$ is given by

$$MIS{E}_x^{*}(h,g)=E\left({\int }_{I_T}{\left({\widetilde{PD}}_{h,g}^{*}\left(t\right|x)-{\widetilde{PD}}_{r,s}\left(t\right|x)\right)}^{2}dt\right),$$

where ${\widetilde{PD}}_{r,s}\left(t\right|x)$ is the smoothed Beran’s PD estimation with pilot bandwidths $(r,s)\in I_1\times I_2$ using the sample ${\left\{(X_i,Z_i,{\delta }_i)\right\}}_{i=1}^n$ and ${\widetilde{PD}}_{h,g}^{*}\left(t\right|x)$ is the bootstrap estimation of $PD$ with bandwidths $(h,g)$, using the bootstrap resample ${\left\{(X_i^{*},Z_i^{*},{\delta }_i^{*})\right\}}_{i=1}^n$. Since the sampling distribution of ${\widetilde{PD}}_{h,g}^{*}\left(t\right|x)$ is unknown, the Monte Carlo method gives the following approximation

$${MISE}_x^{*}(h,g)\simeq {\displaystyle \frac{1}{B}}\sum _{k=1}^B\left({\int }_{I_T}{\left({\widetilde{PD}}_{h,g}^{*,k}\left(t\right|x)-{\widetilde{PD}}_{r,s}\left(t\right|x)\right)}^{2}dt\right),$$

(10)

based on the empirical distribution of ${\widetilde{PD}}_{h,g}^{*}\left(t\right|x)$ obtained from B bootstrap resamples. The integral is approximated by a Riemann sum.

2.2.2. Algorithm for Bootstrap Bandwidth Selector Based on the Smoothed Beran’s Estimator

Let x be a fixed value of the covariate, $t\in I_T$ and $(r,s)\in I_1\times I_2$:

• Compute ${\widetilde{PD}}_{r,s}\left(t\right|x)$ from the original sample ${\left\{(X_i,Z_i,{\delta }_i)\right\}}_{i=1}^n$.
• Obtain B bootstrap resamples of the form ${\left\{(X_i^{*,k},Z_i^{*,k},{\delta }_i^{*,k})\right\}}_{i=1}^n$ with $k=1,...,B$ using the doubly smoothed bootstrap and calculate ${\widetilde{PD}}_{h,g}^{*,k}\left(t\right|x)$ for each of them.
• Approximate $MIS{E}_x^{*}\left(h\right)$ according to (10).
• Repeat Steps 1–3 for pairs of values $(h,g)$ in a grid of $I_1\times I_2$.
• Obtain the pair $(h,g)$ that provides the smallest $MIS{E}_x^{*}(h,g)$ as the bootstrap bandwidth $(h^{*},g^{*})$.

The auxiliary bandwidth $r\in I_1$ was defined in (6). The pilot bandwidth $s\in I_2$ for the time variable smoothing is

$$s={\displaystyle \frac{3}{4}}\left(Q_Z\left(0.975\right)-Q_Z\left(0.025\right)\right){\left(\sum _{i=1}^n{\delta }_i\right)}^{-1/7}$$

(11)

where $Q_Z\left(u\right)$ is the u quantile of the sample ${\left\{Z_i\right\}}_{i=1}^n$.

3. Simulation Study for Bandwidth Selection

A simulation study was conducted in order to show the behavior of bootstrap bandwidth selectors for Beran’s and smoothed Beran’s estimators proposed in Section 2. Two models are considered, one with Weibull lifetime and censoring time distributions and another one with exponential distributions.

Model 1 considers a $U(0,1)$ distribution for X. The time to occurrence of the event conditional to the covariate, ${T|}_{X=x}$, follows a Weibull distribution with parameters $d=2$ and $\Gamma {\left(x\right)}^{-1/d}$ where $\Gamma \left(x\right)=1+5x$, and the censoring time conditional to the covariate, ${C|}_{X=x}$, follows a Weibull distribution with parameters $d=2$ and $\Delta {\left(x\right)}^{-1/d}$ where $\Delta \left(x\right)=10+d_{1}x+20x^2$. In this case, the conditional survival function and the censoring conditional probability are given by

$$S\left(t\right|x)=e^{-\Gamma \left(x\right)t^d},$$

$$P(\delta =0|X=x)={\displaystyle \frac{\Delta \left(x\right)}{\Gamma \left(x\right)+\Delta \left(x\right)}}.$$

Having set the value of the covariate $x=0.6$, the value of $d_1$ is chosen so that the censoring conditional probability is $0.2$ and $0.5$. These values are $d_1=-27$ and $d_1=-22$, respectively. The conditional survival function for this model is estimated in a time grid of size $n_T$, $0<t_1<\cdots <t_{n_T}$, where $t_{n_T}+b=F^{-1}\left(0.95\right|x)=0.8654$ and $b=0.15$, i.e., about $20\%$ of the time grid range for the value of the covariate $x=0.6$. Therefore, in this case, $I_T=(0,0.8654)$.

Model 2 also considers a $U(0,1)$ distribution for X. The time to occurrence of the event conditional to the covariate, ${T|}_{X=x}$, follows an exponential distribution with parameter $E\left(x\right)=2+58x-160x^2+107x^3$, and the censoring time conditional to the covariate, ${C|}_{X=x}$, follows an exponential distribution with parameter $\Theta \left(x\right)=10+c_{1}x+20x^2$. In this scenario, the conditional survival function and the censoring conditional probability are the following:

$$S\left(t\right|x)=e^{-E\left(x\right)t},$$

$$P(\delta =0|X=x)={\displaystyle \frac{\Theta \left(x\right)}{E\left(x\right)+\Theta \left(x\right)}}.$$

Having set the value of the covariate, $x=0.8$, the value of $c_1$ is chosen so that the censoring conditional probability is $0.2$ and $0.5$. These values are $c_1=-113/4$ and $c_1=-55/2$, respectively. The conditional survival function is estimated in a time grid of size $n_T$, $0<t_1<\cdots <t_{n_T}$, where $t_{n_T}+b=F^{-1}\left(0.95\right|x)=3.8211$ and $b=0.7$, i.e., about $20\%$ of the time grid range for the value of the covariate $x=0.8$. Therefore, in this case, $I_T=(0,3.8211)$.

It can be proved that Model 1 is close to a proportional hazards model, while Model 2 moves away from this parametric model. These two models were used in the simulation study carried out by [14].

The boundary effect is corrected using the reflexion principle, and the truncated Gaussian kernel with a truncation range $(-50,50)$ is considered. The size of the lifetime grid is $n_T=100$. The sample size is $n=400$. The simulation study is carried out with software developed in R by the authors themselves. In order to minimize the MISE error function without increasing CPU time more than necessary, a limited-memory algorithm for solving large nonlinear optimization problems is used, L-BFGS-B. It was proposed by [27] for solving optimization problems subject to simple bounds on the variables in which information on the Hessian matrix is difficult to obtain. Results of numerical studies about this method are shown in [27]. It is available at the stats package from the Comprehensive R Archive Network (CRAN) using Fortran 77 subroutines (see [28]).

3.1. Simulation Study for Beran’s Estimator

In this subsection, the behavior of the bootstrap bandwidth selector for Beran’s estimator is shown. For each model, the estimation error function $MIS{E}_x\left(h\right)$ is approximated via Monte Carlo using 300 simulated samples. The bandwidth that minimises $MIS{E}_x\left(h\right)$ is obtained and denoted by $h_{MISE}$. The values of $h_{MISE}$ and $MIS{E}_x\left(h_{MISE}\right)$ are used as a benchmark.

In the simulation study, $N=300$ simulated samples are used. For each sample, $B=500$ bootstrap resamples are obtained to approximate the bootstrap MISE function, $MIS{E}_x^{*}\left(h\right)$, and obtain the bootstrap bandwidth associated to each simulated sample $h_j^{*}$, $j=1,2,\dots ,N$. The mean value of the N bootstrap bandwidths and the standard deviation are defined as follows

$$\overline{h^{*}}={\displaystyle \frac{1}{N}}\sum _{j=1}^N{h}_j^{*},sd\left(h^{*}\right)=\sqrt{{\displaystyle \frac{1}{N}}\sum _{j=1}^N{\left(h_j^{*}-\overline{h^{*}}\right)}^2}$$

As a relative measure of the difference between the bootstrap bandwidth and the optimal one, we compute

$$H_j^{*}={\displaystyle \frac{h_j^{*}-h_{MISE}}{h_{MISE}}}$$

with $j=1,\dots ,N$. The mean of the absolute value of these relative deviations, $\overline{H^{*}}=\frac{1}{N}{\sum }_{j=1}^N\left|H_j^{*}\right|$, is a good measure of how close the bootstrap bandwidth is to the optimal one.

For each sample, the estimation error committed by Beran’s estimator with the corresponding bootstrap bandwidth,

$${\displaystyle MIS{E}_x\left(h_j^{*}\right)=E\left({\int }_{I_T}{\left({\widehat{PD}}_{h_j^{*}}\left(t\right|x)-PD\left(t\right|x)\right)}^{2}dt\right),}$$

and its squared root, $RMIS{E}_x\left(h_j^{*}\right)$, are approximated via Monte Carlo using 300 simulated samples. The mean of these estimation errors given by

$$\overline{RMIS{E}_x\left(h^{*}\right)}={\displaystyle \frac{1}{N}}\sum _{j=1}^{N}RMIS{E}_x\left(h_j^{*}\right)$$

is used as a measure of the estimation error made by the bootstrap bandwidth, when compared with the estimation error made by the MISE bandwidth.

As a relative measure of the difference between the estimation errors using the bootstrap and the MISE bandwidths, the following ratios are defined:

$$R_j^{*}={\displaystyle \frac{RMIS{E}_x\left(h_j^{*}\right)-RMISE\left(h_{MISE}\right)}{RMISE\left(h_{MISE}\right)}}$$

satisfying $R_j^{*}\ge 0$ for all $j=1,\dots ,N$. The mean of the $R_j^{*}$ values with $j=1,\dots ,N$ is denoted by $\overline{R^{*}}=\frac{1}{N}{\sum }_{j=1}^N{R}_j^{*}$. Small values (close to zero) of $\overline{H^{*}}$ and $\overline{R^{*}}$ indicate good behavior of the bootstrap bandwidth. Values of the bootstrap bandwidths, estimation errors and relative measures for Models 1 and 2 are included in

Table 1. MISE, average bootstrap bandwidths and estimation errors of Beran’s PD estimator in each level of censoring conditional probability for Models 1 and 2. Numbers within brackets are standard deviations.

	Model 1		Model 2
$P(\delta =0|X=x)$	$0.2$	$0.5$	$0.2$	$0.5$
$h_{MISE}$	0.37576	0.35909	0.09494	0.10959
$RMIS{E}_x\left(h_{MISE}\right)$	0.05520	0.11144	0.27942	0.49991
$\overline{h^{*}}\left(sd\right)$	0.27856 (0.092)	0.30892 (0.110)	0.21763 (0.041)	0.23091 (0.068)
$\overline{H^{*}}$	0.31431	0.29306	1.29211	1.10692
$\overline{RMIS{E}_x\left(h^{*}\right)}$	0.05700	0.11405	0.29671	0.50824
$\overline{R^{*}}$	0.03260	0.02336	0.06188	0.01666

. The results show a good performance of the proposed bootstrap selector.

Figure 1 and Figure 2 show the function $MIS{E}_x\left(h\right)$ along with the Monte Carlo approximations of $MIS{E}_x^{*}\left(h\right)$ for some simulated samples and the boxplots of $H_j^{*}$ and $R_j^{*}$ with $j=1,\dots ,N$ for Models 1 and 2. The method tends to slightly underestimate the value of $h^{*}$ with respect to $h_{MISE}$ in Model 1 and overestimate its value in Model 2, which is reflected in the boxplots of $H_j^{*}$. Nevertheless, these figures show that the $MIS{E}_x\left(h\right)$ curve is fairly flat and variations in the selection of h do not imply an important increase in the estimation error.

In order to illustrate the results, Figure 3 shows the theoretical probability of default function and Beran’s estimation with the MISE and bootstrap bandwidths drawn for one sample from Model 1 and 2 when the conditional probability of censoring is $0.5$. For large values of time, the performance of the estimator becomes worse, due to the fact that in that region there are few data, most of them censored, and therefore offering poor information.

3.2. Simulation Study for the Smoothed Beran’s Estimator

In this section, a simulation study on the bootstrap bandwidth selector of the smoothed Beran’s estimator in (8) is carried out. The resampling technique and Monte Carlo approximation of the MISE presented in Section 2.2 are used.

For each model, the error function $MIS{E}_x(h,g)$ is approximated via Monte Carlo from 300 simulated samples, and the bivariate bandwidth that minimizes $MIS{E}_x(h,g)$ is obtained and denoted by $(h_{MISE},g_{MISE})$. The values of $(h_{MISE},g_{MISE})$ and $MIS{E}_x(h_{MISE},g_{MISE})$ are used as a benchmark.

In the study, $N=300$ samples are simulated. For each simulated sample, the corresponding bootstrap bandwidths are approximated from $B=500$ resamples, obtaining $(h_j^{*},g_j^{*})$ with $j=1,\dots ,N$. The mean value of the N bootstrap bandwidths and the standard deviation are the following:

$$(\overline{h^{*}},\overline{g^{*}})=\left({\displaystyle \frac{1}{N}}\sum _{j=1}^N{h}_j^{*},{\displaystyle \frac{1}{N}}\sum _{j=1}^N{g}_j^{*}\right),$$

$$sd\left(h^{*}\right)=\sqrt{{\displaystyle \frac{1}{N}}\sum _{j=1}^N{\left(h_j^{*}-\overline{h^{*}}\right)}^2},sd\left(g^{*}\right)=\sqrt{{\displaystyle \frac{1}{N}}{\sum }_{j=1}^N{\left(g_j^{*}-\overline{g^{*}}\right)}^2}.$$

In order to measure the distance of the bootstrap bidimensional bandwidth of the j-th sample, $(h_j^{*},g_j^{*})$, from the corresponding MISE bandwidth, $(h_{MISE},g_{MISE})$, consider the vector

$$D_j^{*}=\left({\displaystyle \frac{h_j^{*}-h_{MISE}}{h_{MISE}}},{\displaystyle \frac{g_j^{*}-g_{MISE}}{g_{MISE}}}\right)\in {\mathbb{R}}^2.$$

and its Euclidean norm denoted by $H_j^{*}={\parallel D_j^{*}\parallel }_2$ with $j=1,\dots ,N$. The mean value, $\overline{H^{*}}=\frac{1}{N}{\sum }_{j=1}^N{H}_j^{*}$ is a measure of how close the bootstrap bandwidths are to the MISE one.

For each sample, the estimation error committed by the smoothed Beran’s estimator with the corresponding bootstrap bandwidth,

$${\displaystyle MIS{E}_x(h_j^{*},g_j^{*})=E\left({\int }_{I_T}{\left({\widetilde{PD}}_{h_j^{*},g_j^{*}}\left(t\right|x)-PD\left(t\right|x)\right)}^{2}dt\right),}$$

and its squared root, $RMIS{E}_x(h_j^{*},g_j^{*})$, are approximated via Monte Carlo using 300 simulated samples. The mean of these estimation errors given by

$$\overline{RMIS{E}_x(h^{*},g^{*})}={\displaystyle \frac{1}{N}}\sum _{j=1}^{N}RMIS{E}_x\left(h_j^{*},g_j^{*}\right)$$

is used as a measure of the estimation error committed by the bootstrap bidimensional bandwidth in the model.

The ratio

$$R_j^{*}={\displaystyle \frac{RMIS{E}_x(h_j^{*},g_j^{*})-RMIS{E}_x(h_{MISE},g_{MISE})}{RMIS{E}_x(h_{MISE},g_{MISE})}}$$

is defined as a relative measure of the difference between the error committed by the estimator with bootstrap bandwidth and MISE bandwidth. The mean of the positive values $R_j^{*}$ with $j=1,\dots ,N$ is denoted by $\overline{R^{*}}=\frac{1}{N}{\sum }_{j=1}^N{R}_j^{*}$. Values of the bootstrap bivariate bandwidths, estimation errors and relative measures for Models 1 and 2 are included in

Table 2. MISE, average bootstrap bandwidths and estimation errors of the smoothed Beran’s PD estimator in each level of censoring conditional probability for Models 1 and 2. Numbers within brackets are standard deviations.

	Model 1		Model 2
$P(\delta =0|X=x)$	$0.2$	$0.5$	$0.2$	$0.5$
$h_{MISE}$	0.21633	0.16735	0.11122	0.37551
$g_{MISE}$	0.09286	0.14612	1.27755	1.68878
$RMIS{E}_x(h_{MISE},g_{MISE})$	0.03710	0.05094	0.09829	0.12322
$\overline{h^{*}}\left(sd\right)$	0.11736 (0.051)	0.11219 (0.057)	0.19813 (0.180)	0.16593 (0.218)
$\overline{g^{*}}\left(sd\right)$	0.12647 (0.039)	0.19671 (0.054)	0.60005 (0.375)	1.45428 (0.711)
$\overline{H^{*}}$	0.68121	0.63604	1.22609	0.89164
$\overline{RMIS{E}_x(h^{*},g^{*})}$	0.04620	0.06793	0.22135	0.28342
$\overline{R^{*}}$	0.24517	0.33357	1.25199	1.30003

.

Figure 4 shows the $MIS{E}_x(h,g)$ function of the smoothed Beran’s estimator and its bootstrap approximation for one sample of both Models 1 and 2 when the conditional probability of censoring is $0.5$. It is approximated on a grid of 50 values of h and 50 values of g. Note that both $MIS{E}_x(h,g)$ and $MIS{E}_x^{*}(h,g)$ curves for each fixed h value are quite similar in the region close to the minimum value of $MIS{E}_x^{*}(h,g)$. Thus, the influence of covariate smoothing parameter h is weak when estimating the PD using values of bandwidth g close to the optimal one.

Figure 5 and Figure 6 show the boxplots of $H_j^{*}$ and $R_j^{*}$ with $j=1,\dots ,N$. In general, the selector tends to underestimate the value of the bandwidths. Due to the behavior of the $MIS{E}_x(h,g)$ curves mentioned above, this does not lead to a significant increase in the estimation error.

Figure 7 shows the theoretical probability of the default function and Beran’s estimation with MISE and bootstrap bandwidths for one sample from Models 1 and 2 when the conditional probability of censoring is $0.5$. Comparing this figure with the equivalent one for Beran’s estimator shown in Figure 3, the improvement in estimation due to the double smoothing is remarkable.

The results showed in Table 1 and Table 2 are summarized in

Table 3. Comparative table of the estimation error of Beran’s estimator and the smoothed Beran’s estimator in Models 1 and 2.

		Model 1		Model 2
$P(\delta =0|X=x)$		$0.2$	$0.5$	$0.2$	$0.5$
BERAN	$\overline{RMIS{E}_x\left(h^{*}\right)}$	0.05579	0.11206	0.28593	0.49916
SBERAN	$\overline{RMISE(h^{*},g^{*})}$	0.04629	0.07216	0.20007	0.27611
$R_S$		0.20523	0.55294	0.42915	0.80783

to compare the behavior of Beran (BERAN) and the smoothed Beran’s (SBERAN) estimators for the PD and to evaluate whether the improvement that smoothing in the time variable provides for PD estimation is preserved when approximating the smoothing parameters by resampling techniques. Table 3 shows the estimation errors committed by Beran and the smoothed Beran’s estimators of the probability of default using bootstrap bandwidths. In order to measure the increase in estimation error resulting from using Beran’s estimator, the following ratio is defined:

$$R_S={\displaystyle \frac{\overline{RMIS{E}_x\left(h^{*}\right)}-\overline{RMISE(h^{*},g^{*})}}{\overline{RMISE(h^{*},g^{*})}}}$$

and included in Table 3.

In Model 1, the estimation error committed by Beran’s estimator is $20\%$ larger than the error committed by the smoothed Beran’s estimator when the conditional probability of censoring is 0.2 and by $50\%$ when the conditional probability of censoring is 0.5. In Model 2, these differences are even more significant: the estimation error increases up to $80\%$ when using Beran’s estimator with bootstrap bandwidth instead of the doubly smoothed Beran’s estimator.

4. Confidence Regions Using Beran and Smoothed Beran’s Estimators

Let $x\in I$ be a fixed value of the covariate and consider $PD\left(t\right|x)$ the probability of default curve with $t\in I_T$. The curve $PD\left(t\right|x)$ belongs to the function space $\mathcal{F}\left(I_T\right)$ whose elements are real-valued functions with domain $I_T$. From the sample $\{(X_i,Z_i,{\delta }_i),$$i=1,...,n\}$, Beran’s estimation of $PD\left(t\right|x)$, ${\widehat{PD}}_h\left(t\right|x)$, is obtained and a confidence region of $PD\left(t\right|x)$ at $1-\alpha$ confidence level associated to Beran’s estimator can be constructed. A similar construction is done for the smoothed Beran’s estimator. This confidence region of $PD\left(t\right|x)$ is a random subset of $\mathcal{F}\left(I_T\right)$ denoted by $R_{\alpha }$ that satisfies

$$P(PD\left(t\right|x)\in R_{\alpha },\forall t\in I_T)=1-\alpha .$$

In this section, a method for constructing confidence regions, $R_{\alpha }$, based on Beran and the smoothed Beran’s estimator is developed.

First, Beran’s estimator of the probability of default, ${\widehat{PD}}_h\left(t\right|x)$, given in (3) is used. This method follows the ideas of [29] to obtain prediction regions. It is based on finding the value of ${\lambda }_{\alpha }\in {\mathbb{R}}^{+}$ such that

$$P\left(|{\widehat{PD}}_h\left(t\right|x)-PD\left(t\right|x)|<{\lambda }_{\alpha }\sigma \left(t\right),\forall t\in I_T\right)=1-\alpha$$

with ${\sigma }^2\left(t\right)=Var\left({\widehat{PD}}_h\left(t\right|x)\right)$. Thus, the theoretical confidence region is defined by

$$R_{\alpha }=\left\{\left({\widehat{PD}}_h\left(t\right|x)-{\lambda }_{\alpha }\sigma \left(t\right),{\widehat{PD}}_h\left(t\right|x)+{\lambda }_{\alpha }\sigma \left(t\right)\right):t\in I_T\right\}.$$

Since ${\lambda }_{\alpha }$ and $\sigma \left(t\right)$ are unknown, they are approximated by means of a bootstrap technique. The bootstrap confidence region is defined as follows:

$$R_{\alpha }^{*}=\left\{\left({\widehat{PD}}_h^{*}\left(t\right|x)-{\lambda }_{\alpha }^{*}{\sigma }^{*}\left(t\right),{\widehat{PD}}_h^{*}\left(t\right|x)+{\lambda }_{\alpha }^{*}{\sigma }^{*}\left(t\right)\right):t\in I_T\right\}.$$

where ${\widehat{PD}}_h^{*}\left(t\right|x)$ is the bootstrap estimation of $PD$ with bandwidth h and ${\lambda }_{\alpha }^{*}$ and ${\sigma }^{*}\left(t\right)$ are the bootstrap analogue of ${\lambda }_{\alpha }$ and $\sigma \left(t\right)$. The confidence region $R_{\alpha }^{*}$ satisfies

$$p\left({\lambda }_{\alpha }^{*}\right)=P\left({\widehat{PD}}_r\left(t\right|x)\in R_{\alpha }^{*},\forall t\in I_T\right)=1-\alpha .$$

(12)

From the original sample ${\left\{(X_i,Z_i,{\delta }_i)\right\}}_{i=1}^n$, Beran’s estimator of $PD\left(t\right|x)$ is obtained with appropriate bandwidth h, ${\widehat{PD}}_h\left(t\right|x)$. The algorithm to obtain the bootstrap confidence region for $PD\left(t\right|x)$ at confidence level $1-\alpha$ associated to ${\widehat{PD}}_h\left(t\right|x)$ is explained below. The Monte Carlo method is used to approximate ${\sigma }^{*}\left(t\right)$, and an iterative method is used to approximate the value of ${\lambda }_{\alpha }^{*}$ so that the confidence region has a confidence level approximately equal to $1-\alpha$.

4.1. Confidence Region Based on Beran’s Estimator

• Compute Beran’s estimator ${\widehat{PD}}_r\left(t\right|x)$ from the original sample ${\left\{(X_i,Z_i,{\delta }_i)\right\}}_{i=1}^n$ and pilot bandwidth $r\in I_1$.
• Generate B bootstrap resamples of the form ${\left\{(X_i^{*,k},Z_i^{*,k},{\delta }_i^{*,k})\right\}}_{i=1}^n$ by means of the resampling algorithm presented in SubSection 2.1 and pilot bandwidth r.
• For $k=1,\dots ,B$, compute ${\widehat{PD}}_h^{*,k}\left(t\right|x)$ with the k-th bootstrap resample and bandwidth h, obtaining ${\left\{{\widehat{PD}}_h^{*,k}\left(t\right|x)\right\}}_{k=1}^B$.
• Approximate the standard deviation of ${\widehat{PD}}_h^{*}\left(t\right|x)$ by

  $${\sigma }^{*}\left(t\right)\simeq {\left({\displaystyle \frac{1}{B}}\sum _{k=1}^B{\left({\widehat{PD}}_h^{*,k}\left(t\right|x)-{\displaystyle \frac{1}{B}}\sum _{\mathrm{l}=1}^B{\widehat{PD}}_h^{*,\mathrm{l}}\left(t\right|x)\right)}^2\right)}^{1/2},t\in I_T.$$
• Use an iterative method to obtain an approximation of the value ${\lambda }_{\alpha }^{*}$ defined in (12).
• The confidence region is given by

  $$R_{\alpha }=\left\{\left({\widehat{PD}}_h\left(t\right|x)-{\lambda }_{\alpha }^{*}{\sigma }^{*}\left(t\right),{\widehat{PD}}_h\left(t\right|x)+{\lambda }_{\alpha }^{*}{\sigma }^{*}\left(t\right)\right):t\in I_T\right\}.$$

4.2. Iterative Method to Approximate ${\lambda }_{\alpha }^{*}$

The iterative method to approximate the value of ${\lambda }_{\alpha }^{*}\in {\mathbb{R}}^{+}$ so that the confidence region $R_{\alpha }^{*}$ has a confidence level approximately equal to $1-\alpha$ is explained below. This algorithm allows the parameter ${\lambda }_{\alpha }^{*}$ to be approximated quickly and efficiently.

Let ${\left\{{\widehat{PD}}_h^{*,k}\left(t\right|x)\right\}}_{k=1}^B$ be the Beran’s estimations of the PD with bandwidth h over a set of B bootstrap resamples of ${\left\{(X_i,Z_i,{\delta }_i)\right\}}_{i=1}^n$. Define the Monte Carlo approximation of $p\left(\lambda \right)$ in (12), for any $\lambda \in {\mathbb{R}}^{+}$, as follows:

$$p\left(\lambda \right)\simeq {\displaystyle \frac{1}{B}}\sum _{k=1}^{B}I\left\{{\widehat{PD}}_r\left(t\right|x)\in \left({\widehat{PD}}_h^{*,k}\left(t\right|x)-\lambda {\sigma }^{*}\left(t\right),{\widehat{PD}}_h^{*,k}\left(t\right|x)+\lambda {\sigma }^{*}\left(t\right)\right),\forall t\in I_T\right\}.$$

(13)

Let ${\lambda }_L,{\lambda }_H\in {\mathbb{R}}^{+}$ be such that $p\left({\lambda }_L\right)\le 1-\alpha \ge p\left({\lambda }_H\right)$ and let $\zeta >0$ be a tolerance, for example, $\zeta ={10}^{-4}$.

• Obtain ${\lambda }_M={\displaystyle \frac{{\lambda }_L+{\lambda }_H}{2}}$ and compute Monte Carlo approximations of $p\left({\lambda }_L\right)$, $p\left({\lambda }_M\right)$ and $p\left({\lambda }_H\right)$ according to (13).
• If $p\left({\lambda }_M\right)=1-\alpha$ or $p\left({\lambda }_H\right)-p\left({\lambda }_L\right)<\zeta$, then ${\lambda }_{\alpha }^{*}={\lambda }_M$. Otherwise,

  2.1

  If $1-\alpha <p\left({\lambda }_M\right)$, then ${\lambda }_H={\lambda }_M$ and return to Step 1.

  2.2

  If $p\left({\lambda }_M\right)<1-\alpha$, then ${\lambda }_L={\lambda }_M$ and return to Step 1.

A preliminary analysis not shown here suggests the following pilot bandwidth:

$$r={\displaystyle \frac{3}{4}}\left(Q_X\left(0.975\right)-Q_X\left(0.025\right)\right){\left(\sum _{i=1}^n{\delta }_i\right)}^{-1/3}.$$

This method to obtain confidence regions for the curve $PD\left(t\right|x)$ for fixed $x\in I$ and t covering $I_T$ based on Beran’s estimator can be adapted to obtain confidence regions using the doubly smoothed Beran’s estimator. Simply replace Beran’s estimator ${\widehat{PD}}_h\left(t\right|x)$ by the smoothed Beran’s estimator ${\widetilde{PD}}_{h,g}\left(t\right|x)$ given in (8) where necessary, and obtain the analogous bootstrap approximations of ${\lambda }_{\alpha }$ and $\sigma \left(t\right)$. The confidence region is given by

$$R_{\alpha }=\left\{\left({\widetilde{PD}}_{h,g}\left(t\right|x)-{\lambda }_{\alpha }^{*}{\sigma }^{*}\left(t\right),{\widetilde{PD}}_{h,g}\left(t\right|x)+{\lambda }_{\alpha }^{*}{\sigma }^{*}\left(t\right)\right):t\in I_T\right\}.$$

Denote the lower and upper bounds of the confidence region by $l(t,x)$ and $u(t,x)$, respectively. It may happen that the lower bound of the confidence region is less than 0 or the upper bound is greater than 1 for some points $(t_0,x_0)$. When this happens, we set $l(t_0,x_0)=0$ or $u(t_0,x_0)=1$, as appropriate.

The pilot bandwidths defined in (6) and (11) are used for the confidence region algorithm based on both Beran and smoothed Beran’s estimators.

5. Simulation Study for Confidence Regions

A simulation study is carried out to test the performance of bootstrap confidence regions proposed. Models 1 and 2 described in Section 3 are considered in this study, with identical features. The methods shown in Section 4 are used for this purpose with both Beran and smoothed Beran’s estimators. When Beran’s estimator is used, the bandwidth that minimizes the mean integrated squared error, $h=h_{MISE}$, is used. Similarly, if the smoothed Beran’s estimator is used, the two-dimensional bandwidth that minimizes the mean integrated squared error, $(h,g)=(h_{MISE},g_{MISE})$, is used. These bandwidths are unknown in practice, but they allow a fair comparison of methods in the simulation study.

The simulation set-up is the one explained in Section 3. Two conditional probabilities of censoring are considered for each model: $P(\delta =0|x)=0.2$ and $P(\delta =0|x)=0.5$. The number of bootstrap resamples of each samples is $B=500$, and $N=300$ simulated samples of each model are obtained. The sample size is $n=400$. The confidence level is $1-\alpha$ with $\alpha =0.05$.

Figure 8 shows Beran’s estimations of the PD for $B=500$ resamples from one sample of Model 2 when the conditional probability of censoring is 0.5. The theoretical probability of default is also plotted in the figure. The PD is estimated on a time grid $t_1=0<t_2<\cdots <t_{n_T}$ such that $t_{n_T}+b=F^{-1}\left(0.95\right|x)$. The information provided by the data in the right tail of such a time distribution is sparse due to high censoring. The method results in extremely wide confidence regions or degeneration to zero as in the case of Model 2 (see Figure 8). Therefore, the time grid in this section is restricted to the interval where sufficient information is available.

For this section, we consider the problem of obtaining the bootstrap confidence region for the probability of default in a time grid $t_1=0<t_2<\cdots <t_{n_T}$ such that with $t_k\in I_T\subseteq {\mathbb{R}}^{+}$ for all $k=1,\dots ,n_T$ and $t_{n_T}+b=F^{-1}\left(0.70\right|x)$, with b being approximately equal to $20\%$ of the grid length. For Model 1, having set the value of the covariate, $x=0.6$, the time horizon is $b=0.1$ ($20\%$ of the time range) and $t_{n_T}+b=F^{-1}\left(0.70\right|x=0.6)=0.55$. For Model 2, having set the value of the covariate, $x=0.8$, the time horizon is $b=0.3$ ($20\%$ of the time range) and $t_{n_T}+b=F^{-1}\left(0.70\right|x=0.8)=1.55$.

Table 4. MISE bandwidths and RMISE of Beran and smoothed Beran’s estimators in each level of censoring conditional probability for Models 1 and 2 when $t_{n_T}+b=F^{-1}\left(0.70\right|x)$.

		Model 1		Model 2
$P(\delta =0|X=x)$		$0.2$	$0.5$	$0.2$	$0.5$
BERAN	$h_{MISE}$	0.375510	0.320408	0.041837	0.057755
	$RMIS{E}_x\left(h_{MISE}\right)$	0.019403	0.025943	0.193334	0.220733
SBERAN	$h_{MISE}$	0.230612	0.196939	0.094286	0.154490
	$g_{MISE}$	0.073673	0.083469	0.908163	1.071429
	$RMIS{E}_x(h_{MISE},g_{MISE})$	0.013658	0.018165	0.026161	0.029007

contains the bandwidths that minimize the MISE function for Beran’s estimator and the smoothed Beran’s estimator along this new time grid.

For each model, the confidence region is obtained according to the method explained in Section 4 using both Beran’s estimator and the smoothed Beran’s estimator. The criteria for comparing the methods are set out below.

A confidence region performs well if its coverage is close to the nominal one, in this case $1-\alpha =0.95$, and has a small area or average width. Denoting $l(t,x)={\widehat{PD}}_h\left(t\right|x)-{\lambda }_{\alpha }^{*}{\sigma }^{*}\left(t\right)$ and $u(t,x)={\widehat{PD}}_h\left(t\right|x)+{\lambda }_{\alpha }^{*}{\sigma }^{*}\left(t\right)$ when using Beran’s estimator or $l(t,x)={\widetilde{PD}}_{h,g}\left(t\right|x)-{\lambda }_{\alpha }^{*}{\sigma }^{*}\left(t\right)$ and $u(t,x)={\widetilde{PD}}_{h,g}\left(t\right|x)+{\lambda }_{\alpha }^{*}{\sigma }^{*}\left(t\right)$ when using the smoothed Beran’s estimator, the following values measure the performance of the confidence region and allow for the comparison of results.

Coverage is the percentage of bootstrap regions that contain the whole theoretical probability of default curve and it is defined as follows

$${\displaystyle \frac{1}{N}}\sum _{j=1}^{N}I\left\{PD\left(t_k\right|x)\in \left(l(t_k,x),u(t_k,x)\right),\forall k=1,...,n_T\right\}.$$

The mean pointwise coverage is the mean of the proportion of time grid values for which the confidence region contains the theoretical probability of default curve. It is given by

$${\displaystyle \frac{1}{N}}\sum _{j=1}^N\left({\displaystyle \frac{1}{n_T}}\sum _{k=1}^{n_T}I\left\{PD\left(t_k\right|x)\in \left(l(t_k,x),u(t_k,x)\right)\right\}\right).$$

Average width of the bootstrap confidence region is defined by

$${\displaystyle \frac{1}{N}}\sum _{j=1}^N\left({\displaystyle \frac{1}{n_T}}\sum _{k=1}^{n_T}\left(u(t_k,x)-l(t_k,x)\right)\right).$$

Winkler score (see [30]) is also used to compare the behavior of the methods. For classical confidence or prediction intervals, it is defined as the length of the interval plus a penalty if the theoretical value is outside the interval. Thus, it combines width and coverage. For values that fall within the interval, the Winkler score is simply the length of the interval. So low scores are associated with narrow intervals. When the theoretical value falls outside the interval, the penalty is proportional to how far the observation is from the interval. The formula of the Winkler score (WS) as a function of the time and covariate variables is as follows:

$$\begin{array}{ccc}\hfill \mathrm{WS}(t,x)& =& u(t,x)-l(t,x)+{\displaystyle \frac{2}{\alpha }}(l(t,x)-S\left(t\right|x))I\left(S\left(t\right|x)<l(t,x)\right)\hfill \\ \hfill & & +{\displaystyle \frac{2}{\alpha }}(S\left(t\right|x)-u(t,x))I\left(S\left(t\right|x)>u(t,x)\right).\hfill \end{array}$$

Since we are working with confidence regions for fixed $x\in I$ and t varying over the interval $I_T$, the integrated Winkle score is proposed as a criteria for the comparison of the confidence regions. It is defined by

$$IWS\left(x\right)={\int }_{I_T}WS(t,x)dt.$$

and the lower the value of IWS, the better the performance of the confidence region.

The results obtained are shown in

Table 5. Coverage, average width and IWS of the $95\%$ confidence regions by means of Beran’s and the smoothed Beran’s estimators using $N=300$ simulated samples from Models 1 and 2.

	Model 1
$P(\delta =0|X=x)$	0.2		0.5
Estimator	BERAN	SBERAN	BERAN	SBERAN
Coverage (%)	96.33	90.67	90.00	85.33
Mean pointwisecoverage (%)	99.94	98.05	99.63	96.85
Mean width	0.21997	0.09539	0.24827	0.10937
IWS	0.09869	0.04537	0.11218	0.05571
	Model 2
$P(\delta =0|X=x)$	0.2		0.5
Estimator	BERAN	SBERAN	BERAN	SBERAN
Coverage (%)	97.33	83.00	91.46	98.00
Mean pointwisecoverage (%)	99.88	98.53	99.65	99.85
Mean width	0.50514	0.17969	0.55581	0.33033
IWS	0.62590	0.22825	0.71009	0.40845

. The high values of pointwise coverage in all scenarios are remarkable. Furthermore, these coverage percentages are preserved when using double smoothing, while the average width of the confidence regions is halved. This is reflected in the IWS, which presents much larger values in the Beran’s estimator-based confidence regions.

This analysis is also illustrated in Figure 9 and Figure 10, where the confidence region for the probability of default of one sample from Models 1 and 2 is shown. These graphs show the higher variability of the Beran’s estimations in the resamples with respect to the smoothed Beran’s estimations. This leads to much wider confidence regions, especially at the right tail of the time distribution.

6. Application to Real Data

In this section, bandwidth selectors for Beran’s and the smoothed Beran’s estimators are applied to the German Credit dataset, and the confidence region of the probability of default is obtained. This dataset is publicly available on the webpage http://archive.ics.uci.edu/ml/datasets/Statlog+(German+Credit+Data) (accessed on 15 September 2021) and was previously analyzed in [31]. It includes information about 1000 credits, from which 293 were classified as bad credits and 707 as good credits. Then, the censoring ratio of this dataset is $70.7\%$. The duration of the credits in months is available along with the credit amount, checking account, savings amount and time of employment, among others.

The duration of the credits is set as the time to default, Z, the bad/defaulted credits are denoted by $\delta =1$ and the good credits by $\delta =0$. Let us denote the credit scoring with $X=(1,{\theta }_2,{\theta }_3,{\theta }_4)(X_1,X_2,X_3,X_4)$.

Since some of the original covariates are ordinal (interval) variables, we change them into numerical variables by following the criteria explained in [31]: $X_1$ is already a continuous variable denoting amount of credit in DM, $X_2\in \{-0.05,0.01,0.25,0\}$ denotes the amount of money in the checking account in thousands of DM, $X_3\in \{0,0.05,0.25,0.75,1.25\}$ denotes the savings amount in thousands of DM, and $X_4\in \{0,0.5,2.5,5.5,8.5\}$ denotes the years of employment. The single-index method proposed by [31] is used to estimate $(1,{\theta }_2,{\theta }_3,{\theta }_4)$. The credit scoring is obtained as follows:

$$X=X_1+3.2091X_2+0.2312X_3+2.1891X_4$$

Figure 11 shows the scatter plot between credit scoring and follow-up time variable by distinguishing between the censored and uncensored (and therefore, defaulted) credits. A dependency relationship between the two variables can be identified in the plot.

The probability of default, $PD\left(t\right|x)$, is estimated when $x=0.85$, which is a close value to the sample mean of the credit scoring, and $t\in [0,60]$. The bandwidth selector presented in Section 2.1 is used to approximate the optimal bandwidth for Beran’s estimator, obtaining $h^{*}=0.500$. The bandwidth selector presented in Section 2.2 gives the bootstrap approximation of the optimal bivariate bandwidth for the smoothed Beran’s estimator, $(h^{*},g^{*})=(0.102,13.614)$. The estimations of the conditional survival function and the probability of default by means of Beran’s and the smoothed Beran’s estimator with the corresponding bootstrap bandwidths are shown in Figure 12. The poor behavior of Beran’s PD estimator for large values of time is evident. The results obtained by the doubly smoothed Beran’s estimator seem to be more appropriate, since the roughness of the Beran’s estimaton is not expected in this type of curve. Supporting the conclusion of our real data analysis, the estimation of the probability of default over time for the assessment of risk in portfolios and bond rating obtained in [32,33] have shapes similar to those obtained here.

Finally, the confidence region methods proposed in Section 4 are applied. Since the MISE bandwidths are unknown in this context, bootstrap bandwidths are used. The bootstrap resamples and the resulting confidence regions at confidence level $95\%$ using each estimator are shown in Figure 13. The average width of the confidence region based on Beran’s estimator is $0.5581$, and the average width of the one based on the smoothed Beran’s estimator is $0.1438$. Note that the confidence region of $PD\left(t\right|x)$ is computed over the time interval $[0,40]$. Since the information provided in the right tail of the time distribution is sparse, Beran’s estimator performs very poorly, leading to extremely wide confidence regions. However, this problem is not as severe for the doubly smoothed Beran’s estimator, so the confidence region is computable for higher values of time. Figure 14 shows the confidence region of $PD\left(t\right|x)$ based on the smoothed Beran’s estimator with $t\in [0,60]$. The average width of this confidence region is $0.2398$.

In practice, the financial institution measures different features of its clients, such as age, amount of money in the bank account, salary, years of employment, etc. They summarize, usually by logistic regression, these covariates into the single variable credit scoring. Subsequently, techniques such as those shown in this paper allow the calculation of the probability of default at horizon b for all of them. The curve $PD\left(t\right|x)$ provides the probability that the client will default after a certain period of time b.

7. Conclusions and Future Lines

This article proposes an automatic bandwidth selector and confidence regions algorithms based on bootstrap selectors for Beran’s estimator and the smoothed Beran’s estimator of the probability of default in credit risk. The proposed resampling methods and bootstrap selectors allow to approximate the MISE bandwidths corresponding to each estimator: the covariate-smoothing bandwidth in the case of the Beran estimator and the two-dimensional covariate and time-smoothing bandwidth in the case of the doubly-smoothed estimator. In view of the simulation study carried out, it can be concluded that the bandwidth selectors work properly. The doubly smoothed Beran’s estimator with bootstrap bandwidths commits smaller estimation errors than Beran’s estimator. The simulation results also show the good behavior of the confidence regions, especially those based on the doubly smoothed Beran’s estimator. They have a lower average width, reducing the uncertainty about where the true probability of default curve is located, while preserving a high coverage.

The main limitation of proposed methods is their high computational cost. Approximating the bootstrap bandwidth or computing a confidence region by 500 resamples from a sample of size 100 requires one minute, while a sample of size 500 requires 25 min to obtain the result. These times are similar for both estimators. They seem to increase quadratically as the sample size grows, which may lead to prohibitive times for very large sample sizes. Using subsampling techniques is an appealing idea to be considered in the future for optimizing these methods.

In a financial context, credit scoring typically summarizes several interesting features of clients in order to measure their creditworthiness. However, this work could be extended to the case of having a multidimensional covariate $(X_1,\dots ,X_q)$, where each $X_i$ is a feature of the individual. Methods such as single-index can be useful for this purpose to avoid the curse of dimensionality. An approach along the lines similar to [31] can be used.