Phase transition in the Bayesian estimation of the default portfolio

TL;DR

This paper explores a phase transition phenomenon in Bayesian default probability estimation, showing how temporal correlation decay affects convergence and providing conditions for reliable estimation with empirical data.

Contribution

It introduces a hierarchical Bayesian method for PD estimation considering temporal correlations and characterizes the phase transition related to correlation decay.

Findings

Phase transition occurs at power index = 1 for correlation decay.

PD estimator convergence depends on the power index of decay.

Empirical data suggests long memory with a power index greater than one.

Abstract

The probability of default (PD) estimation is an important process for financial institutions. The difficulty of the estimation depends on the correlations between borrowers. In this paper, we introduce a hierarchical Bayesian estimation method using the beta binomial distribution and consider a multi-year case with a temporal correlation. A phase transition occurs when the temporal correlation decays by power decay. When the power index is less than one, the PD estimator does not converge. It is difficult to estimate the PD with limited historical data. Conversely, when the power index is greater than one, the convergence is the same as that of the binomial distribution. We provide a condition for the estimation of the PD and discuss the universality class of the phase transition. We investigate the empirical default data history of rating agencies and their Fourier transformations to confirm the form of the correlation decay. The power spectrum of the decay history seems to be 1/f, which corresponds to a long memory. But the estimated power index is much greater than one. If we collect adequate historical data,the parameters can be estimated correctly.

Tables

Table 1. Table 1: Asymptotic behavior of C ​ ( t ) 𝐶 𝑡 C(t) , and the scaling functions f τ ​ ( ξ t ) subscript 𝑓 𝜏 subscript 𝜉 𝑡 f_{\tau}(\xi_{t}) , f ξ ​ ( ξ t ) subscript 𝑓 𝜉 subscript 𝜉 𝑡 f_{\xi}(\xi_{t}) , and ξ t subscript 𝜉 𝑡 \xi_{t} . The assumed asymptotic form of C ​ ( t ) 𝐶 𝑡 C(t) is given in the second column. The second and the third columns provide the scaling functions. The last column contains the limit values of ξ ​ ( t ) / t 𝜉 𝑡 𝑡 \xi(t)/t .

No.	Asymptotic behavior	$f_{\tau }\left({\xi }_t\right)={lim}_{t\to \infty }\frac{\tau \left(st\right)}{\tau \left(t\right)}$	$f_{\xi }\left({\xi }_t\right)={lim}_{t\to \infty }\frac{\xi \left(st\right)}{\xi \left(t\right)}$	${\xi }_t={lim}_{t\to \infty }\xi \left(t\right)/t$
1	$C\left(t\right)\simeq c+\Delta C\left(t\right),c>0$	$s$	$s$	$1/\sqrt{3}$
2		$s^{1-\delta }=s^{\frac{2{\left(\xi /t\right)}^2}{1-{\left(\xi /t\right)}^2}}$	$s$	$\sqrt{\frac{1-\delta }{3-\delta }}$
3		$1$	$s^{\left(3-\delta \right)/2}$	$0$
4	$C\left(t\right)\propto t^{-\delta },\delta \ge 3$	$1$	$1$	$0$

Table 2. Table 2: MAP estimation of the parameters for the exponential and power decay models.

			Exponential decay			Power decay
No.	Model	$\theta$	${\rho }_D$	$r$	$\theta$	${\rho }_D$	$\gamma$
1	Moody’s 1920-2017	0.96%	1.9%	0.044	0.95%	2.0%	4.7
2	Moody’s 1920-2017 SG	2.37%	3.9%	0.044	2.35%	4.1%	4.7
3	Moody’s 1981- 2017	1.49%	0.7%	0.023	1.46%	0.7%	5.9
4	Moody’s 1990-2017	1.65%	0.7%	0.006	1.70%	0.8%	7.0
5	Moody’s 1981-2017 SG	4.25%	1.8%	0.020	4.29%	1.8%	6.0
6	S&P 1981-2017	1.54%	0.8%	0.024	1.54%	0.8%	5.7
7	S&P 1990-2017	1.72%	0.8%	0.006	1.72%	0.8%	7.5
8	S&P 1990-2017 SG	4.21%	2.0%	0.024	4.17%	1.9%	5.7

Equations

$$P(\theta ,{\rho }_D|X=k)=\frac{P(\theta ,{\rho }_D,X=k)}{P(X=k)}=\frac{P(X=k|\theta ,{\rho }_D)f(\theta ,{\rho }_D)}{P(X=k)},$$

$$P(\theta ,{\rho }_D|X=k)$$

$$\frac{\partial P(\theta ,{\rho }_D|X=k)}{\partial \theta }$$

$$\frac{\partial P(\theta ,z|X=k)}{\partial z}$$

$$P(\theta ,{\rho }_D|k_1,k_2)=\frac{P(k_2|\theta ,{\rho }_D,k_1)}{P(k_2)}\frac{P(k_1|\theta ,{\rho }_D)f(\theta ,{\rho }_D)}{P(k_1)}.$$

$$X(t+1)={\text{1}}_{U_{t+1}\le Z_d(t)},$$

$$Z_d(t)\equiv \frac{\alpha +{\sum }_{s=n_i}^{t}X(s)+{\sum }_{j=1}^i{d}_{i-j}{k}_j}{\alpha +\beta +(t-n_i)+{\sum }_{j=1}^i{d}_{i-j}{n}_j}.$$

$$V_{i+1}$$

