Climate-Related Default Probabilities

Abstract

Climate risk refers to the risks associated with climate change and has already started to impact various sectors of the economy. In this work, we focus on the impact of physical risk on the probability of default for a firm in the agribusiness sector. The probability of default is estimated based on the Merton model, where the firm defaults when its asset value falls below the threshold defined by its liabilities. We study the relationship between the stock value of the firm and global surface temperature anomalies, observing that an increase in temperature negatively affects the stock value and, consequently, the asset value of the firm. A decrease in the asset value of the firm translates into an increase in its probability of default. We also propose a model to assess the exposure of the firm to transition risk.

1. Introduction

Climate risk refers to the risks associated with climate change and has already begun to impact various sectors of the economy, including agriculture, fisheries, and tourism, to name just a few. There are two types of risks linked to climate change: physical risk, which is related to extreme weather events, and transition risk, which arises from the shift from a carbon-intensive economy to a low-carbon economy.

Given the intricacy of the subject, an extensive body of research has been developed in the last two decades on the relationship between climate change exposure on the one side and firm credit risk, pricing of financial assets, and investment portfolio allocation on the other side (see Gianfrate 2018; Gianfrate and Peri 2019; Oikonomou et al. 2014).

Bauer and Hann (2010) analyzed environmental management and its implications for bond investors. It has been proven that poor environmental practices influence the credit merit of borrowing firms through the legal, reputational, and regulatory risks associated with environmental incidents. In contrast, firms with proactive environmental engagement benefit from a lower cost of debt financing. Using a panel least squares regression, Capasso et al. (2020) investigated the relationship between climate-related risks and firm creditworthiness, measured by the distance to default. Companies with larger carbon footprints are relatively more exposed to progressively stricter climate-related regulations. Consequently, their future cash flows and, hence, firm asset values, are likely to be influenced to a larger extent than those of companies with smaller carbon footprints. Therefore, companies with high carbon footprints are perceived by the market as more likely to default. These findings clearly indicate that financial markets are already factoring in the climate change exposure of listed companies and that exposure to climate risks affects the creditworthiness of loans and bonds issued by corporations. Ilhan et al. (2020) estimated the effects of carbon emissions on downside risk as reflected in the options market. The authors found that higher carbon emissions increase downside risk and that the cost of option protection against downside tail risks is greater for firms with more carbon-intensive business models.

The literature related to the link between climate and financial risk is rather recent and some of it is controversial. In a recent review, Chakrabarty and Nag (2023) analyzed many articles on the estimation of carbon risk and hedging strategies based on Markowitz’s mean-variance formulation or capital asset pricing model (CAPM) approaches. Risk measures are classified into two broad groups, i.e., Fama–French-extended (FFE) models and linear risk models (LRMs). The conclusion is that there is a lack of consensus about the effects of carbon risk on stock prices. Campiglio et al. (2023) developed a critical analysis of the literature concerning the impact of climate risks on financial asset prices. They discussed studies showing that both physical and transition risks can trigger a revaluation of financial assets through multiple channels. Their analysis suggests that climate-related risks can have important implications on financial stability. Battiston et al. (2017) analyzed the impact of climate change on the value of assets held by banks and financial companies. By adopting a network approach, the authors examined how climate policy risk can propagate through the financial system. Thus, climate change risk may potentially pose systemic threats to global financial stability.

Recent studies have proposed new approaches to climate financial risk; for example, see the natural capital analysis and climate VaR by Dietz et al. (2016). Garcia-Jorcano and Sanchis-Marco (2024) analyzed the connection between global and regional mean sea level rise, used as proxies for climate risk, and financial market risk using VaR and a coherent risk measure based on quantile and expectile regression methods. Their findings indicate that measures of sea level rise have contrasting effects on different financial and economic sectors. In particular, the insurance sector exhibits the highest risk premium, while the oil and gas sector bears the highest risk cost.

Among the risk measurement procedures currently used by institutions and supervisors (policymakers, regulators, and investors) to mitigate future physical risk, some commonly used practices include risk scores, scenario analysis, stress testing, and sensitivity analysis. Finally, we should mention that the increase in severe weather events predicted by most climate scientists has a significant impact on the insurance industry as well. Insurers must adapt pricing strategies, use innovative technologies, and develop new models to measure, predict, and mitigate risks.

In our work, we put forward a methodology to calculate the impact of an increase in the global surface temperature on the probability of default (PD) of a company from the agriculture sector. As pointed out in Campiglio et al. (2023) and the references therein (Anttila 2016; Choi et al. 2020; Griffin et al. 2019), extreme temperatures have a negative impact on asset prices in all sectors. Moreover, Cuculiza et al. (2023) measured a firm sensitivity to temperature changes by regressing a stock’s excess return on the market factor and the temperature anomaly variable.

We performed a regression analysis of the firm’s stock value using the S&P 500 (SPX) and temperature anomalies (TAs) as predictors. The results indicate that an increase in TA has a negative impact on the stock value of the firm considered in this study. When modeling the TA time series, it is crucial to identify points in time where significant changes occur in the data behavior. These points are called breakpoints (BPs). We modeled the TA time series using segmented linear regression (SLR), with BPs estimated through wavelet analysis. We further investigated whether quantile regression provides a better model fit than wavelet analysis, but concluded that this was not the case for modeling TA data series. We calibrated the Merton model using real stock data from the company and estimated the PD based on the likelihood that the firm’s asset value falls below its liabilities. A sensitivity analysis is given, facilitated by the closed-form PD formula. We used the regression model to forecast future stock values influenced by predicted TA and estimated a new PD. It is worth noting that the research method presented in this work is limited to publicly listed companies, as it relies on the firm’s stock values quoted in the market. For non-listed companies, a similar methodology could be used based on the firm’s accountability. Finally, transition risk is studied and actions on climate-related risks are recommended.

This paper is organized as follows. In Section 2, we present the econometric model as well as the wavelet analysis on TA. Section 3 is devoted to the estimation of a firm’s PD based on the Merton model. In Section 4, we propose a model to study to what extent a listed firm is exposed to transition risk. Section 5 gives some recommendations to policymakers, risk managers, and investors in regard to the analysis carried out. Finally, Section 6 concludes the article.

2. Econometric Model and Temperature Anomaly Modeling

