The evaluation of the economic impact of global warming requires a balance of future costs against present costs [1]. The most basic assumption considers an exponential with a constant discount rate. The choice of discount rate is part on current economic debates on the urgency of the response to global warming (see [2][3][4][5]). Most normative approaches attempts to derive the discount from axiomatic principles of justice, or from utility theory and assumptions about growth [6,7]. Discounting includes impatience, economic growth, and declining marginal utility and it can embedded in the Ramsey formula [8][9][10]). This approach makes however difficult to make quantitative estimates because factors involved are difficult to be measured empirically. Here we take a positive approach assuming that discounting is equivalent to bond pricing. However, time to maturity of bonds is not generally higher than 30 years, but climate change time horizons need the discount 100 years or even more into the future. Thus, inferring long maturity bonds prices from empirical data on short maturity bonds is fundamental and long historical time series for having a reliable statistical inference are needed. This data is most often not available and then it will be required to consider a reasonable model for real interest rates at different maturities. There are two effects influencing long term rates that must be considered. The first of these is risk aversion. Longer term bonds bear greater risk than shorter ones, and higher risk should imply higher interest rates. The second effect is more subtle, and is due to the fact that interest rates are uncertain and highly persistent. In the environmental context, this effect was originally pointed out by Weitzman [11] and Gollier et al. [12] (see also [13]). However, the effect was also pointed out in a general context in [14] and it has been known much longer in the context of bond pricing [15]. The long run rate is thus dominated by the lowest value, since asymptotically all the other discount factors will be negligible in comparison. To see this concretely, consider two possible rates, r 1 and r 2 , with r 1 > r 2 , and assume they have equal probability. Then the average discount factor at time t in the future is D(t) = (e -r1t + e -r2t )/2. Note that since the sum of two exponentials is not an exponential, the discounting function is no longer an exponential. But for t sufficiently large D(t) ≈ e -r2t /2, i.e. the discount function becomes approximately exponential with the lower interest rate. This illustrates that when interest rates are uncertain and persistent, lower interest rates tend to dominate. We are here building on earlier empirical work. Stochastic interest rates provided by Litterman et al. [18], Newell and Pizer [19], and Groom et al. [20] are more realistic processes than those provided by Weitzman and Gollier. Groom et al. [20] also noted that the drop in rates depends on uncertainty and persistence. All these authors calibrated the models suggested against long historical time series of 10 year bonds. We hare add to these previous works in several ways. We study century-long records of the prices of 3 month bonds, 10 year bonds, and inflation, for both the US and the UK. Real interest rates are found to become negative more than one third of the time, often by substantial amounts. To our knowledge only few authors have addressed this issue (e.g. Freeman et al [21]) while the usual approach has been ignoring this fact, and instead forced historical interest rates to be positive in order to be consistent with their models. We instead take the view that negative real interest rates are a fact that cannot be ignored. This leads us to choose the Ornstein-Uhlenbeck (OU) model, which is compatible with negative real interest rates. The model provides real interest rates that can be negative but they are constantly pulled back to certain rate. The formula for the yield curve for the OU model has been derived in the finance literature using several approaches (see, for instance, [22][23][24]; and the review in [25]). We derive this by a somewhat simpler way using Fourier transform methods. The previous work in environmental discounting cited above [20] was based on numerical simulations. In contrast, we take advantage of the existence of closed form solutions, which allows us to better estimate the long term discount rate. When the OU model is used in the finance literature, the nominal interest rate and the mean reversion parameter is taken to be so strong that nominal interest rate goes negative with probability close to zero. In a different approach, Davidson et al. [26] examined a square root Ornstein Uhlenbeck model which can never go negative, is pulled down toward 0 as it is an absorbing state. Davidson et al. [26] has solved the asymptotic (long) rate by using the Feynman-Kac functional, which is quite different from our approach. The prevalence of negative real rates makes many of the standard nominal interest rate models inappropriate vehicles for studying real rates, and the OU becomes an December 29, 2023 obvious choice. In addition to our treatment of negative rates, we also give an original contribution regarding properly taking risk aversion into account. The works described above [19,20] were calibrated against a single maturity bond (10 years). Risk aversion is a well-accepted notion in bond pricing, and one would expect models that neglect it to underestimate long term interest rates. We fix this by taking risk aversion into account and fit the resulting model to both 3 month and 10 year bond yields, which provide us with two points on the yield curve. In summary, we herein present a complete empirical study on real interest rates of US and UK for which we have century-long records and that are considered ones of the most economically stable countries as far inflation and economic growth are concerned. Our focus is based on the interdisciplinary approach of complex systems developments.
We have also presented elsewhere a similar study for other less stable countries but with minor statistical reliability due to scarcity of data and without considering risk aversion (i.e., the market price of risk) [27]. Several different theoretical aspects of discounting have been recently discussed by some of us mainly regarding the effects on discounting of extreme market situations [28] or the introduction of resettings in real interest rates [29] (see also the recent review [30]). We now derive the form of the yield curve for the OU model. In the continuous limit the discount function is