$$V\sim \sum _{i=1}^{T}pq{n}_i+\sum _{i=1}^{T}pq{n}_i(n_i-1){\rho }_D+2pq{\rho }_D\sum _{i>j}^T{n}_i{n}_j{d}_{i-j}.$$

$$Corr\sim\rho_{D}\left(\begin{array}[]{ccccc}1&d_{1}&d_{2}&\cdots&d_{T}\\ d_{1}&1&d_{1}&\ddots&\vdots\\ \ddots&\ddots&\ddots&\ddots&\ddots\\ \vdots&\ddots&\ddots&\ddots&d_{1}\\ d_{T}&\cdots&d_{2}&d_{1}&1\\ \end{array}\right).$$

$$X(t+1)={\text{1}}_{U_{t+1}\le Z_d(t)}.$$

$$Z_d(t)\equiv \frac{\alpha +{\sum }_{s=1}^{t}X(s)d_{t-s}}{\alpha +\beta +{\sum }_{s=1}^t{d}_{t-s}}.$$

$$Z(t)\equiv \sum _{s=1}^{t}X(s)/t.$$

$$c_1(t)$$

$${\Delta }_t=2c_1(t)-t.$$

$${\Delta }_t$$

$$Y_{\tau }$$

$$P(y,\tau )$$

$$\text{d}{Y}_{\tau }=\frac{\alpha -\beta +{\int }_{\sigma =1}^{\tau }d(\tau -\sigma )\text{d}{Y}_{\sigma }}{\alpha +\beta +{\int }_1^{\tau }d(\tau -\sigma )\text{d}\sigma }\text{d}\tau +\sqrt{ϵ},$$

$$Y_{\infty }=\bar{v}\tau ,$$

$$\bar{v}=\frac{\alpha -\beta +\bar{v}\hat{T}}{\alpha +\beta +\hat{T}},$$

$$C(t)\equiv \text{mbox}E(X(t+1)|X(1)=1)-\text{mbox}E(X(t+1)|X(1)=0)=\frac{\text{mbox}Cov(X(1),X(t+1))}{\text{mbox}V(X(1))}.$$