In order to provide a real-world example, we model the stock of a company with activities in the agricultural sector. For the sake of anonymity, we refer to the company by the acronym ‘ABF’, which stems from an agribusiness firm. The agricultural sector is an important component of the economy and it might be impacted by climate change. The explanatory variables of our regression model are SPX and TA. We use monthly data for temperature anomalies, measured in degrees Celsius, gathered from www.ncei.noaa.gov/ (accessed on 2 November 2023) (the National Centers for Environmental Information), as in Cuculiza et al. (2023). Monthly financial data for the quoted SPX and ABF stocks were obtained from Yahoo Finance (accessed on 2 November 2023).

In Figure 1, we plot the global land and ocean TA corresponding to the period from January 1850 to September 2023. The concept of temperature anomaly refers to a departure from a long-term average taken as a reference value. If the value of the anomaly is positive, then it indicates that the observed temperature is higher than the value taken as a reference. If the value of the anomaly is negative, then the observed temperature is below the reference value. The global time series was produced from the Smith and Reynolds land and ocean dataset (see Smith et al. 2008). The dataset consists of monthly average temperature anomalies measured on a 5° × 5° grid across ocean surfaces and land. The grid boxes are averaged to provide a mean global temperature anomaly. The anomalies are provided with respect to the base period 1901–2000.

It is worth remarking that we consider the global TA data series since ABF is a firm headquartered in the U.S. and spread all around the world. However, for more localized companies, it would be more accurate to consider the TA data series of that particular region.

The econometric model considers data within the period from January 2006 to September 2023,

$$y_i={\beta }_0+{\beta }_1{x}_{1i}+{\beta }_2{x}_{2i}+{ϵ}_i,i=1,\dots ,n,$$

(1)

where $y_i$ represents ABF equity data, $x_{1i},x_{2i}$ represent SPX and TA, respectively, and ${ϵ}_i$ denotes statistical noise. The estimated coefficients are summarized in

Table 1. Estimation summary of the econometric model (1).

Coefficients	Estimate	Std. Error	p-Value
Intercept (${\beta }_0$)	83.57971	4.79491	<2.0 $\times {10}^{-16}$
SPX (${\beta }_1$)	0.01009	0.00129	2.48 $\times {10}^{-13}$
TA (${\beta }_2$)	−41.37281	7.24330	3.79 $\times {10}^{-8}$

.

As we can observe, an increase in TA has a negative impact on the stock value of ABF. We model in the next section the time series of TA through SLR and forecast its future values to predict the impact of an increase of TA in ABF stock values.

2.1. Temperature Anomaly Modeling with SLR

In this section, we model the time series of TA through SLR, which is also employed in Mudelsee (2019), where the author estimates a unique breakpoint by moving block bootstrap resampling. Segmented regression makes sense when there are meaningful breakpoints for TA, allowing us to analyze changes in trends (and slopes when it comes to segmented linear regression) as well as handle the presence of structural breaks in non-stationary data. We will estimate the breakpoints using the Haar wavelets theory. Nonparametric regression with wavelets is explained in Section 3.1 of Abramovich et al. (2000). The sample size must be a power of 2, although this restriction can be overcome by considering some modifications (see Abramovich et al. (2000) and the references therein). In our case, the largest power of two samples contained has a length of 2048, corresponding to the period February 1853–September 2023. For the sake of completeness, we provide some details on the regression via Haar wavelets (for more details, refer to Abramovich et al. (2000)).

We consider a nonparametric regression,

$$y_i=g\left(t_i\right)+{ϵ}_i,i=1,\dots ,n,$$

where ${ϵ}_i$ are independent random variables with mean zero and variance ${\sigma }^2$. We aim at recovering the unknown function g from the data $y_i$ without assuming any particular a priori parametric structure for g.

A classical approach in regression analysis involves considering the function g expanded as a Fourier series and estimating the Fourier coefficients from the given data. The selection of the basis for the expansion is a crucial step. Ideally, the basis should be parsimonious in the sense that a wide set of potential response functions can be well approximated with only a few terms of the expansion. Wavelet series have remarkable approximation properties and allow a parsimonious expansion for a wide variety of functions. The simplest wavelet basis for $L^2\left(\mathbb{R}\right)$ is the Haar basis given by the following:

$$\varphi \left(t\right)=\left\{\begin{array}{cc}& 1,\mathrm{if}t\in [0,1],\hfill \\ & 0,\mathrm{otherwise},\hfill \end{array}\right.\psi \left(t\right)=\left\{\begin{array}{cc}& 1,\mathrm{if}t\in [0,1/2),\hfill \\ & -1,\mathrm{if}t\in [1/2,1],\hfill \\ & 0,\mathrm{otherwise},\hfill \end{array}\right.$$

where $\varphi$ is called the father wavelet or scaling function, and $\psi$ is called the mother wavelet. Other wavelets in the basis are generated by dilations and translations of $\psi$, that is, ${\psi }_{jk}\left(t\right)=2^{j/2}\psi (2^{j}t-k)$, called wavelet functions. This is the orthonormal wavelet basis for functions $g\in L^2\left(\mathbb{R}\right)$. We plot in Figure 2 the scaling function $\varphi$ and the mother wavelet $\psi$.

From now on, and without loss of generality, we can assume that the points $t_i={\displaystyle \frac{i}{n}}$ are equally spaced within the unit interval where the sample size n is a power of 2, that is, $n=2^J$ for some positive integer J.

In this work, we consider the so-called linear wavelet estimators of g (in contrast, nonlinear estimators can be used, as explained in Section 3.1.2 of Abramovich et al. (2000)). Suppose that g is expanded as a wavelet series on the interval $[0,1]$,

$$g\left(t\right)=c_0\varphi \left(t\right)+\sum _{j=0}^{\infty }\sum _{k=0}^{2^j-1}{w}_{jk}{\psi }_{jk}\left(t\right),$$

with $c_0=\langle g,\varphi \rangle$ and $w_{jk}=\langle g,{\psi }_{jk}\rangle$, where $\langle f,h\rangle :={\int }_0^{1}f\left(t\right)h\left(t\right)dt$ represents the $L^2$ inner product. Then, $c_0$ is called the scaling coefficient, while $w_{jk}$ denotes wavelet coefficients. Since we cannot estimate an infinite set of coefficients $w_{jk}$, we assume that g is well approximated by a finite set of basis functions,

$$g\left(t\right)\approx c_0\varphi \left(t\right)+\sum _{j=0}^M\sum _{k=0}^{2^j-1}{w}_{jk}{\psi }_{jk}\left(t\right),$$

for some $M<J$, called the level of approximation. The corresponding wavelet estimator ${\widehat{g}}_M\left(t\right)$ is of the following form:

$${\widehat{g}}_M\left(t\right)={\widehat{c}}_0\varphi \left(t\right)+\sum _{j=0}^M\sum _{k=0}^{2^j-1}{\widehat{w}}_{jk}{\psi }_{jk}\left(t\right),$$

where the sample estimates of the scaling coefficient and the wavelet coefficients are given by the following:

$${\widehat{c}}_0={\displaystyle \frac{1}{n}}\sum _{i=1}^n\varphi \left(t_i\right)y_i,{\widehat{w}}_{jk}={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\psi }_{jk}\left(t_i\right)y_i.$$

(2)

As pointed out by Abramovich et al. (2000), the performance of the estimator ${\widehat{g}}_M$ relies on the appropriate choice of level M. Intuitively, the optimal choice of M is related to the regularity of the response function g. A small value of M is associated with an over-smoothed estimator, while $M=J-1$ would simply reproduce the data.

We aim to estimate the breakpoints of the TA series by computing the wavelet coefficients of expression (2). A change in the size of these coefficients measured in absolute value indicates the presence of a jump in the time series. The intuition behind this fact is that, since wavelets ${\psi }_{jk}\left(t\right)$ are supported on the interval $\left[{\displaystyle \frac{k}{2^j}},{\displaystyle \frac{k+1}{2^j}}\right]$, each coefficient ${\widehat{w}}_{jk}$ of expression (2) is a weighted average of values $y_i$. These values are weighted by $2^j$ when $t_i\in \left[{\displaystyle \frac{k}{2^j}},{\displaystyle \frac{k+1/2}{2^j}}\right]$, and they are weighted by $-2^j$ when $t_i\in \left[{\displaystyle \frac{k+1/2}{2^j}},{\displaystyle \frac{k+1}{2^j}}\right]$.

We illustrate this fact in Figure 3 with the TA series of sample size $n=2048=2^{11}$, where $J=11$ and $M=0,\dots ,10$. Observe that at level $M=0$, there is only one coefficient, while at level $M=10$, the number of coefficients is 1024. In our case, we can identify four breakpoints by observing the largest coefficients within the first four approximation levels, and they correspond to the years 1938, 1981, 2002, and 2013. Some milestones stated in Weart (2008) confirm the trend changes during the following periods:

• In the 1930s, a global warming trend was reported, which started in the late nineteenth century.
• In 1938, CO₂ greenhouse global warming was underway, which revived interest in the question.
• Since the mid-1970s, strong global warming was reported, with 1981 being the warmest year on record.
• In 2001, warming was observed in ocean basins; the match with computer models gives a clear signal of the effect of greenhouse warming.

We carried out segmented linear regressions by considering different numbers of breakpoints and the corresponding plots are shown in Figure 4. For comparisons, we added a plot with linear regression (LR) and quadratic regression (QR).

We measure the difference between the true and predicted values through the root mean square error (RMSE), that is,

$$\sqrt{{\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}{\left(y_i-{\widehat{g}}_M\left(t_i\right)\right)}^2},$$

and we report these errors in

Table 2. RMSEs corresponding to segmented linear regressions with 1, 2, 3, and 4 BPs.

Method	HW		CE
	Year	RMSE	Year	RMSE
SLR 1 BP	1938	0.1603	1979	0.1506
SLR 2 BP	1938, 1981	0.1384	1937, 1994	0.1410
SLR 3 BP	1938, 1981, 2002	0.1379	1936, 1979, 2001	0.1334
SLR 4 BP	1938, 1981, 2002, 2013	0.1373	1936, 1979, 1997, 2014	0.1315

. Results given by the Haar wavelet (HW) method are compared with the cross-entropy (CE) method (see Priyadarshana and Sofronov 2015). The experiments corresponding to the CE method were performed with the R package breakpoint. The breakpoint obtained with both methods differs in the case of only two segments, while regressions with 1, 2, and 3 breakpoints give similar segmentations. In all cases, the RMSE is close when comparing HW and CE methods. As a reference, We provide the RMSE values of LR and QR (0.2218 and 0.1453, respectively). We observe that SLR values with 2, 3, and 4 breakpoints outperform QR in terms of RMSE. As pointed out in Priyadarshana and Sofronov (2015), the CE method for breakpoint detection is an iterative stochastic optimization method that starts with a parametrized distribution, from which a random sample is generated with respect to the number of breakpoints. The authors state that the overall processing time significantly increases with the increase in sample size. In contrast, HW is a nonparametric method and it simply relies on the computation of wavelet coefficients ${\widehat{w}}_{jk}$ given by expression (2); therefore, it is not based on a simulation. Further, since breakpoints are associated with coefficients at different levels, the estimated breakpoints remain when increasing the number of segments (when using the CE method, the breakpoints change depending on the number of segments considered).

We compute the coefficient of determination $r^2$ corresponding to the last segment for each SLR. We note that the last segment will be used in Section 3 to forecast TA values. The results obtained are shown in

Table 3. $r^2$ corresponding to the last segment.

Breakpoints	HW	CE
1	0.72	0.74
2	0.75	0.63
3	0.52	0.51
4	0.17	0.01

. We observe that the HW method gives, in general, better results than the CE method. The $r^2$ of the last segment, when considering four breakpoints, is much smaller for CE than for the HW method (this last segment only contains 10 years of monthly data).

Finally, in

Table 4. RMSEs for the last segment of the SLR with 4 BPs and the QRs for the same segment, along with the minimum and maximum values for a 1-year forecast using both models.

Model	RMSE	Min	Max
QR	0.0416	0.93	0.94
SLR 4 BP	0.0358	1.02	1.04

, we present the RMSE corresponding to the HW method for the last segment of SLR with four breakpoints and the RMSE of the quadratic regression on the same segment. We also present the minimum and maximum values of TA forecasted for the next 12 months (that is, from October 2023 until September 2024).

The most recent values of the TA series (from October 2023 to July 2024) are 1.37, 1.42, 1.38, 1.29, 1.41, 1.36, 1.29, 1.18, 1.23, and 1.21, which are more aligned with the forecast given by the SLR method with four breakpoints.

In this section, we show that segmented linear regression outperforms linear and quadratic regression when it comes to handling structural changes in the data. The next section is devoted to investigating whether quantile regression offers a better model fit than ordinary least squares regression when it comes to modeling TA time series.

2.2. Temperature Anomaly Modeling with Quantile Regression

Quantile regression consists of modeling the relation between two variables in the tails of the distribution. While classical linear regression minimizes the residuals in the least squares sense, quantile regression minimizes the loss function given by the following:

$$\underset{a\left(\tau \right),b\left(\tau \right)}{min}\sum _{i=1}^n{\rho }_{\tau }\left(y_i-a\left(\tau \right)-b\left(\tau \right)t_i\right),$$

where ${\rho }_{\tau }$ is defined as follows:

$${\rho }_{\tau }\left(x\right)=\left\{\begin{array}{cc}& -(1-\tau )x,x\le 0,\hfill \\ & \tau x,x>0,\hfill \end{array}\right.$$

and $\tau \in (0,1)$ indicates the level of quantile. We carried out the regressions corresponding to $0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,\mathrm{and}0.99$ quantiles and plot them in Figure 5. The slopes increase from quantile $0.3$ onward. This can be clearly observed in Figure 6, where we can also see that quantile regression estimates are not within the bounds of the linear regression estimates, suggesting a statistically significant difference between the two models.

For this reason, we perform a model comparison.

Table 5. AIC values corresponding to regressions for the $0.1$, $0.2$, $0.3$, $0.4$, $0.5$, $0.6$, $0.7$, $0.8$, $0.9$, and $0.99$ quantiles.

Quantile	AIC
0.1	255.58550
0.2	−14.25708
0.3	−129.24424
0.4	−139.04305
0.5	−76.57060
0.6	19.03958
0.7	161.60019
0.8	375.65011
0.9	744.52075
0.99	2014.62376

shows the AIC values corresponding to all quantiles represented in Figure 6. The AIC value for the ordinary least squares regression is −350.0208.

Based on AIC values, quantile regression at level 0.4 provides a better fit compared to the rest of the quantiles. However, the AIC value for linear regression is even lower, suggesting a better fit than quantile regression. Further, if we compute the mean absolute error for these two models, we have 0.1858 in the case of quantile regression while the error is 0.1815 for linear regression.

Finally, we conclude that the SLR put forward in this work provides a better model fit than linear and quantile regressions for TA data.

3. Calibrating the Default Probability of a Firm Based on Its Asset Values

The estimation of a firm’s default probability is a major concern for banks, investors, and the financial system as a whole, as it serves as a measure of credit quality. Supervisory authorities require banks to estimate probabilities of default to compute regulatory and economic capital for credit risk, among other risks. A complete list of publications by the Basel Committee on Banking Supervision (BCBS) is available at www.bis.org (accessed on 21 November 2023). The BCBS consists of senior representatives of bank supervisory authorities and central banks from Argentina, Australia, Belgium, Brazil, Canada, China, France, Germany, Hong Kong SAR, India, Indonesia, Italy, Japan, Korea, Luxembourg, Mexico, the Netherlands, Russia, Saudi Arabia, Singapore, South Africa, Spain, Sweden, Switzerland, Turkey, the United Kingdom, and the United States. The committee usually meets at the Bank for International Settlements (BIS) in Basel, Switzerland, where its permanent Secretariat is located. The PD is a key parameter for capital requirements and may be affected by climate risks, making its accurate estimation of the utmost importance (see BCBS 2021). Likewise, PD is crucial information for investors who may want to buy debt issued by firms.

In this section, we estimate the parameters of the Merton model with real stock data from ABF as well as its PD. Then, we use the econometric model put forward in Section 2, along with a forecast of the TA series, to obtain the climate PD.

The literature on estimating the probability of default (PD) of a firm distinguishes between two families of models: structural models and reduced-form models. While structural models explicitly relate to the firm’s equity, assets, and credit quality, reduced-form models rely on an exogenous intensity process to specify the default event. The dynamics of this process are calibrated to market data (e.g., credit spreads), but the economic conditions behind the process are unknown (see Spangler (2018) and the references therein for a discussion on the two types of models). In this work, we measure the impact on the PD of a firm via its equity values and, for this reason, it seems more natural to consider a structural model of default.

Merton’s model belongs to the class of structural default models and laid the foundation for factor models that have become the standard in credit risk measurement. In the Merton model, the event of default occurs when a firm’s asset value falls below a threshold representing its liabilities. This concept, known as the distance to default, was used by Nguyen et al. (2023) in the context of carbon emissions and transition risk. The Merton model provides a closed-form formula for the PD, making its computation straightforward A limitation of the basic Merton model is that it does not account for jumps in asset value; this issue can be addressed by using the Merton model with jumps. For the sake of completeness, we provide an overview of the seminal work on the Merton model Merton (1974) that can also be found in Section 3 of Lutkebohmert (2009) and Section 10.3 of McNeil et al. (2015).

The Merton model assumes that the asset value of a company is given by a stochastic process, which we denote by $\left(V_t\right)$. In the market, it is assumed that the firm has only equity (used interchangeably with the term stock) and debt. Further, it is assumed that the equity pays no dividends and the firm cannot issue new debt. The debt of the company is modeled by a zero-coupon bond with the value of B at a future maturity time T. The firm defaults when the value of its assets is less than the promised debt repayment at expiry T. We denote by $S_t$ and $B_t$ the values at time t of equity and debt, respectively. If we assume that there are no taxes or transaction costs, the value of the firm’s assets is given by $V_t=S_t+B_t,0\le t\le T$. At maturity time, there are only two possibilities:

(i)

$V_T>B$: The value of the firm is above the value of the debt. In this case, the owners of the zero-coupon bonds receive $B_T=B$; the shareholders receive the residual value $S_T=V_T-B$; and there is no default.

(ii)

$V_T\le B$: The value of the firm’s assets is below the debt value. Hence, the firm cannot meet its financial obligations and defaults. In this second case, the shareholders are left with nothing because the owners of the debt take control of the company, so we have $B_T=V_T,S_T=0$.

Combining the results, the payoff at maturity T from the point of view of the shareholders is given by the following:

$$S_T=max(V_T-B,0)={[V_T-B]}^{+},$$

(3)

and the owners of the zero-coupon bonds receive the following:

$$B_T=min(V_T,B)=B-{[B-V_T]}^{+}.$$

(4)

This shows that the value of the firm’s equity at expiry $t=T$ can be represented as the payoff of a European call option on the firm’s assets, with a strike price equal to the promised debt B at maturity. By using the put-call parity argument, the firm’s debt is equivalent to a risk-free bond that guarantees the payment of B, plus a short European put option on the firm’s assets with a strike price of B. The Merton model considers the asset value, $V_t$, as the underlying asset and assumes that under the real-world (or physical) probability measure P, the dynamics of the asset value $\left(V_t\right)$ are given by a geometric Brownian motion of the following form:

$$d{V}_t={\mu }_V{V}_{t}dt+{\sigma }_V{V}_{t}d{W}_t,0\le t\le T,$$

(5)

for constants ${\mu }_V\in \mathbb{R}$, ${\sigma }_V>0$, and a standard Brownian motion $\left(W_t\right)$. The solution at time $t=T$ of the stochastic differential Equation (5) with initial value $V_0$ can be obtained explicitly and is given by the following:

$$V_T=V_0{e}^{\left({\mu }_V-{\displaystyle \frac{1}{2}}{\sigma }_V^2\right)T+{\sigma }_V{W}_T}.$$

Since $W_T\sim N(0,T)$, we have the following:

$$ln{V}_T\sim N\left(ln{V}_0+\left({\mu }_V-{\displaystyle \frac{1}{2}}{\sigma }_V^2\right)T,{\sigma }_V^{2}T\right),$$

and then $ln{V}_T$ follows a normal distribution with a mean of $ln{V}_0+\left({\mu }_V-{\displaystyle \frac{1}{2}}{\sigma }_V^2\right)T$ and standard deviation of ${\sigma }_V\sqrt{T}$. The probability of default of the company by time $t=T$ is the probability that shareholders decide not to exercise the call option to buy the company’s assets at the price B at time T, that is, it is the probability that the call option expires out of the money. It can be computed as follows:

$$P(V_T\le B)=P(ln{V}_T\le lnB)=\Phi \left({\displaystyle \frac{ln\frac{B}{V_0}-({\mu }_V-\frac{1}{2}{\sigma }_V^2)T}{{\sigma }_V\sqrt{T}}}\right),$$

(6)

where $\Phi$ denotes the cumulative standard normal distribution. In the context of the Merton model, we can price the firm’s debt and the firm’s equity provided that we make the following assumptions:

(i)

The risk-free interest rate r is deterministic and positive.

(ii)

The firm’s asset value process $\left(V_t\right)$ is independent of the way the company is financed and it is also independent of the level of the debt B (this assumption makes sense thanks to the Modigliani–Miller theorem in Modigliani and Miller (1958)).

(iii)

The asset value $\left(V_t\right)$ can be traded on the market without friction, and the dynamics of the asset is given by the geometric Brownian motion (5).

Under these assumptions, the risk-neutral pricing theory then yields that the market value of equity at time $t<T$ can be computed as the discounted expectation of the payoff function (3), that is, we have the following:

$$S_t={\mathbb{E}}_{\mathbb{Q}}\left[e^{-r(T-t)}{[V_T-B]}^{+}\right],$$

given by the following:

$$S_t=V_t·\Phi \left(d_{t,1}\right)-B·e^{-r(T-t)}·\Phi \left(d_{t,2}\right),$$

where

$$d_{t,1}={\displaystyle \frac{ln(V_t/B)+(r+\frac{1}{2}{\sigma }_V^2)(T-t)}{{\sigma }_V\sqrt{T-t}}},d_{t,2}=d_{t,1}-{\sigma }_V\sqrt{T-t}.$$

Also, according to Equation (4), we can value the firm’s debt at time $t\le T$ as follows:

$$\begin{array}{cc}\hfill B_t& ={\mathbb{E}}_{\mathbb{Q}}\left[e^{-r(T-t)}(B-{[B-V_T]}^{+})\right]\hfill \\ \hfill & =B{e}^{-r(T-t)}-\left(B{e}^{-r(T-t)}\Phi (-d_{t,2})-V_t\Phi (-d_{t,1})\right).\hfill \end{array}$$

(7)

Note that expression (7) can be written in an equivalent form as follows:

$$B_t=V_t\left(1-\Phi \left(d_{t,1}\right)\right)+B{e}^{-r(T-t)}\Phi \left(d_{t,2}\right).$$

(8)

Next, we compare real and risk-neutral probabilities of the default in Merton’s model. A basic result from financial theory states that—under the risk-neutral measure Q—the process $\left(V_t\right)$ satisfies the stochastic differential equation, as follows:

$$d{V}_t=r{V}_{t}dt+{\sigma }_V{V}_{t}d{\tilde{W}}_t,$$

(9)

for a standard Q-Brownian motion $\left(\tilde{W}\right)$. Observe how the value of the drift ${\mu }_V$ in (5) was replaced by the risk-free interest rate r. Then, the risk-neutral default probability is given by Formula (6), evaluated with ${\mu }_V=r$, as follows:

$$q=Q(V_T\le B)=\Phi \left({\displaystyle \frac{ln\frac{B}{V_0}-(r-\frac{1}{2}{\sigma }_V^2)T}{{\sigma }_V\sqrt{T}}}\right).$$

(10)

If we compare this with the real default probability $p=P(V_T\le B)$ as given in (6), we end up with the following relation:

$$p=\Phi \left({\Phi }^{-1}\left(q\right)-{\displaystyle \frac{{\mu }_V-r}{{\sigma }_V}}\sqrt{T}\right).$$

(11)

We will use Formula (11) to go from risk-neutral default probabilities to physical probabilities.

3.1. Sensitivity Analysis

In this section, we perform a mathematical analysis to show how the probability of default is affected when either the initial asset value of the firm $V_0$ or the asset volatility ${\sigma }_V$ changes. Parameters affecting the asset value are of paramount importance since the asset value depends on the equity of the firm which, in turn, is impacted by changes in temperature anomalies. The expression (6) for the physical PD of the firm by time $t=T$ is given by the following:

$$p(V_0,{\sigma }_V)=\Phi \left({\displaystyle \frac{ln\frac{B}{V_0}-({\mu }_V-\frac{1}{2}{\sigma }_V^2)T}{{\sigma }_V\sqrt{T}}}\right),$$

where we want to emphasize the dependence of the PD on parameters $V_0$ and ${\sigma }_V$. Thanks to the above closed-form formula, we can easily calculate the PD derivative with respect to the asset value $V_0$, as follows:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial p}{\partial V_0}}& =-{\displaystyle \frac{1}{{\sigma }_V\sqrt{T}{V}_0}}{\Phi }^{\prime }\left({\displaystyle \frac{ln\frac{B}{V_0}-({\mu }_V-\frac{{\sigma }_V^2}{2})T}{{\sigma }_V\sqrt{T}}}\right)\hfill \\ \hfill & =-{\displaystyle \frac{1}{{\sigma }_V\sqrt{2\pi T}{V}_0}}exp\left(-{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{ln\frac{B}{V_0}-({\mu }_V-\frac{{\sigma }_V^2}{2})T}{{\sigma }_V\sqrt{T}}}\right)}^2\right),\hfill \end{array}$$