where t 0 is the initial time and the expectation E[•] is an average over all possible real rate trajectories r(t ′ ) (t 0 ≤ t ′ ≤ t) up to time t. The mathematical expression is formally identical to the bond prices. The price B(t 0 |t 0 + t) of a zero-coupon bond issued at time t 0 with unit payoff and time maturity t 0 + t (t ≥ 0) reads

which depends on the nominal rate n(t). Following the same strategy as for bonds, we compute D(t) via a stochastic process model for the time evolution of r(t). Simplest and most common assumption is to consider r(t) as a diffusion process. The stochastic differential equation is

where f (r) gives reason of the drift while g(r) > 0 quantifies the noise intensity dW (t) = ξ(t)dt (a Wiener process) characterized with

and where δ(•) is Dirac's delta function. The discount function is therefore given by

December 29, 2023

where the cumulative process x(t) is defined as

Therefore, the discount function can written in the following form:

The function p(x, r, t|x 0 , r 0 , t 0 ) corresponds to the probability density function (PDF) that depends on (x(t), r(t)). The measure p is often referred to as the data generating measure. We note that when there are uncertainties in the value of certain parameters which may appear in the estimation of the PDF from empirical data then the average in Eq (6) has to be extended in order to include these additional uncertainties (see, e.g. [31]). Equations ( 3) and ( 5) together refer to a bidimensional process defined in terms of a pair of stochastic differential equations

These stochastic different equations imply that the joint density follows the Fokker-Planck equation (FPE)

The initial condition of this equation are fixed and are p(x, r, t 0 |x 0 , r 0 , t 0 ) = δ(x)δ(r -r 0 ). (10) The process r(t) is time homogeneous because f (r) and g(r) do not depend on time explicitly. The process is therefore time invariant and, without loss of generality, t 0 = 0. We now make explicit choices for the functions f and g. For the reasons given in the introduction, we focus on the OU process, which is a standard model from finance that allows negative rates. In finance the OU process was originally proposed by Vasicek [15] and is sometimes referred to as the Vasicek model. It is a diffusion process with linear drift and constant noise intensity

The process is thus governed by Eq (3). Letting r 0 = r(0) be the initial return, in the large time limit the probability density function p(r, t|r 0 ) has mean m and variance σ 2 = k 2 /2α. The Local Expectation Hypothesis (LEH) [32,33] assumes that prices can reasonably be modeled based on the data generating measure p. This results is attributed to the fact that inverstors are assumed to be risk neutral. However, if investors are sensitive to risk, bonds are priced with an artificial probability density function, p * , called either the risk-neutral probability measure or the risk-correcting measure (see Appendix B). The two measures p and p * are related by the market price of risk (see below), which is the extra return per unit of risk that investors demand to bear risk. Following a standard procedure for bond pricing [15,22,24,30] we take risk into account by adjusting the drift term according to