$$\text{mbox}V(Z(t))=\text{mbox}{E}_{X(1)}(\text{mbox}V(Z(t)|X(1)))+\text{mbox}{E}_{X(1)}((\text{mbox}E(Z(t)|X(1))-\text{mbox}E(Z(t)){)}^2,$$

$$\text{mbox}{E}_{X(1)}((\text{mbox}E(Z(t)|X(1))-\text{mbox}E(Z(t)){)}^2)=\frac{1}{t^2}\frac{\alpha \beta }{(\alpha +\beta {)}^2}{\left(\sum _{s=0}^{t-1}C(s)\right)}^2.$$

$$\text{mbox}E(X(t+1)|X(1)=x)=\frac{\alpha +{\sum }_{s=1}^t\text{mbox}E(X(s)|X(1)=x)d_{t-s}}{\alpha +\beta +{\sum }_{s=1}^t{d}_{t-s}}.$$

$$C(t)=\frac{{\sum }_{s=1}^{t}C(s-1)d_{t-s}}{\alpha +\beta +{\sum }_{s=1}^t{d}_{t-s}}.$$

$$\sum _{s=1}^{t-1}C(s-1)r^{t-s}=r\sum _{s=1}^{t-1}C(s-1)r^{t-1-s}=r\cdot C(t-1)(\alpha +\beta +\sum _{s=1}^{t-1}{r}^{t-1-s}).$$

$$C(t)=\frac{1+r(\alpha +\beta +{\sum }_{s=1}^{t-1}{r}^{t-1-s})}{\alpha +\beta +{\sum }_{s=1}^t{r}^{t-s}}C(t-1).$$

$$r_{eff}\equiv \underset{t\to \infty }{\lim }C(t)/C(t-1)=r+\frac{1-r}{(\alpha +\beta )(1-r)+1}<1,$$

$$M_n(t)$$

Full text

I Introduction

Anomalous diffusion is an emerging subject in many fields Bro ; W2 ; G ; M . The models describing such phenomena depend on long memory. These are related to the phase transition, which has received considerable interest in sociophysics galam ; galam2 and econophysics Man . In previous papers, we investigated voting models that were similar to the Keynesian beauty contest Mori ; Hisakado2 ; Hisakado3 ; Hisakado4 ; Hisakado5 . This model has two kinds of phase transitions. One is the information cascade transition, which is similar to the phase transition of the Ising model Hisakado3 . The other is the convergence transition of super-normal diffusion Hod ; Hisakado2 .

Estimations of the probability of default (PD) and default correlation have been obtained from empirical studies on the historical data from credit events. These two parameters are important for pricing financial products such as synthetic CDOs M2010 ; M2008 ; Sch . Also called “long run PDs”, these parameters are important to financial institutions for portfolio management. If the number of defaults is minimal, it is not easy to estimate these parameters Tas ; FSA .

In this paper, we introduce a Bayesian estimation method using the beta-binomial distribution Hisakado ; Wit . For the usual cases the Merton model, which incorporates the default correlation by the correlation of the asset price movements (asset correlation), is used to estimate the PD and default correlation Mer . Monte Carlo simulations are necessary to estimate these parameters, except in the limit of large homogeneous portfolios, where the Merton model is used Sch . In the beta-binomial case, default correlation, instead of asset correlation, is used Hisakado . Moreover, we consider a multi-year case with temporal correlation, which refers to a time-dependent correlationTas ; FSA .

A phase transition occurs when the temporal correlation decays by power-law. A power-law decay implies that the PD has a long memory compared to that of exponential decay Man . When the power index is less than one, the estimator distribution of the PD does not converge to the delta function. Alternatively, when the power index is greater than one, the convergence is the same as that of the normal case. When the distribution does not converge, it is difficult to estimate the PD with limited data. The required condition for estimating the PD is clarified. The critical exponents for the power-law decay of the correlation function depend on microscopic features of the model. The universality class of the phase transition is different from those of the nonlinear Pólya urns Mori5 ; Mori6 .

To confirm the decay form of the temporal correlation, we investigate the empirical default data history using Fourier transformations. We determine whether the power spectrum of the default history follows 1/f Kes ; Man . When this condition is satisfied, it corresponds to the correlation of the PD with long memory where a phase transition of the convergence exists. However, it is difficult to accurately confirm a 1/f power spectrum when the estimation of the power index is much greater than one. It follows that when there is adequate historical data, parameters such as PD, default correlation, and temporal correlation can be estimated correctly.

The remainder of this paper is organized as follows. In section 2, we introduce a hierarchical Bayesian estimation method using the beta-binomial distribution. In section 3, we consider the convergence of the PD estimator. In section 4, we study the phase transition of the Pólya urn with a discount factor using an analytic method and a finite-size scaling analysis. In section 5, we apply the Bayesian estimation to the empirical data of default history. Finally, the conclusions are presented in section 6.

II Bayesian Estimation using Beta-binomial distribution

We denote the PD estimation as $\theta$ and default correlation as $\rho_{D}$, where $0\leq\theta\leq 1$ and $0\leq\rho_{D}\leq 1$. The distribution of $\theta$ and $\rho_{D}$ is $P(\theta,\rho_{D})$. The number of obligors in the portfolio is $n$. $\theta$ and $\rho_{D}$ are estimated using a Bayesian estimation. We consider the Bernoulli random variables $X_{i}(i=1,2,\cdots,n)$ that take the values $1$ or [math] . When the obligor, $i$, is the default (non-default), $X_{i}=1(0)$. We define $X=\sum_{j=1}^{n}X_{j}$ and consider a default correlation for $X_{i}$, and not an asset correlation.

When the number of defaults is $k$, the Bayes formula for the posterior distribution $P(\theta,\rho_{D}|X=k)$ is

[TABLE]

where $f(\theta,\rho_{D})$ is the prior distribution.

We use the beta-binomial distribution for $P(X=k|\theta,\rho_{D})$. The posterior distribution is given by

[TABLE]

where, $\theta=\frac{\alpha}{\alpha+\beta}$ and $\rho_{D}=\frac{1}{\alpha+\beta+1}$. Hence, we obtain the relations $\alpha=\theta\frac{1-\rho_{D}}{\rho_{D}}$ and $\beta=(1-\theta)\frac{1-\rho_{D}}{\rho_{D}}$. Here, we use the beta function $B(\alpha,\beta)=\Gamma(\alpha)\Gamma(\beta)/\Gamma(\alpha+\beta)$.

We consider the maximum a posteriori (MAP) estimation of Eq. (2). When the prior function $f(\theta,\rho_{D})$ is a constant function, the maximum point is

[TABLE]

where $\varphi(x)$ is the digamma function. The summation from $i=1$ to $k$ is a monotonously decreasing function of $\theta$ because $\alpha$ increases, while the second summation in Eq. 3 is a monotonously increasing function about $\theta$ because $\beta$ decreases. When $\theta\sim 0$, the difference of the two summations is positive. Conversely, when $\theta\sim 1$, the difference of the two summations becomes negative. Hence, the function $P(\theta|X=k,\rho_{D})$ has one peak in the range $0<\theta<1$. The multi-term case is provided in Appendix A.

Next, we consider the variable $\rho_{D}$. The maximum point is

[TABLE]

where $z=(1-\rho_{D})/\rho_{D}$.

All the summations in the last term with parenthesis are monotonously decreasing functions about $z$. When $z\sim 0$, the last term becomes positive. Conversely, when $z>>1$, the last term becomes 0. When $(k-1)/n\leq\theta\leq k/n$ or is adequately close to this condition, the last term becomes positive. $(k-1)/\theta\leq n-1$ and $(n-k-1)/(1-\theta)\leq n-1$ become $(k-1)/(n-1)\leq\theta\leq k/(n-1)$. In this case, the last term increases monotonously, and the peak is $z=\infty$ and $\rho_{D}=0$. This implies that the optimization of $\rho_{D}$ is zero for the single term model. When $\theta$ is not adequately close to $(k-1)/n\leq\theta\leq k/n$, the last term changes from positive to negative as $z$ increases. Therefore, one peak occurs in $P(\theta,z|X=k)$.

We extend this method to the multi-year case. There are $n_{i}$ obligors in year $i$ and $k_{i}$ defaults occur. The prior distribution for the second year is the posterior distribution, which is calculated from the first year’s data. In this way, the posterior distribution is updated every year. We write the posterior distribution $P(\theta,\rho_{D}|k_{1},k_{2})$ as

[TABLE]

It is natural to assume that the number of defaults of the current year is affected by the number of defaults in previous years, thus the defaults have a temporal correlation. When the default rate is high (low), it is reasonable to assume that the default rate will be high (low) in the next year. This is similar to volatility clustering, which has a long memory AR ; GA as well. We confirm this using empirical data in the following sections.

We introduce the temporal correlation by adjusting $\alpha$ and $\beta$, and consider the $j$ th year. The number of obligors and defaults in the $j$ th year are $n_{j}$ and $k_{j}$. In the same year, the correlation is $\rho_{D}$. We set the temporal correlation parameters between the $i$ and $j$ th years; $d_{i-j}$ and $j<i$. $\alpha$ and $\beta$ are adjusted to $\alpha+\sum_{j=1}^{i-1}d_{i-j}k_{j}$ and $\beta+\sum_{j=1}^{i-1}d_{i-j}(n_{j}-k_{j})$ Mori18 . This implies that the previous years’ data affects the present defaults. It is easy to confirm that $d_{i}=1$ indicates that all the data is correlated to $\rho_{D}$. When $d_{i}=0$, the data is independent each year.

III Correlation Decay

In the previous section, $d_{i}$ was introduced to represent the temporal correlation. In this section, to clarify the behavior of the parameter $d_{i}$, where $i=1,2,\cdots,T$ and $d_{0}=1$, the variance of the stochastic process is considered. In each year the diffusion has $n_{i}$ steps and $k_{i}$ defaults, where $i=1,2,\cdots,T$.

The adjustments related to parameters $\alpha$ and $\beta$ are the effects of the temporal correlation from the previous conclusions. We shrink the previous years’ conclusions and add them to the initial parameters for the adjustment process. The shrinking ratio for the interval $i$ is $d_{i}$

The two term model is examined first. We consider the relation between the first and second years. $n_{1}$ and $k_{1}$ are the number of obligors and defaults, respectively, in the first year. The second year’s parameters become $\alpha+d_{1}k_{1}$ and $\beta+d_{1}(n_{1}-k_{1})$. We consider the shrinking processes from $\alpha$ to $\alpha+d_{1}k_{1}$ and $\beta$ to $\beta+d_{1}(n_{1}-k_{1})$. The variance of the second term of process is $n_{1}d_{1}pq+d_{1}n_{1}(n_{1}-1)pq\rho_{D}$, where $q=1-p$; that is, we approximate $d_{1}B_{\alpha,\beta}(k_{1},n_{1}-k_{1})\sim B_{\alpha,\beta}(d_{1}k_{1},d_{1}n_{1}-d_{1}k_{1})$ where $B_{\alpha,\beta}$ is the beta-binomial distribution with parameters $\alpha$ and $\beta$. We approximate this variance by $n_{1}d_{1}pq+d_{1}n_{1}(d_{1}n_{1}-1)pq\rho_{D}$, and the difference becomes $n_{1}^{2}pq\rho_{D}d_{1}(1-d_{1})\geq 0$. Hence, the approximation is exact when $d_{1}=0,1$ or $\rho_{D}=0$. However, if $d_{1}\sim 0,1$ or $\rho_{D}\sim 0$, this approximation can be used. In other cases, the real variance is larger than the approximation. We use this approximation to study the meaning of this process.

For the defaults of the obligors, the hypothesis $d_{1}\sim 0$ or $1$ and $\rho_{D}\sim 0$ is given. In other words, the temporal correlation is either a high or low case, or a low correlation case. Hereafter, we use this approximation to calculate the variance of this process.

We extend the stochastic process to the multi-year case. Let $\{U_{t};t\geq 1\}$ be an independent and identically distributed (i.i.d.) sequence that is uniformly distributed on [0,1]. The discrete dynamics of the process is described by:

[TABLE]

when $n_{i}+1\leq t\leq n_{i+1}$. Here $Z_{d}(t)$ is given by

[TABLE]

The expectation value of $X(t)$ is $\mbox{E}(X(t))=\alpha/(\alpha+\beta).$ When $d_{i}=1$, the process is beta-binomial.

We consider the relationship between the year $i$ and $i+1$. The distribution of year $i$ is a beta-binomial distribution. Hence, the conditional variance, $V_{i+1}$, of the year $i+1$ can be evaluated, using the above approximation, as

[TABLE]

Therefore, the difference of the summations $\sum_{j=1}^{i+1}n_{j}d_{i+1-j}$ and $\sum_{j=1}^{i}n_{j}d_{i+1-j}$ correspond to the the variance of $(i+1)^{th}$ step. Therefore, using this approximation, the correlation between the $i^{th}$ and $j^{th}$ years is approximated by $\rho_{D}d_{i-j}$. The term $d_{i-j}$ plays the role of a discount factor in the correlation $\rho_{D}$. It can be seen that as time progresses, the correlation is discounted. It is reasonable to assume a monotonically decreasing function for $d_{i}$ because the effects decrease as the distance between $i$ and $j$ increase.

The total variance for the diffusion is approximated by

[TABLE]

The first, second, and third terms correspond to the variance for binomial distribution, constant correlation $\rho_{D}$ in the portfolio, and temporal correlation, respectively.

In summary, when $d_{i}\sim 0,1$ or $\rho_{D}\sim 0$, the correlation between year $i$ and year $j$ is approximated by

[TABLE]

The average PD, correlation of the Bernoulli random variables, and temporal correlation using this approximation are $p$, $\rho_{D}$, and $d_{i}$, respectively.

In the Bayesian estimation, if the scaled variance converges as the data increases, these parameters can be estimated correctly. Conversely, if the variance does not converge, the parameters cannot be estimated. It should also be considered whether the process has a stationary solution, which will be discussed regarding the spectrum analysis in the following sections.

It is difficult to estimate all the $d_{i}$ values due to limited data. By introducing a prior distribution for $d_{i}$, the estimation becomes a hierarchical Bayesian estimation. It is reasonable to assume that the prior distribution is a monotonically decreasing function. Therefore, we considered two hyperprior distributions, an exponential and power decay, to have long memory.

IV Phase transition in the estimation of PD

In this section we determine whether the PD in the Bayesian estimation converges. To simplify the model, we set $n_{j}=1,j\geq 1$ in Eq. (7). This does not affect the outcome of the PD estimation. Let $\{U_{t};t\geq 1\}$ be an independent and identically distributed (i.i.d.) sequence that is uniformly distributed on [0,1]. The discrete dynamics of the process is described by:

[TABLE]

Here, $Z_{d}(t)$ is the weighted sum of $X(s),s\leq t$ with the discount factor $d_{t-s}$,

[TABLE]

This is the Pólya urn modelPolya with a discount factor $\{d_{i}\}$.

The expectation value of $X(t)$ is $\mbox{E}(X(t))=\alpha/(\alpha+\beta).$ The PD estimator is $Z(t)$,

[TABLE]

The success of the PD estimation depends on the the behavior of the variance of $Z(t)$. More specifically, if the variance of $Z(t)$ converges, then the PD can be estimated.

IV.1 Stochastic differential equation

First, the stochastic process is rewritten using $c_{1}(t)=\sum_{s=1}^{t}X(s)$;

[TABLE]

where $P_{k,t}$ and $Q_{k,t}$ are the process probabilities. The sum of $P_{k,t}$ and $Q_{k,t}$ is $1$.

For convenience, we define a new variable $\Delta_{t}$ such that

[TABLE]

We change the variables from $k$ to $\Delta_{t}$ and $X(s)$ to $x_{s}=2X(s)-1$. Given $\Delta_{t}=u$, we obtain a random walk model:

[TABLE]

We now consider the continuous limit $\epsilon\rightarrow 0$,

[TABLE]

where $\tau=t/\epsilon$ and $y=\Delta_{t}/\epsilon$. On approaching the continuous limit, we obtain the following stochastic partial differential equation:

[TABLE]

where $d(\tau)$ is the continuous function of $d_{t}$, the discount factor, and $dY_{\tau}=\epsilon x_{[t/\epsilon]}$.

We are interested in the behavior of $Y_{\tau}$ in the limit $\tau\rightarrow\infty$. We assume that the stationary solution is

[TABLE]

where $\bar{v}$ is a constant. Substituting Eq. (15)into Eq. (14), we obtain

[TABLE]

where $\hat{T}=\lim_{\tau\rightarrow\infty}\int_{1}^{\tau}d(\tau-\sigma)\textrm{d}\sigma$.

Eq. (16) is a self-consistent equation. When $\hat{T}<\infty$, Eq. (16) is solved when $\bar{v}=(\alpha-\beta)/(\alpha+\beta)$. The process converges to the average point. On the other hand, when $\hat{T}\rightarrow\infty$, we can obtain the identity equation $\bar{v}=\bar{v}$, suggesting that the process does not converge to the delta function. The expected value of $Y_{s}$ is $(\alpha-\beta)/(\alpha+\beta)$. Hence, the phase transition at the point $\hat{T}$ diverges to infinity. When the distribution does not converge, we cannot estimate the parameters correctly, even if the amount of data increases. This is a critical issue when using the Bayesian estimation. In other words, $\hat{T}<\infty$ is a compulsory condition for parameter estimation.

IV.2 Correlation function and finite size scaling analysis

To understand the phase transition, we investigated the correlation function, $C(t)$. $C(t)$ is defined as the correlation between $X(1)$ and $X(t)$ such that

[TABLE]

The function $C(t)$ represents the propagation of the memory of $X(1)$ to later variables $X(t+1)$. To understand the relationship between the variances of $Z(t)$ and $C(t)$, the variance of $Z(t)$ can be written as

[TABLE]

where $\mbox{V}(Z(t)|X(1))$ is the conditional variance of $Z(t)$ on $X(1)$. The expectation value of $x$ is $\mbox{E}_{X(1)}(x)$ and the probability function is $P(X(1))$. The second term on the right-hand side of Eq. (18) represents the variance of $\mbox{E}(Z(t)|X(1))$ from the dependence on $X(1)$. In Eq. (18), the second term is related to $C(t)$ as it originates from the dependence of $\mbox{E}(Z(t)|X(1))$ on $X(1)$. We write the second term of $C(t)$ as

[TABLE]

If $c=\lim_{t\to\infty}C(t)>0$, $\lim_{t\to\infty}\mbox{V}(Z(t))>0$ and $Z(t)$ does not converge.

Using $X(t+1)=\textbf{1}_{U_{t+1}\leq Z_{d}(t)}$, we obtain the next relation for the conditional expectation value of $X(t+1)$ with the condition $X(1)=x$, $\mbox{E}(X(t+1)|X(1)=x)$, as

[TABLE]

As $C(t)=\mbox{E}(X(t+1)|1)-\mbox{E}(X(t+1)|0)$, we obtain the following recursive relation for $C(t)$ as

[TABLE]

This recursive relation contains all information regarding the asymptotic behavior of $C(t)$. If one assumes a functional form for $d_{i}$ with the initial condition $C(0)=1$, we can estimate $C(t)$ for $t\geq 1$.

IV.2.1 Exponential decay case

We consider the exponential decay case, $d_{i}=r^{i},r\leq 1$. $\hat{T}$ is finite and there is no phase transition. We decompose the numerator of Eq. (19) as $C(t-1)+\sum_{s=1}^{t-1}C(s-1)r^{t-s}$. We rewrite the second term using Eq. (19) for $t-1$ as

[TABLE]

We then obtain the next recursive relation for $C(t)$:

[TABLE]

As we are interested in the asymptotic behavior of $C(t)$, we estimate the decay rate, $r_{eff}$, with $C(t)\sim r_{eff}^{t}$, which gives

[TABLE]

where $r_{eff}<1$ for $r<1$, and $C(t)$ decays exponentially.

Numerical studies of the system were performed. To estimate $C(t)$, the recursive relation of Eq. (19) is solved for $t\leq 2\times 10^{5}$. A Monte Carlo sampling procedure is adopted for the variance of $Z(t)$. We obtained $10^{4}$ sample sequences for $\{X(t)\},t=1,\cdots,2\times 10^{5}$ and estimated the variance of $Z(t)$. Figure 1 (a) shows the plot of $C(t)$ vs. $t$. It is clearly shown that $C(t)$ decays exponentially. Figure 1 (b) shows the plot of V$(Z(t))$ vs. $t$. For all $r<1\in\{0.8,0.9,0.99\}$, $V(Z(t))$ decays as $1/t$. When $r=1$, the $Z(t)$ distribution converges to the beta distribution. Hence, there is no phase transition for $r<1$.

IV.2.2 Power-law decay case

For the case of power-law decay, namely $d_{i}=\frac{1}{(1+i)^{\gamma}}$, when $\gamma>1$ and $\hat{T}<\infty$, the process converges to the delta function. On the other hand, when $\gamma\leq 1$ and $\hat{T}$ goes to infinity, the process does not converge.

The behaviors of $C(t)$ and $V(Z(t))$ were investigated by the numerical method, in the same manner as the exponential decay case. Fig. 2 (a) shows the double logarithmic plot of $C(t)$ vs. $t$. It can be seen that $C(t)$ decays with a power-law form for $\gamma\in\{1.5,2,3\}$. For small $\gamma$, such as $\gamma=0.5,0.1$, the slope is extremely small. Fig. 2 (b) shows the double logarithmic plot of $V(Z(t))$ vs. $t$. For $r=3.0,2.0,1.5$, $V(Z(t))$ decays as $1/t$. At $\gamma=1$, the slope of the decay is less than one. For $r<1$, the curve is concave down. These results suggest the validity of the self-consistent equation analysis.

To investigate the phase transition, we apply finite-size scaling (FSS) analysis Mori5 . We define the relaxation and second-moment correlation times, $\tau(t)$ and $\xi(t)$, respectively, using the $n^{th}$ moment of $C(t)$ as

[TABLE]

For FSS, we assume that the scaling function, $\lim_{t\to\infty}A(st)/A(t)$, for some observable, $A(t)$, with a scale factor, $s$, is expressed as a function of $\xi_{t}\equiv\lim_{t\to\infty}\xi(t)/t$ such that

[TABLE]

We assume the following asymptotic forms for $C(t)$;

[TABLE]

Here, $c=\lim_{t\to\infty}C(t)$ is the order parameter of the phase transition and $c^{\prime}$ is a constant. Using the asymptotic forms, we can classify the behavior of the scaling functions. We show the results for $f_{\tau}(\xi_{t}),f_{\xi}(\xi_{t})$ and $\xi_{t}$ in Table I. (In detail, see Appendix B)

Figure 3 shows the numerical estimations of $\xi(2t)/\xi(t)$ and $\tau(2t)/\tau(t)$ vs. $\xi(t)/t$ with $t=10^{5}$. The symbols show the fixed points under the renormalization transformation $t\to 2t$. There are two stable fixed points at $\xi_{t}=0$ and $\xi_{t}=1/\sqrt{3}$, and one unstable fixed point at $\xi_{t}=\sqrt{(1-\delta)/(3-\delta)}\simeq 0.4073\equiv\xi_{t}^{c}$. If $\xi_{t}>\xi_{t}^{c}$, then $\xi(2t)/\xi(t)>2$ and $\xi(t)/t$ moves to $1/\sqrt{3}$ under the transformation $t\to 2^{n}t$ and $n\to\infty$. $\xi(t)$ diverges linearly with the system size, $t$, at the fixed point, which reflects the memory of $X(1)$ that retains. If $\xi_{t}<\xi_{t}^{c}$, $\xi(2t)/\xi(t)<2$ and $\xi(t)/t$ moves to [math]. $\lim_{t\to\infty}\xi(t)<\infty$ and the memory of $X(1)$ is lost for sufficiently large $t$. At the stable fixed points of $\xi_{t}=1/\sqrt{3}$ and at $\xi_{t}=0$, $\tau(2t)/\tau(t)$ becomes 2 and 1, respectively. From the unstable fixed point at $\xi_{t}=\xi_{t}^{c}$, we can estimate $\delta$ using $f_{\tau}(\xi_{t}^{c})=2^{1-\delta}\simeq 1.3174$. This estimation is in accordance with the estimation from $\xi_{t}^{c}=\sqrt{(1-\delta)/(3-\delta)}\simeq 0.4073$. These results support the phase transition between the two phases, $C(t)\simeq c+\Delta C(t),c>0$ and $C(t)\propto t^{-\delta},\delta>1$, in the limit $t\rightarrow\infty$. At the critical point $\gamma=1$, $\xi_{t}=\xi_{t}^{c}$ and $C(t)\propto t^{-\delta}$ with $0<\delta<1$.

We estimate $c$ by $C(2t)$ and $\delta$ by $\log_{2}C(t)/C(2t)$ with $t=10^{3}$ and $10^{5}$. By comparing the values for $t=10^{3}$ and $10^{5}$, one can anticipate the limit behavior $t\to\infty$. The results are shown in Fig. 4. Figure 4 (a) shows $C(2t)$ vs. $\gamma$. For $\gamma>1$, $C(2t)$ is almost zero. For $\gamma<1$, $C(2t)$ is positive. The derivative of $c$ at $\gamma=1$ is seemingly continuous. Figure 4 (b) shows $\delta$ vs $\gamma$ with $(\alpha,\beta)=(1,1)$ and $(1,4)$. For $\gamma>1$, one can anticipate that $\delta=\gamma$ by observing the change from $t=10^{3}$ to $10^{5}$. For $\gamma<1$, $\delta=0$ which suggests that $c>0$. At the critical point $\gamma=1$, $\delta$ depends on $(\alpha,\beta)$.

Next, we investigated $\delta$ at the critical point $\gamma=1$. We assume that $C(t)\propto t^{-\delta}$. Eq. (19) can be approximated in the continuous limit as

[TABLE]

By the following change of variables, $(t+1)\mu=s$, we obtain

[TABLE]

We see that $\delta$ depends on $\alpha$ and $\beta$ through the combination $\alpha+\beta$. In the limit $t\rightarrow\infty$, when $\delta=1$ and [math], $\alpha+\beta=0$ and $\alpha+\beta\rightarrow\infty$, respectively. The critical exponent $\delta$ is in the range $\delta\in(0,1)$. Figure 4 (c) shows $\delta$ vs. $\alpha+\beta$ for $\gamma=1$. We adopt two cases $\alpha:\beta=1:1$ and $1:4$. The symbols show the results of the numerical estimation, and the solid line shows the results by numerically solving Eq. (24). The results for $\alpha:\beta=1:1$ and $1:4$ collapses onto the same curve vs. $\alpha+\beta$, which confirms that $\delta$ depends on $\alpha$ and $\beta$ through $\alpha+\beta$.

V Is the temporal correlation decay exponential or power?

In this section, we use three data sets from the default data. Two sets are rating agency data, and the other is from a Japanese company.

V.1 Standard & Poor’s data

As discussed in the previous section, temporal correlation is a critical issue for determining whether there is an exponential or a power decay. This affects whether the parameters are estimated correctly. In this section we investigate the temporal correlation using empirical data. First, the S&P default data from 1981 to 2017 Data1 are used. The average PD is 1.58 $\%$ for all ratings and 3.09 $\%$ for speculative ratings. A speculative grade rating represents the rating under BBB-(Baa3). In Fig. 5 (a) we show the historical default rate. The solid and dotted lines correspond to all the samples and the speculative grade, respectively, below BBB+(Baa3).

The autocorrelation is shown in Fig. 6 (a). The x-axis represents the year. The exponential decay and cyclical increase are confirmed. This represents the cyclical bubbles and their collapse in recent years. However, it is difficult to confirm whether the decay is exponential or power-law from the autocorrelation data alone. Therefore, a Fourier transformation was applied to the PD data in Fig. 7 (a), but it was still difficult to obtain confirmation because the data was annual, and its size was not very large.

V.2 Moody’s data

Next, we used Moody’s default data from 1920 to 2017 for 98 years Data2 . It includes the Great Depression in 1929 and Great Recession in 2008. It is one of the longest sets of default data Hec . The average default rate is 1.56$\%$ for all the ratings and $3.87\%$ for the speculative ratings. In Fig. 5 (b), we show the historical default rate.

The autocorrelation is shown in Fig. 6 (b). The x-axis represents the year. The exponential decay is confirmed for a short time. Over the long historical data, we cannot confirm the cyclical trend that was observed in recent years. We applied a Fourier transformation to the default ratio data in Fig. 7 (b). as it is difficult to confirm whether the decay is exponential or power-law from the autocorrelation alone.

V.3 Risk Data Bank data

Next, we apply our data to the risk data bank (RDB) data Data3 . The data covers all of the enterprise data without individual owner-managers in Japan. The data is monthly from 2001 to 2017 and the seasonal effects were adjusted. The historical data and autocorrelation are shown in Fig. 8 (a), which is different from the previous two samples. The slow decay of the correlation was confirmed. In Fig. 8 (b) 1/f fluctuations were confirmed. This corresponds to the power decay of the Wiener-Khinchin theorem, which shows the relationship between the autocorrelation and power spectrum by a Fourier transformation. In Fig. 9 we show the power spectrum for each sector, namely, construction, wholesale, real estate, retail sales, other services, and manufacturing. The solid line represents the trend. We can conclude that the temporal correlation may contain a long memory for this data. However, it is difficult to confirm a strict power law.

VI Estimation of parameters

We estimate the long run probability of default, $\theta$, and the default correlation, $\rho_{D}$, for S&P and Moody’s data by the MAP estimation. We use a uniform distribution for the prior distribution $f(\theta,\rho_{D})$. As discussed in the previous section, the exponential and power decays are used for the temporal correlation. The conclusions are listed in Table 2 for the exponential and power decay models. We confirmed a small $r$ value that represents the small temporal correlation. The parameter $\gamma$ for the power decay is greater than the phase transition point, $\gamma=1$. The PD and default correlation are almost the same as the estimations by the exponential and power decay models. The reason is that the power index, $\gamma$, is adequately large and there is only a small difference between the exponential and power decay models. The first- and second-year’s temporal correlations, $d_{1}$ and $d_{2}$, respectively, are important for representing the data.

The parameters depend on the data terms. In the recent past, the default and temporal correlations have become minimal. This may depend on the smooth financial operations of governments and central banks. Alternatively, the long history data of  100 years have long correlations that are less than the phase transitions. This depends on the old data before the 1980s. For the RDB data, we can estimate $\gamma=2$, which is in the normal convergence phase. Hence, we can estimate the PD by the Bayesian formula, which we introduced.

VII Concluding Remarks

In this paper, we introduced a hierarchical Bayesian estimation method using the beta-binomial distribution to estimate the parameters, probability of default (PD), and default correlation. Moreover, we considered a multi-year case with temporal correlation. We confirmed phase transitions when the temporal correlation decayed by a power curve, which means that the correlation had a long memory. Conversely, for the case of exponential decay, there was no phase transition. When the power index, $\gamma$, was above or equal to one, the estimator distribution of the PD converged. Conversely, when the power index was less than 1, the distribution did not converge. The critical exponent $0<\delta<1$ depended on the microscopic feature of the model and the universality class of the phase transition differed from those of the nonlinear Pólya urn. We call this phase transition a ”short memory-long memory transition”. In summary, the condition for the estimation of parameters is $\hat{T}=\lim_{\tau\rightarrow\infty}\int_{1}^{\tau}d(\tau-s)\textrm{d}s<\infty$.

To confirm the form of the decay, we investigated the empirical default history data using a Fourier transformation. We determined that the power spectrum of the default history was seemingly 1/f, which implies that the correlation had a long memory for the RDB monthly data. We applied this method to the historical data and estimated the parameters. The region of the power index provided normal convergence. We have demonstrated that, for adequate data collection, these parameters can be estimated correctly.

Acknowledgements.

This work is supported by JPSJ KAKENHI [Grant No. 17K00347].

Appendix A MAP estimation for Multi-year case

We extend the maximum a posteriori (MAP) estimation, which we discussed in section 2 for the multi-year case. The number of obligors and defaults in the $j^{th}$ year are $n_{j}$ and $k_{j}$, respectively. When a prior function, $f(\theta,\rho_{D})$, is a constant function, the maximum point is

[TABLE]

where $\varphi(x)$ is the digamma function. $\alpha_{j}$ and $\beta_{j}$ are the adjusted $\alpha$ and $\beta$. $\alpha_{j}=\alpha+\sum_{l=1}^{j-1}d_{j-l}k_{l}$ and $\beta_{j}=\beta+\sum_{l=1}^{j-1}d_{j-l}(n_{l}-k_{l})$. The first term in the last set of parentheses in Eq. (25) is a monotonously decreasing function about $\theta$, because $\alpha$ increases. The second term in the last set of parentheses is a monotonously increasing function about $\theta$, because $\beta$ decreases. When $\theta\sim 0$, the difference of the two terms is positive because $\alpha_{1}=\alpha$. In contrast, when $\theta\sim 1$, the difference of the two terms becomes negative because $\beta_{1}=\beta$. Hence, the function, $P(\theta|X=k,\rho_{D})$, has one peak in the range $0<\theta<1$.

Appendix B Scaling functions $f_{\xi}(\xi_{t})$ and $f_{\tau}(\xi_{t})$

We define the relaxation and second-moment correlation times, $\tau(t)$ and $\xi(t)$, respectively, using the $n^{th}$ moment of $C(t)$ as in Eq. (LABEL:moment). If we assume that $C(t)\propto t^{-\delta}$, $M_{n}(t)$ behaves as

[TABLE]

Using the asymptotic behavior of $M_{n}(t)$, we find $\tau(t)$ behaves as

[TABLE]

$\xi(t)$ behaves as

[TABLE]

The scaling function for $\tau$ is defined as $f_{\tau}(\xi_{t})\equiv\lim_{t\to\infty}\frac{\tau(st)}{\tau(t)},s>1$. From the asymptotic behavior of $\tau(t)$, we have

[TABLE]

For $\delta<1$, $\xi_{t}\equiv\lim_{t\to\infty}\xi(t)/t=\lim\sqrt{(1-\delta)/(3-\delta)}$ and the scaling function is given in terms of $\xi_{t}$ as

[TABLE]

$\xi_{t}=1/\sqrt{3}$ and $f_{\tau}(\xi_{t})=2$ in the limit $\delta\to 0$.

The scaling function for $\xi$ is defined as $f_{\xi}(\xi_{t})\equiv\lim_{t\to\infty}\frac{\xi(st)}{\xi(t)}$. We have

[TABLE]

By the renormalization transformation $t\to s^{n}t$, $\lim_{n\to\infty}\xi(s^{n}t)/s^{n}=\xi(t)$ for $\delta\leq 1$. For $\delta>1$, $\xi(s^{n}t)/s^{n}=0$. The critical state of the system exists at $\delta<1$.

We assume $C(t)\simeq c+\Delta C(t),c>0$ and $\Delta C(t)$ rapidly decays to zero. $\lim_{t\to\infty}\tau(t)=ct$ and $\xi_{t}=1/\sqrt{3}$. $f_{\xi}(\xi_{t})\equiv\lim_{t\to\infty}\xi(st)/\xi(t)=s$ and $f_{\tau}(\xi_{t})\equiv\lim_{t\to\infty}\tau(st)/\tau(t)=s$ holds.