where ${\Phi }^{\prime }$ denotes the derivative of $\Phi$ with respect to $V_0$. As we can observe, ${\displaystyle \frac{\partial p}{\partial V_0}}$ is strictly negative, that is, p is a decreasing function of $V_0$. Therefore, an increase in the equity value implies an increase in the firm’s initial asset value $V_0$, which, in turn, implies a decrease in the firm’s PD. Conversely, an increase in temperature anomalies has a detrimental effect on equity (shown by the econometric model (1)) and asset values and a potential increase in PD.

If we consider the derivative with respect to the asset value volatility ${\sigma }_V$, we end up with the following:

$${\displaystyle \frac{\partial p}{\partial {\sigma }_V}}={\Phi }^{\prime }\left({\displaystyle \frac{ln\frac{B}{V_0}-({\mu }_V-\frac{{\sigma }_V^2}{2})T}{{\sigma }_V\sqrt{T}}}\right)\left[{\displaystyle \frac{ln{V}_0-lnB}{{\sigma }_V^2\sqrt{T}}}+\left({\mu }_V+{\displaystyle \frac{{\sigma }_V^2}{2}}\right){\displaystyle \frac{\sqrt{T}}{{\sigma }_V^2}}\right].$$

This function is strictly positive when $V_0>Bexp(-({\mu }_V+{\displaystyle \frac{{\sigma }_V^2}{2}})T)$, that is, p is an increasing function of ${\sigma }_V$, and it is decreasing otherwise. A well-known result derived from the Itô formula states that the equity volatility ${\sigma }_E$ is proportional to the asset value volatility,