where in our case g(r) = k > 0 and q = q(r) ≥ 0 is the market price of risk. which, in general, may also depend on the current time, q = q(r, t), but we here assume stationarity. The introduction of the market price of risk raises the effective interest rate, which in the context of bond pricing means that investors are compensated for taking increased risk. The most common assumption consists in taking q(r, t) independent of r and t, that is q(r, t) = q = constant. In such a case the adjusted drift can be rewritten as

where

The result is that the effective mean interest rate m * is increased relative to the historically observed interest rate m by a constant amount that depends on the volatility, the reversion rate, and the market price of risk. For the case of constant price of risk, the joint PDF is found by solving the Fokker-Planck equation (9) , which now becomes (cf. Eqs (13))

with the initial condition provided in Eq. (10). The equation problem can be addressed with the characteristic function (Fourier transform)

e -iω2r p(x, r, t|r 0 )dr. ( 16) Equation ( 15) can therefore take a simpler form:

with (cf. Eq (10))

The solution is given by

where the expressions for the functions A(ω 1 , t), B(ω 1 , t), and C(ω 1 , t) are obtained in Appendix A. Comparing Eqs ( 6) and ( 16) we obtain

Finally, it can be shown in Appendix A that D(t) = exp{-C(-i, t)} which is a tractable expression for the discount rate (cf. Eq (36)). We now make a change of notation.
We have so far assumed that the discount functions is computed for a time t in the future, starting at a fixed initial time t 0 = 0 with initial interest rate r 0 for a time t in the future (let us remind that we have set t 0 = 0; had we kept t 0 ̸ = 0 then time t would have been replaced by t -t 0 ). However in what follows we will want to evaluate the discount starting from different initial times. Recalling that the whole process is invariant under time translations, we can thus perform the change of notation such that t 0 → t, r 0 → r(t) and t -t 0 → τ . Note that now t is the present time, r(t) = r is the present interest rate and the maturity time is t + τ , that is to say, τ is the "time to maturity". The discount function is then (cf. Eq (36

When the time to maturity τ is small we can approximate D(τ ) by expanding the exponentials to first order using a Taylor series approximation and this yields ln D(τ ) ≈ rτ (τ → 0), as expected. In the limit τ → ∞ we see that

and the discount function becomes independent of the initial rate. It decays exponentially and can be written in the form

where

is the long-run rate. In general the long-run rate may be defined by

as long as the limit exists. Note that Eq (22) can be readily obtained after substituting Eq (20) into Eq (23). From Eq (22) we see that the long-run discount rate depends on the historical rate m, but this is shifted by two terms. The first shift raises the long-run rate due to the market price of risk. The second shift lowers it by an amount given by the ratio of uncertainty (as measured by k) and persistence (as measured by α). We can trivially rewrite the equation above as

This makes it clear that whether or not the overall shift in the long-run discount rate is positive or negative depends on the size of the market price of risk in relation to the ratio of the volatility parameter and the reversion rate. It is not surprising that the market price of risk raises the long term rate, but it is not so obvious that uncertainty and persistence can lower it. Indeed for the OU process it can make it arbitrarily small. For any given mean interest rate m, by varying k and α, the long-run discount rate r ∞ can take on any value less than m, including negative values, while at the same time the standard deviation σ can also be made to take on any arbitrary positive value. In particular, by choosing the appropriate (k, α), we can make r ∞ arbitrarily far below m and σ arbitrarily small. The probability that r(t) < r ∞ can be arbitrarily small, even when r ∞ ≪ m (see Appendix A). Deducing the correct parameters (m, σ) of the stationary distribution of short run interest rates does not determine r ∞ by itself; on the contrary, any r ∞ < m is consistent with them. To infer r ∞ from the data one must also tease out the mean reversion parameter α. Holding the long run distribution (m, σ) constant, by raising the persistence parameter 1/α it is possible to lower r ∞ to any desired level. Of course, this can always be offset by the market price of risk. It is even possible for the long-run rate to be negative. This is due to the amplification of negative real interest rates r(t). Computation of the discount function involves an average over exponentials, rather than the exponential of an average. As a result, periods where interest rates are negative are amplified, and can easily dominate periods where interest rates are large and positive, even if the negative rates are rarer and weaker. It does not take many such periods to substantially reduce the long run interest rate. To summarize, in the OU model the long-run discounting rate can be much lower than the mean, and indeed can correspond to low interest rates that are rarely observed. 3 Empirical results

We estimate real interest rates as nominal rates corrected by inflation which we perform directly through the application of Fisher's equation by subtracting realized inflation from nominal interest rates (see Eq (26) below). This is not, however, the only way of obtaining real rates. Thus, for instance, Freeman et al. [21], among others, pursue an alternative procedure using cointegration methods to tease out real rates. Another way of obtaining real rates would be modeling nominal rates and inflation separately; that is, nominal rates by some positive random process (for example the Feller process as in the CIR model [17,34]) and the OU process for inflation since the latter can assume positive and negative values. Such a procedure -which unfortunately enlarge the number of parameters to be estimated from data-would also need to take into account possible correlations between bond prices and inflation. It is not clear a priori which procedure is better, and the direct procedure we follow has the virtue of being much simpler with lesser number of free parameters to estimate. We transform the annual rates into logarithmic rates and denote the resulting time series by y(t|τ ) (with the maturity time τ equal to either 3 month or 10 years). Nominal rates n(t) are then estimated by n(t) ∼ y(t|τ ) (see Appendix B for details). The inflation rate i(t) is estimated through the Consumer Price Index (CPI) as

where I(t) is the aggregated inflation up to time t, and τ = 10 years (cf. Appendix B). Finally, the real interest rate r(t) is defined by Fisher's equation:

We now estimate the OU model for real interest rates from historical data.
To this end we have collected long historical time series for both short and long run nominal interest rates, as well as inflation, for the United Kingdom and the United States. The properties of the data are summarized in Table 1. For each country we have both three month and ten year interest rates, as well as an inflation index. The recording frequency for each country is either annual or monthly. For ten year government bonds (which pay out over a period of ten years) we smooth inflation rates with a ten-year forward moving average, and subtract the annualized inflation index from the annualized nominal rate to compute the real interest rate. For the three-month bond rates, in contrast, we use the inflation adjustment for the corresponding year (since we do not have inflation adjustments at quarterly frequency). Figure 1 shows nominal rates, inflation and real interest rates for 3 month bonds for the UK and the US, and Figure 2 compares 3 month and 10 year real interest rates. This procedure assumes that people correctly forecast inflation, thus, in absence of any knowledge of behavioral bias, we are assuming perfect rationality. One of the most striking features of these time series is that real interest rates are often negative, in some cases by substantial amounts and for long periods of time (see Figures 1 and2). This is evident looking at the histograms shown in Figure 3. Real 3 month interest rates for the UK are negative 32 % of the time, and there are four distinct periods where they drop nearly reach -10 % and even below. 10 year real interest rates for the UK are negative 38 % of the time. For the US real 3 month interest rates are also negative 32 % of the time, dropping below -10 % during World War II, and 10 year real interest rates are negative 30 % of the time. Given that real interest rates are negative about a third of the time, this makes it clear that models such as the log-normal process or the CIR model [17,34] that assume rates to be non-negative are far from being appropriate. We therefore confine our empirical work to The curve compares to a normal distribution. This illustrates that negative real interest rates are common and that the distribution of interest rates is heavy-tailed relative to a normal distribution. the Ornstein-Uhlenbeck model. From Figure 3 it is also evident that the distribution of interest rates is heavy tailed. This is particularly true for the short term interest rates. Heavy tails are apparent because of the excess in the center of the distribution but also because the observations in the tail that exceed the normal distribution. However these deviations are not extreme and the OU process (which has a normal distribution) is at least a reasonable first approximation. Another striking feature is that the yield curve is intermittently inverted, that is, the 10 year real interest rate is at intervals lower than the 3 month rate (see Figure 1). For the UK the yield curve is inverted slightly more than 50% of the time and for the US it is inverted 32 % of the time. Inversions of the yield curve are obviously important for understanding very long term rates. The OU model with constant price of risk has four parameters to be estimated: m, k, α and q. The 3 month rates are much less sensitive to the risk parameter q than the 10 year rates. Therefore, by making the approximation that the 3 month is equivalent to the instantaneous process, we can estimate m, k and α from the 3 month rate time series alone (see below for the estimation of q). The parameters obtained in this way are shown in Table 2 and the values displayed are based on the maximum likelihood estimators derived in Ref. [23] (see Appendix C). To provide an appreciation for the robustness of the estimated parameters, we divide the time series into four blocks of equal size and evaluate the parameters m and k separately for each block, with the exception of the parameter α which is always estimated using the complete data set beacuse the time series in each block are too short for a reliable estimation of α. The quoted uncertainties in α are then simply the standard least square error computed when fitting an exponential autocorrelation function of the real interest time series. The maximum and minimum values obtained for each parameter is listed in Table 2. The variations are large, indicating some combination of autocorrelation, non-stationarity and heavy tails. We remind that m * implicitly depends on the market price of risk q (cf. Eq ( 14)) and, as reported in Table 2, the value of m * for a very short time window (τ = 3 months) can be approximated by m. To estimate q, we thus make use of the τ = 10 year real interest rate series as it is considered to be a large value of τ (much larger than 3 months). The averages of the 3 month and 10 year interest when considering Eq (20) give us two equations with two different values for m * where k and α are hold fixed. And, finally, the market price of risk q can be obtained with Eq (14). For this estimation we are also assuming in Eq (20) that r(t) ≃ m, i.e., that the mean historical interest rate is equal (or very close) to the current rate. Results are based on the annually sampled time series of three-month real interest rates.We use annual units. The Min and Max columns correspond to the minimum and the maximum values of the parameters obtained by splitting the time series into four blocks of equal length and estimating the parameters separately in each block (except α, see main text). For better estimates see Table 3.
The estimates obtained in this manner are slightly distorted not only because the 3 month bond is sampled annually but also because we have treated it as though it were an instantaneous rate. We can estimate the size of this bias by simulating the instantaneous process, which we approximate as having daily frequency. The simulation procedure is standard, we thus generate the instantaneous process r(t) by solving numerically the OU equation with parameters α, m, and k given by the empirical values taken from 3 month rates as given in Table 2 (see more details in Appendix C). We next create a surrogate time series for the 3 month real interest rate time series using Eq (20) with τ = 0.25 year and the empirical r = r(t) as the initial condition at each time t. We then mimic the procedure employed for the real data by estimating the parameters based on the surrogate 3 month rate, sampled at annual frequency. We can thus adjust the parameters of the instantaneous process to roughly match, on average, those observed for the real data (so that the estimated values based on the surrogate 3 month series match those of the observed 3 month series). The parameters with the bias corrected are given in Table 3 (see Appendix C for more details on the procedure we followed in order to correct the bias). The resulting shift in parameters is small, as can be seen by comparing Tables 2 and3. The main difference is in the parameter α, which sets the timescale for mean reversion; this changes by a little more than 10 %. Our numerical experiments indicate that the main source of the bias is the annual sampling and since parameter α sets the timescale for mean reversion it is not surprising that it is affected by this. We now make some comparisons of the simulated OU model with real data (see Table 4). Similarly to the procedure for the 3 month rates, we create simulated 10-year rates using Eq (20), so that rates are defined as (1/τ ) ln D(τ ) with τ = 10 year, where r = r(t) as the initial condition at each time t (recall that maturity time is t + τ ). A comparison shows that the standard deviation of the 3 month rates matches reasonably well, but the standard deviation of the simulated rates is much lower than that of the 10 year rates. We correct this by adding IID normally distributed noise to match the standard deviation of the simulated series with real data. We also neglect the 10 year smoothing of the 10 year inflation data. The simulated result has a lower correlation between 3 month and 10 year rates than real data; for the UK the correlations are 21 % (simulated) vs. 39 % real and for the US 24 % (simulated) vs. 59 % (real). However, the distributions for both 3 month and 10 years agree reasonably well. Figure 4  the simulated 3 month and 10 year interest rate time series, which should be compared to the real data shown in Figure 2. Not surprisingly, since the 3 month simulated rates are normally distributed, they lack the extreme values observed in the real data. The OU model does a good job of capturing the frequency of negative interest rates and yield curve inversions. Table 4 compares the frequency of negative interest rates for the real data and the simulation for both 3 month and 10 year rates. In Figure 5 we present a histogram of yield curve inversions for both the data and the model for the UK. The US histogram is qualitatively similar and we do not present it here. We use the difference between the 10 year real interest rate and the 3 month interest rate as our measure of inversion. The inversions of the data are somewhat more heavy tailed than those of the model, but the agreement is surprisingly good. The real UK yield curve is inverted roughly 50 % of the time and the simulated yield curves are inverted 46 % of the time. Similarly the real US yield curves is inverted 32 % of the time and the simulated yield curves are inverted 41 % of the time. With such a short time series as we have, it is difficult to estimate confidence intervals by methods such as bootstrapping. This is particularly true for α, where the time series properties of the data matter, so that one would need to do a block bootstrap and unfortunately there are not many blocks of sufficient length. However, assuming that the model is well-specified we can at least compute error estimates which are consistent with our model. The width of the resulting confidence intervals can be regarded as lower bounds on the width of the confidence intervals, and provide a perception for the magnitude of the estimating errors. We repeatedly simulate the instantaneous process r(t) using the parameters estimated from the data and generate 3 month and 10 year series as described above. In order to accurately mimic the constraints imposed by the data we sample the simulated series at an annual frequency. We then apply the estimation procedure described above to estimate the four parameters. Doing this 1 000 times allows us to compute the 5 % and 95 % quantiles for each parameter. The results are shown in Table 3. The estimated discount functions, together with their confidence intervals, are shown in Figure 6. The uncertainty intervals are estimated by repeatedly simulating the instantaneous, 3 month and 10 year processes as described above, applying the estimation procedure to the simulated data, and computing the discount function at each time interval. This is repeated 1 000 times to estimate the 5 % and 95 % quantiles. We are finally ready to present our key result.
The long term interest rate r ∞ is computed using Eq (22) based on the values in Table 3 and the final key result of this paper is presented in Table 5. The mean long run rate is r ∞ = 1.69 % for the UK and 2.21 % for the US.
Let us note that because of the scarcity of the data available uncertainties are substantial, with standard deviations of roughly 0.45 % in both cases. The real data are heavier tailed, but the agreement is otherwise reasonably good. Climate change and climate action are widely studied from a variety of perspectives [40]. However, researchers in financial economics have only recently turned their attention to climate change and climate finance is currently a quickly growing research field [41]. Finance academics participating in a recent survey have identified discount as a key topic to reduce climate risks [42]. With this paper, we have taken a muldisciplinary perspective from complex systems science and their related methods to study historical bond prices [43] with the aim to contribute to the need of new approaches to evaluate the climate action urgency [44,45]. More specifically, we have wanted to infer that long term discount rates is not just a trivial matter of extrapolating mean interest rates, but rather one must take several non-trivial factors into account. To begin with, because real interest rates are so often substantially negative, one should use a model that permits negative rates. This leads us to the Ornstein-Uhlenbeck model. While the presence of negative rates in this model may be viewed as a liability for describing nominal rates, for real rates this becomes a virtue. Another factor that should be considered is the market price of risk, which tends to raise longer term rates. Finally one should properly take into account the uncertainty and persistence of interest rates, which tends to lower the long-term discount rate. The use of the OU model accommodates all of these factors. When we estimate the OU model and compare it with real data, we see a good agreement on several essential properties, such as the frequencies of negative rates and yield curve inversions. Our results indicate that the long term interest rate used by Stern [3] is supported by historical data. His value of 1.4 % is less than a standard deviation below the estimated long term rate for the UK of 1.69 %, and just under two standard deviations of the US long term rate of 2.21 %. More recent estimates by Stern indeed pose an scenario with a value below 1 % [44]. In contrast, higher long run rates -see, for instance, [4,5]-seem not to be supported, as they are well above the 95 % confidence intervals of 2.63 for the UK and 3.07 for the US. Our estimates of 1.69 % (UK) and 2.2 % (US) are compatible with the rates recently estimated by ) which use data from UK housing markets during 2004-2013 and Singapore during 1995-2013 to estimate an annual discount rate of 2.6 % for payments more than 100 years in the future. Our results could potentially be improved on in several ways. One would be to acquire more data. This could include data for more countries -as we did recently but without considering market price of risk in a recent publication [27]-as well as longer term bonds or inflation-indexed bonds. Another possible improvement would be to extend the model to better capture the nonstationarity and/or heavy tailed behavior observed in the data [31]. In this regard we suspect that had we were able to take the observed nonstationarity [44] and/or heavy tails into account, mean values would have decreased because of the boosting of the uncertainty/persistence effect. This possibility is under present investigation. In a broader level, both Stern and Stiglitz have recently provided new methods and models which could be further studied from a complex systems science perspective [45]. Stern [44] has indeed put the accent in incorporating stronger multidisciplinary perspective to consider social and behavioral dimensions [46][47][48]. The same author [3] also mentions a definition of social discounts that might include inputs from social dilemmas and predefined behavioral experiments [49]. To conclude, we have demonstrated that historical data indicates that the long-term discount rate is probably not very large. While the error bars remain large, a value of 2 % or less seems plausible, corresponding to a present value of about 14 % for a payment received 100 years in the future. Other higher values do not seem to be fully consistent with the historical data. The need for immediate and substantial spending to combat climate change is thus sustained. We have seen in the main text that when rates are described by the OU process the joint characteristic function p(ω 1 , ω 2 , t|r 0 ) obeys the first-order partial differential equation (cf Eq (17))

with initial condition (t 0 = 0)

Due to the linearity of the OU process and the Gaussian character of the input noise, we may look for a solution of the initial value problem ( 27)- (28) in the form of a Gaussian density:

where A(ω 1 , t), B(ω 1 , t), and C(ω 1 , t) are unknown functions to be consistently determined. Substituting Eq (29) into Eq (27), identifying like powers in ω 2 and taking into account Eq (28), we find that these functions satisfy the following set of differential equations

Equation ( 30) is a first-order linear differential equation that can be readily solved giving

substituting this expression for A(ω 1 , t) into Eq (31) results in another first-order equation for B(ω 1 , t), whose solution reads

Finally, the direct integration of Eq (32) yields the expression for C(ω 1 , t)

From Eq (19) we see that the effective discount is given by the characteristic function, p(ω 1 , ω 2 , t|r 0 ), evaluated at the points ω 1 = -i and ω 2 = 0. Thus from Eqs ( 29) and ( 35) we obtain ln

which, with the change of notation regarding time as explained in the main text, agrees with Eq (20).
Knowing the exact form of the probability distribution of the joint process (x(t), r(t)) -through its joint characteristic function (29)-it is possible study in detail some interesting properties of the rate r(t), as we do next. As pointed out in the main text, a characteristic of the OU model is the possibility of attaining negative values. This probability is given by

where p(r, t|r 0 ) is the probability density function of the rate process. This is given by the marginal density

The probability for r(t) to be negative, Eq (37) then reads

where Erfc(z) is the complementary error function [38],

Note that as time increases (in fact starting from t > α -1 ) the probability ( 39) is well approximated by the stationary probability, defined as

In terms of the dimensionless normal level µ and the dimensionless volatility κ defined by

this probability reduces to

Let us now see the behavior of P (-) s

for the cases (i) µ < κ and (ii) µ > κ. (i) If the normal rate µ is smaller than rate's volatility κ, we use the series expansion [38] Erfc

Hence,

For µ/κ sufficiently small, this probability approaches 1/2. In other words, rates are positive or negative with almost equal probability. Note that this corresponds to the situation in which noise dominates over the mean. In the original units (cf. Eq (41)) this case corresponds to the noise intensity k being larger than mα 1/2 . (ii) When fluctuations around the normal level are smaller than the normal level itself, κ < µ, we use the asymptotic approximation [38] Erfc

and

Therefore, for mild fluctuations around the mean (that is, when k ≪ mα 1/2 ) the probability of negative rates is exponentially small. It is also interesting to know the probability that real rates r(t) are below the long-run rate r ∞ . This is given by

p(r, t|r 0 )dr. In the stationary regime, t → ∞, we have

where p(r) is the stationary PDF. For the OU model p(r) is obtained from Eq (38) after taking the limit t → ∞:

Substituting Eq (46) into Eq (45) we write

e -α(r+m) 2 /k 2 dr, or in terms of the complementary error function

Using the asymptotic estimates of the complementary error function discusse above, we see that this probability is exponentially small if (m -r ∞ ) → ∞ with √ α/k fixed, or if √ α/k → ∞ with a fixed differential of rates (m -r ∞ ). B Real and nominal rates. The market price of risk Recall that real rates are defined as the difference between nominal rates and inflation rates (cf. Eq (26)): r(t) = n(t) -i(t). We now discuss how to estimate n(t) and i(t) from empirical data . Let B(t|t + τ ) be the price at time t of a government bond maturing at time t + τ

It is also useful to define the "yield to maturity" y(t|τ ) as

The parameter m is easily estimated by noting that since the Wiener process has zero mean (stationary) rate E[r(t)] = m.

(61)

The parameters α and k can be estimated in terms of the correlation function

From Eqs (60) and (61) we write

e αt2 E[dW (t 1 )dW (t 2 )]. Taking into account that [39] E[dW (t 1 )dW (t 2 )] = δ(t 1 -t 2 )dt 1 dt 2 ,

where δ(•) is the Dirac delta function, and performing the integration on t 2 , we have

where Θ(•) is the Heaviside step function. In the evaluation of the integral we must take into account whether t > t ′ or t < t ′ . It is a simple matter to see that for either case the final result is

Let us incidentally note that Eq (62) proves that the correlation time of the OU process is given by α -1 . Indeed, the correlation time, τ c , of any stationary random process with correlation function K(τ ) is defined by the time integral of K(τ )/K(0). In our case

The empirical auto-correlation fit to an exponential (cf. Eq (62)) allows to estimate α. The parameter k can finally be obtained from the (empirical) standard deviation

We estimate these quantities for the three month interest rates using the maximum likelihood procedure given in [23]. C.2 Correcting for the bias using 3 month rates sampled at annual frequency

Our parameter estimation process that corrects for the bias introduced by using the 3 month rate as an approximation to the instantaneous rate, and by sampling it at annual frequency, has the following steps:

1. Estimate parameters using the historical 3 month and 10 year data as described in the main text. 2. Simulate the instantaneous process (which we approximate as a daily process) using the parameters inferred in step (1) to generate a simulated time series r(t) whose length matches that of the real data (roughly 100 years for the UK and 80 years for the US). 3. Construct simulated 3 month and 10 year time series based on Eq (20) with τ = 0.25 and τ = 10, using the time series for r(t) from step 1 as the initial condition for each time t.

4. Estimate m, k and α on the simulated 3 month series (sampled at annual frequency). 5. Repeat steps (2-4) for 1000 times and compute the average value of each parameter under the estimation process of step (4). This yields systematic shifts in the parameters relative to those estimated on the historical data, making it clear that the estimation process is biased. 6. Correct for this bias by adjusting the parameters of the instantaneous process by the magnitude of the average shift, so that the estimation process for the simulated 3 month bond times series roughly matches the values estimated from the historical series.
bidimensional process (x(t), r(t)) by the simple relation p(ω 2 , t|r 0 ) = p(ω 1 = 0, ω 2 , t|r 0 ). December 29, 2023

Measured in percent, as well as the 5 % and 95 % quantiles.December 29, 2023

December 29, 2023 27/27