$${\sigma }_E={\displaystyle \frac{\partial S_t}{\partial V_t}}{\displaystyle \frac{V_t}{S_t}}{\sigma }_V.$$

Hence, an increase in the equity volatility ${\sigma }_E$ implies an increase in the firm’s asset value volatility ${\sigma }_V$, which, in turn, implies an increase in the firm’s PD when $V_0>Bexp(-({\mu }_V+{\displaystyle \frac{{\sigma }_V^2}{2}})T)$. Again, structural breaks in the TA data series increase the equity volatility and, therefore, the PD.

3.2. Estimation of PD with Real Market Data

Let $\mathbf{s}=\left(s_1,s_2,\dots ,s_m\right)$ be the stock time series of ABF and ${\mathbf{b}}^{\left(0\right)}=\left(b_1^{\left(0\right)},b_2^{\left(0\right)},\dots ,b_m^{\left(0\right)}\right)$ the corresponding liabilities one year ahead, where bold letters denote vectors. The liabilities were obtained as the product of the debt-to-equity ratio obtained from www.macrotrends.net (accessed on 21 November 2023) and the ABF stock value on the same date. Since the debt-to-equity ratios were only publicly available for the period from March 2009 to September 2023, we consider the stock $\mathbf{s}$ within the dates March 2008 to September 2022 ($m=175$ values of monthly data). Then, we can estimate the probabilities of default from March 2009 to September 2023 on a monthly basis. It is worth remarking that, in a practical situation, a more reliable and accurate estimation of one-year-ahead liabilities can be done by using the internal information of the firm.

We need to first calibrate the asset value process of expression (5) using the real data, $\mathbf{s}$ and ${\mathbf{b}}^{\left(0\right)}$, mentioned before. The calibration of the stochastic process in (5) is carried out by estimating parameters ${\mu }_V$ and ${\sigma }_V$. We initialize the asset value of ABF by considering the following:

$${\mathbf{v}}^{\left(0\right)}=\mathbf{s}+{\mathbf{b}}^{\left(0\right)}{e}^{-rT},$$

(12)

where $r>0$ is the risk-free interest rate that we set equal to 0.04 and the $T=1$ year. Then, ${\mathbf{v}}^{\left(0\right)}$ has components ${\mathbf{v}}^{\left(0\right)}=\left(v_1^{\left(0\right)},v_2^{\left(0\right)},\dots ,v_m^{\left(0\right)}\right)$.

It is well-known that when $\left(V_t\right)$ follows a geometric Brownian motion, the log asset returns $r_{\Delta t}:=ln\left({\displaystyle \frac{V_{t+\Delta t}}{V_t}}\right)$ are normally distributed with expectation $({\mu }_V-{\displaystyle \frac{{\sigma }_V^2}{2}})\Delta t$ and variance ${\sigma }_V^2\Delta t$ for a certain time step $\Delta t$. If we define $\widehat{\mathbb{E}}\left[r_{\Delta t}\right]$ and $\widehat{\mathbb{V}}\left[r_{\Delta t}\right]$ as the sample mean and sample variance, respectively, of empirical log returns, then we can provide estimations ${\widehat{\mu }}_V$ and ${\widehat{\sigma }}_V$ for ${\mu }_V$ and ${\sigma }_V$, respectively, as follows:

$${\widehat{\mu }}_V={\displaystyle \frac{\widehat{\mathbb{E}}\left[r_{\Delta t}\right]}{\Delta t}}+{\displaystyle \frac{\widehat{\mathbb{V}}\left[r_{\Delta t}\right]}{2\Delta t}},{\widehat{\sigma }}_V^2={\displaystyle \frac{\widehat{\mathbb{V}}\left[r_{\Delta t}\right]}{\Delta t}},$$

(13)

where $\Delta t=1/12$ is considered for monthly data.

We can now give an initial estimation ${\widehat{\sigma }}_V^{\left(0\right)}$ of ${\sigma }_V$ through expression (13) with real data given in expression (12). The initial estimation of volatility will be refined by the following procedure:

• Compute a new data vector ${\mathbf{b}}^{(k+1)}$, starting from $k=0$, using expression (8) with the $T-t=1$ year, as follows:

  $$b_i^{(k+1)}=v_i^{\left(k\right)}\left(1-\Phi \left(D_{i,1}^{\left(k\right)}\right)\right)+b_i^{\left(0\right)}{e}^{-r(T-t)}\Phi \left(D_{i,2}^{\left(k\right)}\right),$$

  where

  $$\begin{array}{cc}\hfill D_{i,1}^{\left(k\right)}& ={\displaystyle \frac{ln(v_i^{\left(k\right)}/b_i^{\left(0\right)})+\left(r+\frac{1}{2}{\left({\widehat{\sigma }}_V^{\left(k\right)}\right)}^2\right)(T-t)}{{\widehat{\sigma }}_V^{\left(k\right)}\sqrt{T-t}}},\hfill \\ \hfill D_{i,2}^{\left(k\right)}& =D_{i,1}^{\left(k\right)}-{\widehat{\sigma }}_V^{\left(k\right)}\sqrt{T-t},i=1,\dots ,m.\hfill \end{array}$$
• Compute a new data vector ${\mathbf{v}}^{(k+1)}=\mathbf{s}+{\mathbf{b}}^{(k+1)}.$
• Compute a new estimate of volatility ${\widehat{\sigma }}_V^{(k+1)}$ using expression (13) and the new vector ${\mathbf{v}}^{(k+1)}$.
• Given a tolerance $ϵ$, we iterate the process until

  $$\parallel {\mathbf{v}}^{(k+1)}-{\mathbf{v}}^{\left(k\right)}{\parallel }_2<ϵ,$$

  (14)

  where ${\parallel ·\parallel }_2$ denotes the $L^2$ norm.

Let ${\mathbf{v}}^{\left(l\right)}$ and ${\widehat{\sigma }}_V^{\left(l\right)}$ be the outcomes of the procedure given above, where l is the number of iterations until condition (14) is satisfied. Then, the risk-neutral probabilities can be estimated with the formula of expression (10) by computing the following:

$${\widehat{q}}_i=\Phi \left({\displaystyle \frac{ln\frac{b_i^{\left(0\right)}}{v_i^{\left(l\right)}}-\left(r-\frac{1}{2}{\left({\widehat{\sigma }}_V^{\left(l\right)}\right)}^2\right)T}{{\widehat{\sigma }}_V^{\left(l\right)}\sqrt{T}}}\right),i=1,\dots ,m.$$

(15)

Finally, the physical default probabilities are calculated with formula of expression (11), the risk-neutral probabilities of expression (15), and the estimation ${\widehat{\mu }}_V$ of ${\mu }_V$ through expression (13), that is, we have the following:

$${\widehat{p}}_i=\Phi \left({\Phi }^{-1}\left({\widehat{q}}_i\right)-{\displaystyle \frac{{\widehat{\mu }}_V-r}{{\widehat{\sigma }}_V^{\left(l\right)}}}\sqrt{T}\right),i=1,\dots ,m.$$

(16)

Note that physical default probabilities can be obtained straightforwardly from expression (6), and there is no need to calculate risk-neutral probabilities of expression (15). We compute both to illustrate the so-called Sharpe ratio in expression (16), given by ${\displaystyle \frac{{\widehat{\mu }}_V-r}{{\widehat{\sigma }}_V^{\left(l\right)}}}$. The Sharpe ratio represents the risk premium of an investment versus a safe asset, such as a bond. The overall process is summarized in Algorithm 1.

Algorithm 1: PD estimation

We perform several experiments by varying $ϵ$, where the initial volatility is ${\widehat{\sigma }}_V^{\left(0\right)}=0.291775$. We observed rapid convergence, as ${\sigma }_V^{\left(l\right)}=0.287022$ for all values of $ϵ$ considered in

Table 6. Number of iterations l for estimating the volatility following Algorithm 1 when different tolerance values $ϵ$ are considered.

$\mathit{ϵ}$	l
1 $\times {10}^{-2}$	7
1 $\times {10}^{-4}$	11
1 $\times {10}^{-6}$	15
1 $\times {10}^{-8}$	19
1 $\times {10}^{-10}$	23

.

We plot the physical probabilities of the default in Figure 7 and observe a peak around May 2020. The drift is ${\widehat{\mu }}_V=0.060925$ and the average PD corresponding to the last 12 months is 0.04.

Now, we estimate the climate PD until September 2025. To do this, we consider two years more of ABF stock values, that is, from October 2022 to September 2024. While the first year (October 2022–September 2023) comprises real stock, the stock of the second year (October 2023–September 2024) is estimated through the econometric model of expression (1). The forecast of ABF stock values for the second year is conducted using the same SPX values from the first year and the forecasted TA values obtained with the SLR 4 BP model (we chose not to forecast SPX to isolate the impact of the TA forecast). Similarly, we use the most recent one-year liabilities (October 2022–September 2023) for this two-year period.

If we apply Algorithm 1 to the entire data series of ABF stock (March 2008-September 2024), we obtain initial volatility ${\widehat{\sigma }}_V^{\left(0\right)}=0.286864$ and ${\widehat{\sigma }}_V^{\left(l\right)}=0.282564$ (where $l=22$ for $ϵ$ = 1 $\times {10}^{-10}$). The drift is ${\widehat{\mu }}_V=0.057624$ and the forecasted average climate PD corresponding to the last 12 months is 0.06. The increase in PD might be associated with a downgrade in the company’s rating. Referring to the S&P rating-grade nomenclature and the PD ranges reported in Table 2 of Gordy and Lütkebohmert (2013), a 4% PD would correspond to a BB grade, while a 6% PD would be linked to a B grade.

4. Transition Risk

Transition risk has become increasingly important in the context of climate change. It is associated with the transition to a low-carbon economy. Such a transition requires changes in policy, legislation, technology, and market dynamics. We believe it is necessary to provide a quantitative approach to defining transition risk in the context of standard financial theory. In line with this approach, we propose a model in this section for evaluating whether ABF qualifies as a green company. This assessment is based on statistical metrics and asset pricing rationale.

We consider the following model for the log returns of the ABF stock:

$$ln\left({\displaystyle \frac{y_i}{y_{i-1}}}\right)=a+bln\left({\displaystyle \frac{x_{1,i}}{x_{1,i-1}}}\right)+d\left(ln\left({\displaystyle \frac{G_{e,i}}{G_{e,i-1}}}\right)-ln\left({\displaystyle \frac{G_{a,i}}{G_{a,i-1}}}\right)\right)+{ϵ}_i,$$

(17)

where y represents the ABF stock, $x_1$ represents the SPX, $G_e$ represents an index of low-polluting companies (green stocks), and $G_a$ represents an index of highly polluting companies (gray companies). Our measure of transition risk is given by coefficient d. This coefficient can be interpreted as the sensitivity of the stock return series to the premium of a green stock index compared to an index of highly polluting companies. When $d>0$, this implies that the returns of ABF stock increase when $ln\left({\displaystyle \frac{G_{e,i}}{G_{e,i-1}}}\right)-ln\left({\displaystyle \frac{G_{a,i}}{G_{a,i-1}}}\right)$ increases, suggesting that the company is already green and the transition risk is low. If $d<0$, then the returns of ABF stock decrease when $ln\left({\displaystyle \frac{G_{e,i}}{G_{e,i-1}}}\right)-ln\left({\displaystyle \frac{G_{a,i}}{G_{a,i-1}}}\right)$ increases, implying that the company is not green and that the transition risk is high.

For the $G_e$ index, we use the Vanguard ESG U.S. Stock ETF (ESGV), which specifically excludes stocks of certain companies related to adult entertainment, alcohol, tobacco, cannabis, gambling, chemical and biological weapons, cluster munitions, anti-personnel landmines, nuclear weapons, conventional military weapons, civilian firearms, nuclear power, and coal, oil, or gas. It also excludes stocks of companies that fail to meet certain labor, human rights, environmental, and anti-corruption standards as defined by the United Nations Global Compact Principles. Additionally, it excludes companies that do not meet specific diversity criteria. To the best of our knowledge, there are no suitable index proxies for $G_a$, but as a proxy for our empirical analysis, we use the Exxon Mobile stock (XOM). Exxon Mobile is one of the highest polluting companies worldwide and it is a public company listed in the U.S. The Texas-based firm’s three largest refineries (two in Texas and one in Louisiana) are the nation’s top three emitters of small particulate matter, according to the analysis of the latest tests submitted to regulators by the ten largest refineries in the U.S. The three Exxon refineries together averaged emissions of 80 pounds per hour, eight times the average rate of the seven other refineries on the top-ten list, some of which are larger than Exxon’s plants, as per the analysis available at www.reuters.com (accessed on 21 November 2023).

For this analysis, we consider daily data corresponding to the period 2 January 2019–15 November 2023. The estimated coefficients of model (17) are $a=0.0002144,b=0.8584642$, and $d=-0.3338892$; therefore, we conclude that ABF is not a green company and there exist potential transition risks. The transition risk parameter d is suitable for making comparisons between different companies or sectors of the economy and aims to measure the exposure to transition risk. We provide quantitative tools to support and complement qualitative approaches for assessing transition risk.

5. Actionable Insights for Policymakers, Risk Managers, and Investors

The research carried out in Section 3 and Section 4 reveals that ABF is impacted by physical risks and is exposed to transition risks. As pointed out in the introduction, companies with high carbon footprints are perceived by the market as more likely to default.

A higher PD for a firm makes it more difficult and expensive to issue debt, as investors demand a higher return for taking on additional risk. The more polluting the firm, the higher its cost of capital due to the increased cost of debt. As the cost of capital rises, it becomes more difficult for the firm to finance its operations. This would force the polluting companies to be more efficient in their operations and to pollute less. The higher cost of capital would also impact company valuations. Most fundamental investors use the weighted average cost of capital (WACC) to discount free cash flows and determine company valuations. If the cost of capital increases (due to higher debt costs), the price of equity will decline. This has a desirable effect on investors, as it rewards those who invest in green companies and penalizes those who invest in gray stocks. Another desirable effect of higher PDs related to climate factors is that the increased cost of debt would incentivize companies to issue shares instead of taking on more debt. This, in turn, would reduce the leverage effect of debt in polluting companies, leading them to pay more taxes than their non-polluting counterparts.

In regard to the findings featured in Section 4, policymakers and governments could establish tax and subsidy schemes for each public company with respect to their d value. For a positive d, companies could receive a green bond, while for a negative d, companies should pay more taxes. These green bonds could also be traded in open markets. Risk managers should adjust credit rating evaluations to account for potential increases in PD due to physical risk and rebalance their portfolios more frequently to respond to market changes.

6. Conclusions

We investigated the impact of TA on the PD of a firm in the agriculture sector, as well as its exposure to transition risk. We use an econometric model to estimate the firm’s stock value, with SPX and TA as predictors. The TA time series is modeled using wavelets and SLR. The firm’s future stock value is then predicted using the econometric model. Then, the Merton model is applied to calculate the PD of the firm, using both actual and forecasted stock values. Our experiments show that the initial estimated PD of 4% PD increases to 6% PD due to physical climate risk.

Future research will be devoted to extending the present work in two directions. From the methodological standpoint, we might consider using the Merton model with jumps for the asset value of the firm. A new trend change in the TA series could lead to sudden jumps in the firm’s asset value. While jumps allow for abrupt changes in asset value, incorporating them results in non-closed form formulas for PD, posing a mathematical challenge. From a risk management point of view, it would be worth collecting data from other companies in the agribusiness sector to verify whether the increase in PD is a systematic feature with an impact on the entire sector. Likewise, it would be interesting to consider other factors in the econometric model like, for instance, research and development expenses, which might help to reduce the impact of physical climate risk on the stock value of the firm.

Disclaimer: The probabilities of default for the agribusiness firm are calculated based solely on publicly available data and are intended for academic purposes only. We use a default model that can be fitted with that information, but a more accurate assessment would require additional relevant data as well as the expert judgment of risk managers. Therefore, the results presented in this work should be viewed as a methodological starting point rather than as financial advice.