Statistical analysis and stochastic interest rate modeling for valuing the future with implications in climate change mitigation

Statistical physics have been paying attention to economics and finance by providing new models and analyzing data available [1–3]. Most of the contributions investigate the nature of financial markets based on historical records, even its microstructure (see e.g. [4, 5]) or alternatively from a rather macroscopic and aggregated level (see e.g. [6–16]). However, there are still several issues in which an approach from physics can offer new perspectives and results. This is, for instance, the case of 'discounting' which in economics refers to weighting the future relative to the present [17]. Discounting constitutes the subject of this paper.

The choice of a discounting function has enormous consequences in many aspects of the global economy as, for instance, long-run environmental planning and, more specifically, climate action [18]. In a highly influential report on climate change commissioned by the UK government, Stern [19] uses a discounting rate of 1.4% while Nordhaus [20] argues for a discount rate of 4% and at other times [21] has advocated rates as high as 6%. Both estimates constitute a completely different point of view on how to address climate change. Indeed, while Stern's estimate would imply immediate spending, Nordhaus's figures indicate that immediate and strong action would be unnecessary. The choice of discount rate is, therefore, one of the biggest factors influencing the debate on the urgency of the response to climate change. Although Stern has been widely criticized for using such a low rate [20–25], our estimates are on average much closer to Stern than to Nordhauss and support more substantial immediate spending on climate actions. The Calderon report in July 2014 has also claimed that there is a false dilemma behind the choice between the economy growth and the environmental responsibility [26, 27].

Economists present a variety of reasons for discounting, including impatience, economic growth, and declining marginal utility; these are embedded in the Ramsey formula, which forms the basis for the standard approaches to discounting [28, 29]. Here we adopt the net present value approach, which treats the real interest rate as the measure of the trade-off between consumption today and consumption next year, without delving into the factors influencing the real interest rate.

It is often argued that, based on past trends in economic growth, future technologies will be so powerful compared with present technologies that it is more cost-effective to encourage economic growth—or solving other problems such as AIDS or malaria—than it is to take action against global warming now [25]. Analyses supporting this conclusion typically study discounting by working with an interest rate that is fixed over time, ignoring fluctuations about the average. This is mathematically convenient, but it is also dangerous: in this problem, as in many others, fluctuations play a decisive role.

A proper analysis takes fluctuations in the real interest rate, caused partly by fluctuations in growth, into account [30–32]. When the real interest rate r(t) varies randomly the discounting function becomes [33]

$$\begin{equation}D\left(t\right)=\mathbb{E}\left[\mathrm{exp}\left(-{\int }_{0}^{t}r\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }\right)\right],\end{equation} \tag{ 1 }$$

where the expectation $\mathbb{E}\left[\cdot \right]$ is an average over all possible interest rate paths. The fact that this is an average of exponentials, and not an exponential of an average, implies that the paths with the lowest interest rates dominate. This has been shown in several ways. Early papers analyzed an extreme case in which the annual real rate is unknown today, but starting tomorrow it will be fixed forever at one of a finite number of values [30, 31]. Other papers simulate stochastic interest rate processes out to some horizon, leaving aside the asymptotic behavior of real rates [32, 34–36].

The presence of fluctuations can dramatically alter the functional form of the discounting function. If real interest rates follow a geometric random walk, for example, the discounting function asymptotically may decay as a power law of the form D(t) = At^−1/2 [37] (see section 6). In contrast to the exponential function, this is not integrable on (0, ∞), underscoring how important the effect of persistent fluctuations can be. We have recently analyzed these issues by considering three of the most popular stochastic models for the dynamics of interest rates [33]: Ornstein–Uhlenbeck [38], Feller [39], and lognormal [40] processes, which are also very relevant in statistical physics. The Ornstein–Uhlenbeck (OU) model [38] is the only one that allows for negative rates r < 0 and its asymptotic expression has an exponential decay with a long-run rate r_∞ that differs from historical average interest rates by being substantially smaller, zero or eventually negative. We here want to go one step further and provide empirical estimates to such a discount based on historical data of interest rates from Argentina, Australia, Canada, Chile, Denmark, Germany, Italy, Japan, Netherlands, South Africa, Spain, Sweden, United Kingdom, and the United States. Such a diversity of countries, representing a variety of scenarios, allows us to better explore the intrinsic randomness of the real interest rates and how they lead to different costs of global economy planning such as climate action.

Real interest rates are nominal rates corrected by inflation so we need first of all to study nominal rates and inflation separately. The countries in our sample are: Argentina (ARG, 1864–1960), Australia (AUS, 1861–2012), Canada (CAN, 1913–2012), Chile (CHL, 1925–2012), Denmark (DNK, 1821–2012), Germany (DEU, 1820–2012), Italy (ITA, 1861–2012), Japan (JPN, 1921–2012), Netherlands (NLD, 1813–2012), South Africa (ZAF, 1920–2012), Spain (ESP, 1821–2012), Sweden (SWE, 1868–2012), United Kingdom (GBR, 1694–2012), and the United States (USA, 1820–2012). The details of each sample are reported in table 1.

Nominal rates can be obtained through the 10 years Government Bond Yield (see table 1 for further details). Following the standard procedure provided by the literature (see, for instance, [41]), we transform the annual rate β(t|T), where T = 10 years, into logarithmic rates, and denote the resulting nominal rates time series by

$$\begin{equation*}n\left(t\right)=\mathrm{ln}\left[1+\beta \left(t\vert T\right)\right].\end{equation*}$$

The inflation rate i(t) is estimated through the consumer price index (CPI) C(t) by

$$\begin{equation*}i\left(t\right)=\frac{1}{T}\sum _{j=0}^{T-1}\mathrm{ln}\left[1+C\left(t+j\right)\right],\end{equation*}$$

where T is chosen to be 10 years to be consistent with the 10 years nominal rate. We have, therefore, smoothed inflation rates with a 10 years forward moving average as this is again the standard procedure in these cases.

Finally, the real interest rate r(t) is defined by

$$\begin{equation}r\left(t\right)=n\left(t\right)-i\left(t\right).\end{equation} \tag{ 2 }$$

The recording frequency for each country is either annual or quarterly (see table 1). Some examples of the resulting real interest rates r(t) are plotted in figure 1.

Zoom In Zoom Out Reset image size

A striking feature observed in many epochs for all countries is that real interest rates frequently become negative, often by substantial amounts and for long periods of time (see figure 1 and table 2). This rules out most standard financial models, which assume that interest rates are always positive [41]. We thus focus our attention on one of the three most popular stochastic models and on the only one that allows for negative rates: the Ornstein–Uhlenbeck model [38], also known in the financial and economics literature as the Vasicek model [42] and which is also being used for modeling market volatility [6, 7, 9, 10]. The model can be written as [33]

$$\begin{equation}\mathrm{d}r\left(t\right)=-\alpha \left(r\left(t\right)-m\right)\mathrm{d}t+k\mathrm{d}w\left(t\right),\end{equation} \tag{ 3 }$$

where r(t) is the real interest rate and w(t) is a Wiener process, a Gaussian process with zero mean and unit variance. The parameter m is a mean value to which the process reverts and coincides with the long-term average of the process (3) :

$$\begin{equation}\mathbb{E}\left[r\left(t\right)\right]\simeq m.\end{equation} \tag{ 4 }$$

The parameter k is expressing the amplitude of the fluctuations and it is related to the variance which in the long-term limit reads

$$\begin{equation}\mathrm{Var}\left[r\left(t\right)\right]\simeq \frac{{k}^{2}}{2\alpha }.\end{equation} \tag{ 5 }$$

The parameter α is the strength of the reversion to the mean m. The autocorrelation function in its long-term limit is

$$\begin{equation}K\left(t-{t}^{\prime }\right)=\mathbb{E}\left[\left(r\left(t\right)-m\right)\left(r\left({t}^{\prime }\right)-m\right)\right]\simeq \frac{{k}^{2}}{2\alpha }{\mathrm{e}}^{-\alpha \vert t-{t}^{\prime }\vert },\end{equation} \tag{ 6 }$$

where α⁻¹ is the correlation time τ_c as can be seen from the definition

$$\begin{equation*}{\tau }_{\mathrm{c}}\equiv \frac{1}{K\left(0\right)}{\int }_{0}^{\infty }K\left(\tau \right)\mathrm{d}\tau =\frac{1}{\alpha }.\end{equation*}$$

Recall that the OU model may attain negative rates. Let us quantify this characteristic by evaluating the probability P(r < 0, t|r₀), for r(t) to be negative. In the long-term limit we denote this probability by ${P}_{\mathrm{s}}^{\left(-\right)}$, that is,

$$\begin{equation*}{P}_{\mathrm{s}}^{\left(-\right)}=\underset{t\to \infty }{\mathrm{lim}}P\left(r{< }0,t\vert {r}_{0}\right).\end{equation*}$$

For the OU model we have

$$\begin{equation}{P}_{\mathrm{s}}^{\left(-\right)}=\frac{1}{2}\mathrm{E}\mathrm{r}\mathrm{f}\mathrm{c}\left(\mu /\kappa \right),\end{equation} \tag{ 7 }$$

where Erfc(x) is the complementary error function expressed in terms of

$$\begin{equation}\mu =\frac{m}{\alpha },\quad \quad \kappa =\frac{k}{{\alpha }^{3/2}}.\end{equation} \tag{ 8 }$$

The dimensionless parameters μ and κ are related to the average m and the noise intensity k, respectively. As we will see later, these parameters provide a rather convenient way of describing important features about the discount function D(t). In figure 2, we represent equation (7) and show the different values that the function ${P}_{\mathrm{s}}^{\left(-\right)}$ can attain in terms of μ and κ.

Zoom In Zoom Out Reset image size

Using standard asymptotic expressions of Erfc(x) we can also get the behavior of ${P}_{\mathrm{s}}^{\left(-\right)}$ in the cases (a) μ < κ and (b) μ > κ.

• (a)  
  When the normal rate μ is smaller than the volatility of the rate κ, we can use the series expansion

  $$\begin{equation*}\mathrm{E}\mathrm{r}\mathrm{f}\mathrm{c}\left(z\right)=1-\frac{2}{\sqrt{\pi }}z+O\left({z}^{2}\right).\end{equation*}$$

  Hence,

  $$\begin{equation}{P}_{\mathrm{s}}^{\left(-\right)}=\frac{1}{2}-\frac{1}{\sqrt{\pi }}\left(\mu /\kappa \right)+O\left({\mu }^{2}/{\kappa }^{2}\right).\end{equation} \tag{ 9 }$$

  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 a rather stressed situation in which noise κ dominates over the mean value μ.
• (b)  
  When fluctuations around the normal level are smaller than the normal level itself, κ < μ, we can use the asymptotic approximation

  $$\begin{equation*}\mathrm{E}\mathrm{r}\mathrm{f}\mathrm{c}\left(z\right)\sim \frac{{\mathrm{e}}^{-{z}^{2}}}{\sqrt{\pi }z}\left[1+O\left(\frac{1}{{z}^{2}}\right)\right],\end{equation*}$$

  and

  $$\begin{equation}{P}_{\mathrm{s}}^{\left(-\right)}\sim \frac{1}{2\sqrt{\pi }}\left(\frac{\kappa }{\mu }\right){\mathrm{e}}^{-{\mu }^{2}/{\kappa }^{2}}.\end{equation} \tag{ 10 }$$

  Therefore, for mild fluctuations around the mean, the probability of negative rates is exponentially small.When κ = μ, the probability of negative rates is ${P}_{\mathrm{s}}^{\left(-\right)}=0.079$. Due to the ergodic character of the OU process [43], this means that when noise is balanced by the mean value (that is, κ = μ), one may expect to have negative real rates 7.9% of the time [33].

It is possible to derive the exact expression for the discount function D(t) defined in equation (1) in the case of the time-dependent OU model. As thoroughly described in reference [33], we write this expression in the form

$$\begin{equation}\mathrm{ln}D\left(t\right)=-\left(m-\frac{{k}^{2}}{2{\alpha }^{2}}\right)t+\frac{1}{\alpha }\left[m-{r}_{0}-\frac{{k}^{2}}{4{\alpha }^{2}}\left(3-{\mathrm{e}}^{-\alpha t}\right)\right]\left(1-{\mathrm{e}}^{-\alpha t}\right).\end{equation} \tag{ 11 }$$

The best way to study the discount rate is to work with the dimensionless time unit τ = αt, for afterward focusing on the long-term limit τ ≫ 1 since climate action is primarily interested in this asymptotic value. Thus, as τ → ∞, the exact expression (11) shows at once that the discount function of the OU model decays exponentially⁷

$$\begin{equation}D\left(t\right)\simeq {\mathrm{e}}^{-{r}_{\infty }t},\end{equation} \tag{ 12 }$$

where (see equation (8))

$$\begin{equation}{r}_{\infty }=m-{k}^{2}/2{\alpha }^{2}=\alpha \left(\mu -{\kappa }^{2}/2\right).\end{equation} \tag{ 13 }$$

We see from this expression that the long-run discount rate r_∞ is always lower than the average interest rate m, by an amount that depends on the dimensionless noise parameter κ. The long-run discount rate can therefore be much lower than the mean, and indeed can correspond to low interest rates that are rarely observed. This clearly illustrates the imprudence of assuming that the average real interest rate is the correct long-run discount rate.

The long-run behavior of the discount rate (13) depends on the two dimensionless parameters μ and κ (see equation (8)). The parameter space can be therefore divided into four regions, as shown in figure 3. In the region (1), where μ > κ²/2 (or equivalently m > k²/2α²) and μ > κ, the mean interest rate is large in comparison to the noise and negative rates are very infrequent. The long-run discounting function decays exponentially with rate r_∞ > 0. In the region (2), albeit small, the long-run discounting function still decays exponentially with rate r_∞ > 0 but negative rates are more frequent than 7.9%. Region (3) represents the most catastrophic situation since μ < κ²/2 and thus r_∞ < 0, meaning that the discount function D(t) increases exponentially and negative rates are rather frequent. Region (4) also shows r_∞ < 0 although, in this case, it is mostly because the noise component is very intense and not due to the presence of a relevant frequency of negative return events. Finally, at the boundary μ = κ²/2, the long-run interest rate r_∞ = 0 and the discount function is asymptotically constant.

Zoom In Zoom Out Reset image size

We now estimate the parameters m, k and α together with the dimensionless parameters μ and κ defined in equation (8). We perform such an estimation for each historical series (see table 1) by using a well-established maximum likelihood procedure for the OU model [41]. The resulting estimators $\hat{m}$, $\hat{\alpha }$, and $\hat{{k}^{2}}$ are listed in table 3 along with their standard deviation derived from formulas provided in reference [44]. Table 3 shows that the most inaccurate estimator is $\hat{\alpha }$, a not surprising fact since the estimation of α is quite a challenge in any Ornstein–Uhlenbeck process [44]. The last two columns in table 3 include the long-run interest rate estimator ${\hat{r}}_{\infty }$ and its error calculated through error propagation.

Table 3. Maximum likelihood estimation for the Ornstein–Uhlenbeck process described and the long-run interest rate. Countries have been reordered based on their estimated ${\hat{r}}_{\infty }$. $\hat{m}$ is the estimator of the mean real interest rate in 1/years. $\hat{\alpha }$ is the estimator related to the characteristic reversion time in 1/year. The squared root of the estimator of k² is the volatility of the process and k² is given in terms of 1/(year)³. These estimators are accompanied with the square root of the variance, σ's, of each estimator. ${\hat{r}}_{\infty }$ is the subsequent estimator of the long-run real interest rate in 1/year. Negative values of ${\hat{r}}_{\infty }$ mean the discount function is asymptotically increasing and its standard error is obtained through error propagation. The last two rows show separately the average over all countries, the stable countries with r_∞ > 0 and the unstable countries with r_∞ < 0. In all three rows standard error provided corresponds to the standard deviation of the ${\hat{r}}_{\infty }$ for the different countries.

Country	$\hat{m}$	${\sigma }_{\hat{m}}$	$\hat{\alpha }$	${\sigma }_{\hat{\alpha }}$	$\hat{{k}^{2}}$	${\sigma }_{\hat{{k}^{2}}}$	${\hat{r}}_{\infty }$	${\sigma }_{{\hat{r}}_{\infty }}$
Germany	−0.0945	0.6695	0.0071	0.0089	41.72 × 10−4	2.19 × 10−4	−40.94	2.28
Chile	−0.0579	0.3146	0.0201	0.0227	31.07 × 10−4	2.49 × 10−4	−3.917	0.442
Japan	0.0502	0.2468	0.0053	0.0114	13.96 × 10−5	1.09 × 10−5	−2.431	0.314
Italy	0.0197	0.1595	0.0056	0.0089	11.46 × 10−5	0.68 × 10−5	−1.778	0.192
Spain	0.0671	0.0692	0.0167	0.0137	23.71 × 10−5	1.26 × 10−5	−0.3578	0.0728
Argentina	0.0315	0.0709	0.0228	0.0231	22.40 × 10−5	1.71 × 10−5	−0.1831	0.0727
Australia	0.0397	0.0450	0.0089	0.0112	2.23 × 10−5	0.13 × 10−5	−0.1029	0.0458
South Africa	0.0269	0.0472	0.0154	0.0193	4.35 × 10−5	0.34 × 10−5	−0.0649	0.0477
Canada	0.0266	0.0391	0.0142	0.0178	2.75 × 10−5	0.21 × 10−5	−0.0415	0.0394
Denmark	0.0410	0.0259	0.0161	0.0133	3.15 × 10−5	0.17 × 10−5	−0.0197	0.0261
Sweden	0.0279	0.0166	0.0676	0.0317	16.92 × 10−5	2.06 × 10−5	0.0095	0.0167
USA	0.0319	0.0123	0.0603	0.0257	10.03 × 10−5	1.05 × 10−5	0.0181	0.0124
UK	0.0342	0.0062	0.1635	0.0326	31.37 × 10−5	2.53 × 10−5	0.0283	0.0062
Netherlands	0.0599	0.0078	0.1648	0.0550	17.97 × 10−5	2.43 × 10−5	0.0566	0.0078
All countries	0.0217	0.1236	0.0420	0.0211	63.45 × 10−5	4.31 × 10−5	−3.552	0.255
Stable	0.0385	0.0107	0.1140	0.0362	19.07 × 10−5	2.02 × 10−5	0.0281	0.0108
Unstable	0.0150	0.1686	0.0132	0.0150	81.20 × 10−5	5.23 × 10−5	−4.984	0.353

We can also observe the position $\left(\hat{\kappa },\hat{\mu }\right)$ of each country in figure 3 by considering the results presented in table 4. In any case these results need to be understood as a first-order approximation since the errors behind the estimators (which are evaluated through error propagation) are significant (see table 4). Only four countries show a positive long-run rate, r_∞ > 0, and all of them inside, or very close, to the region defined by μ < κ in which rates are frequently negative. The other ten countries show less stable behavior and are all of them in the exponentially increasing region (region 3), which implies they have long-run negative rates, and are widely scattered. In two cases (Germany and Chile) the average rate m (and its dimensionless version μ) is negative due to at least one period of runaway inflation while two others (Japan and Italy) still have a long-run negative rate r_∞ mostly due to a very small strength of the reversion to the mean given by the parameter α (see equation (3)). These four countries are not plotted in figure 3 because they are out the range of μ and/or κ axis.

Table 4. Dimensionless mean interest rate and fluctuation amplitude for all countries. The dimensionless mean interest rate estimator $\hat{\mu }$ is accompanied with its error obtained through error propagation (see equation (8)) and by considering the parameters estimated and provided in table 3. The dimensionless fluctuation amplitude estimator $\hat{\kappa }$ is accompanied with its standard error obtained through error propagation (see equation (8)) and by considering the parameters estimated and provided in table 3.

Country	$\hat{\mu }$	${\sigma }_{\hat{\mu }}$	$\hat{\kappa }$	${\sigma }_{\hat{\kappa }}$
Germany	−13.22	95.11	106.92	198.75
Chile	−2.89	16.01	19.61	33.26
Japan	9.46	50.79	30.59	98.70
Italy	3.49	28.79	25.23	59.94
Spain	4.02	5.30	7.13	8.79
Argentina	1.38	3.40	4.34	6.58
Australia	4.48	7.61	5.67	10.77
South Africa	1.75	3.77	3.45	6.50
Canada	1.88	3.62	3.10	5.83
Denmark	2.55	2.65	2.75	3.41
Sweden	0.41	0.31	0.74	0.52
USA	0.53	0.30	0.68	0.43
UK	0.21	0.06	0.27	0.08
Netherlands	0.36	0.13	0.20	0.10
All countries	1.03	15.56	15.05	30.99
Stable	0.39	0.20	0.47	0.28
Unstable	1.29	21.71	20.89	43.25

Also note that all fourteen countries but one (Netherlands) are below the identity line, μ = κ, in figure 3 which indicates that negative real interest rates are common (even in the stable countries they occur 20% of the time). It is also worth to mention that only one is above Nordhaus's 4% discounting rate [20] (5.7%, Netherlands) and only two more countries are above the more pessimistic discounting rate (1.4%) provided by Stern [19] (1.8% and 2.8% from USA and United Kingdom, respectively). And more generally, it is important to notice that r_∞ is very much smaller than m in most of the cases. All these statements are robust even when considering values of the estimators with shifts of the size of its standard error (see table 3).

The characteristic (correlation) time (τ_c = 1/α) for each country appears to be very different (see table 3). Some countries must spend more than a century to achieve a stationary level and thus finally attain the long-run discount rate r_∞. Furthermore, this time horizon might be even larger than the time interval we must consider to make a response, from an economic point of view, to any climate change catastrophe. For this reason, it is interesting to investigate how the discount rate defined as −ln(D(t))/t changes over time (see equation (11)).

Figure 4 shows the discount rates for all countries as a function of time by considering initial rate r₀ = 1% which clearly illustrates the dramatic differences between countries. In this way we divide the fourteen countries into roughly four groups. There are two countries (DEU, CHL) that show a very fast and very negative rate. There is a second group still having a monotonic behavior but with a much slower trend to raise negative discount rates (JPN, ITA, ESP and ARG). Non-monotonic behavior is indeed observed in a third group (AUS, ZAF, CAN, DEN). This group is of special interest since it shows how the rates might first grow by finally becoming negative after 20 or 30 years. Stable countries represented in the first inset on the left of figure 4 also show that the asymptotic rate r_∞ is raised very slowly being the country with the highest rate (NED) the one that needs more than a century to attain the stationary level. Figure 5 selects four countries (USA, DEU, ZAF and ESP), one from each of the groups mentioned above, to observe the impact of uncertainty as a function of time. For different values of time, the discount function includes a shadow in gray limited by maximum and minimum values when taking into consideration the standard error of each of the estimators. In all four countries and at any time, maximum discount rate value is always below 2.2%. The inclusion of the statistical uncertainty reinforces the robustness of our results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Let us finally note that these results are in contrast to other treatments of fluctuating rates which assume that short term rates are positive and predict that the decrease in the discounting rate occurs over a much longer timescale, usually measured in hundreds or thousands of years [30, 32, 34–37, 45].

As mentioned above the Ornstein–Uhlenbeck model is the only one among the three most classic models allowing for negative rates. This is the reason why we have excluded both the Feller and the lognormal models from our analysis. Let us nonetheless briefly study what modifications should be carried out in order to use these positive rate models in our analysis.

6.1. The shifted Feller model

The Feller process [39] (see also [43]) has a very similar structure than the Ornstein–Uhlenbeck process except that the noise component depends on the interest rate. The process also has a mean reverting force that makes the process have an autocorrelation function that decays in an exponential manner whose characteristic time scale is 1/α_F. Let us however remind that Feller does not allow for negative rates and these are clearly present in our empirical data. Therefore, to consider the Feller process for estimating the long-run discount rate r_∞ would require to redefine the Feller model that reads

$$\begin{equation}\mathrm{d}y\left(t\right)=-{\alpha }_{\mathrm{F}}\left(y\left(t\right)-{m}_{\mathrm{F}}\right)\mathrm{d}t+{k}_{\mathrm{F}}\sqrt{y\left(t\right)}\mathrm{d}w\left(t\right),\end{equation} \tag{ 14 }$$

where

$$\begin{equation}y=r-{r}_{\mathrm{min}}\end{equation} \tag{ 15 }$$

and r_min < 0 is the minimum value observed in the time series. The estimation through maximum likelihood procedure and its error analysis is also possible [44]. Table 5 compares the Ornstein–Uhlenbeck and the shifted Feller models which have very similar mathematical expressions for estimating $\hat{m}$, $\hat{\alpha }$ and ${\hat{k}}^{2}$ parameters. The discount function now reads (see equations (1)) and (15))

$$\begin{equation*}D\left(t\right)=\mathbb{E}\left[\mathrm{exp}\left(-{\int }_{0}^{t}\left({r}_{\mathrm{min}}+y\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }\right)\right)\right]=\mathrm{exp}\left(-{r}_{\mathrm{min}}t\right)\mathbb{E}\left[\mathrm{exp}\left(-{\int }_{0}^{t}y\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }\right)\right].\end{equation*}$$

The asymptotic value of the remaining average shows an exponential decay [33]

$$\begin{equation*}\mathbb{E}\left[\mathrm{exp}\left(-{\int }_{0}^{t}y\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }\right)\right]\simeq {\mathrm{e}}^{-{y}_{\infty }^{\mathrm{F}}t},\end{equation*}$$

whose long-run discount rate [33] is

$$\begin{equation*}{y}_{\infty }^{\mathrm{F}}=\frac{2{m}_{\mathrm{F}}}{1+\sqrt{1+2{k}_{\mathrm{F}}^{2}/{\alpha }_{\mathrm{F}}^{2}}},\end{equation*}$$

so that

$$\begin{equation*}D\left(t\right)={\mathrm{e}}^{-\left({r}_{\mathrm{min}}+{y}_{\infty }^{\mathrm{F}}\right)t}={\mathrm{e}}^{-{r}_{\infty }^{\mathrm{F}}t}.\end{equation*}$$

We observe that (as in the Ornstein–Uhlenbeck process) the log-run rate is smaller than the average rate, ${r}_{\infty }^{\mathrm{F}}{< }m$. However, the shifted Feller process leads us to obtain a slightly larger estimation but within the statistical error range (3.49% versus 1.81%, see table 5). The value is similar than Nordhaus's 4% discounting rate if one considers the statistical error and in any case lower than 6% [20, 21].

Table 5. Maximum likelihood estimation for the three different models described and the long-run interest rate for the case of United States. These estimators are taking years as a basic time unit and they are all accompanied with the square root of the variance, σ's, of each estimator. ${\hat{r}}_{\infty }$ is the subsequent estimator of the long-run real interest rate in 1/year. For the Feller and lognormal cases we have provided a modified version of the model (see equations (14) and (16)). The shifted Feller and lognormal models takes r_min = −0.0415, which is its minimum value in the historical time series, and the estimated ${\hat{r}}_{\infty }$ is corrected by adding r_min but this has been impossible to be done in the lognormal case since the asymptotic value goes to a constant.

Ornstein–Uhlenbeck	$\hat{m}$	${\sigma }_{\hat{m}}$	$\hat{\alpha }$	${\sigma }_{\hat{\alpha }}$	${\hat{k}}^{2}$	${\sigma }_{{\hat{k}}^{2}}$	${\hat{r}}_{\infty }$	${\sigma }_{{\hat{r}}_{\infty }}$
Shifted Feller	0.0319	0.0123	0.0603	0.0257	10.03 × 10−5	1.05 × 10−5	0.0181	0.0124
	${\hat{m}}_{\mathrm{F}}$	${\sigma }_{{\hat{m}}_{\mathrm{F}}}$	${\hat{\alpha }}_{\mathrm{F}}$	${\sigma }_{{\hat{\alpha }}_{\mathrm{F}}}$	${\hat{k}}_{\mathrm{F}}^{2}$	${\sigma }_{{\hat{{k}^{2}}}_{\mathrm{F}}}$	${\hat{r}}_{\infty }^{\mathrm{F}}$	${\sigma }_{{\hat{r}}_{\infty }^{\mathrm{F}}}$
Shifted lognormal	0.0864	0.0041	0.0599	0.0057	12.56 × 10−5	1.31 × 10−5	0.0349	0.0072
	${\hat{m}}_{\mathrm{L}}$	${\sigma }_{{\hat{m}}_{\mathrm{L}}}$	${\hat{k}}_{\mathrm{L}}^{2}$	${\sigma }_{{\hat{k}}_{\mathrm{L}}^{2}}$	${\hat{m}}_{\mathrm{L}}-{\hat{k}}_{\mathrm{L}}^{2}/2$	${\sigma }_{{\hat{m}}_{\mathrm{L}}-{\hat{k}}_{\mathrm{L}}^{2}/2}$	Asymp
	0.0130	0.0163	0.0309	0.0066	−0.0024	0.0130	Constant

6.2. The shifted lognormal model

Another alternative to still allow for negative rates is to consider a modified version of the lognormal process by considering new variable y = r − r_min (where r_min < 0 is the minimum value observed in the time series) and the following stochastic dynamics:

$$\begin{equation}\mathrm{d}y\left(t\right)={m}_{\mathrm{L}}y\left(t\right)\mathrm{d}t+{k}_{\mathrm{L}}y\left(t\right)\mathrm{d}w\left(t\right),\end{equation} \tag{ 16 }$$

whose long-run discount function can lead to three different asymptotic expressions [33]:

$$\begin{equation*}\mathbb{E}\left[\mathrm{exp}\left(-{\int }_{0}^{t}y\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }\right)\right]\sim \left\{\begin{array}{cc}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t}\hfill & \quad {m}_{\mathrm{L}}{< }{k}_{\mathrm{L}}^{2}/2,\hfill \\ {\mathrm{e}}^{-{y}_{\infty }^{L}t}\hfill & {m}_{\mathrm{L}}{ >}{k}_{\mathrm{L}}^{2}/2.\hfill \\ {t}^{-1/2}\hfill & {m}_{\mathrm{L}}={k}_{\mathrm{L}}^{2}/2,\hfill \end{array}\right.\end{equation*}$$

For the exponential case the long-run discount rate reads

$$\begin{equation*}{y}_{\infty }^{\mathrm{L}}=\frac{{m}_{\mathrm{L}}-{k}_{\mathrm{L}}^{2}/2}{\psi \left(2{m}_{\mathrm{L}}/{k}_{\mathrm{L}}^{2}\right)+1/\left(2{m}_{\mathrm{L}}/{k}_{\mathrm{L}}^{2}-1\right)},\end{equation*}$$

where ψ(⋅) is the digamma function. The lognormal process does not show any reversion trend to a certain level and its average grows (or decreases) in an exponential manner

$$\begin{equation*}\mathbb{E}\left[r\left(t\right)\vert r\left(0\right)={r}_{0}\right]=\left({r}_{0}+{r}_{\mathrm{min}}\right){\mathrm{e}}^{{m}_{\mathrm{L}}t}-{r}_{\mathrm{min}},\end{equation*}$$

a result that it is in contradiction with the times series provided in figure 1. We can however also estimate the parameters via maximum likelihood procedures. The results are again provided in table 5 and they show us that the asymptotic discount falls into the constant case since ${m}_{\mathrm{L}}{< }{k}_{\mathrm{L}}^{2}/2$ although the error analysis warn us that we cannot discard the exponential case (being not greater 4%–5%) nor the hyperbolic slow decay.

6.3. Extending the Ornstein–Uhlenbeck process

One can argue that the results presented can change under different historical conditions or periods. To exemplify this issue, we have also estimated these values in the case of Germany once the World War II was over (from March 1946). Parameters are in that case μ = 0.62, κ = 0.32 with now a positive long-run rate r_∞ = 3.4% which is in any case smaller than Nordhaus estimates for valuing climate action [21]. Germany certainly is a quite volatile situation challenging the model which assumes constant (i.e., stationary) parameters. A possible way out is to extend the model with an additional dimension under the form of a 'matrioshka doll' by considering m as a stochastic process following an additional Ornstein–Uhlenbeck process [9]⁸

$$\begin{equation}\begin{array}{cc}\hfill \mathrm{d}r& =-\alpha \left(r-m\right)\mathrm{d}t+k\mathrm{d}w\left(t\right)\hfill \\ \mathrm{d}m& =-{\alpha }_{0}\left(m-{m}_{0}\right)\mathrm{d}t+{k}_{0}\mathrm{d}{w}_{0}\left(t\right),\hfill \end{array}\end{equation} \tag{ 17 }$$

where the Wiener processes w(t) and w₀(t) are both zero mean, have unit variance and are independent from each other implying that $\mathbb{E}\left[\mathrm{d}w\left(t\right)\mathrm{d}{w}_{0}\left(t\right)\right]=0$. We also assume that α > α₀ > 0 thus showing a slower mean reverting force for the subordinated process m₀ than for m. A similar extension has been used in other financial contexts to model stochastic volatility [9] by adding a longer mean reversion process which allows for a slow decaying memory for the volatility process while still preserving a much shorter memory for the so-called leverage effect [7, 10] (see also reference [16] for another setting that could represent an alternative approach to the extended Ornstein–Uhlenbeck model given by equation (17)). In the long-run, the process reads [9]:

$$\begin{align}\hfill m\left(t\right)& ={m}_{0}+{k}_{0}{\int }_{-\infty }^{t}{\mathrm{e}}^{-{\alpha }_{0}\left(t-{t}^{\prime }\right)}\mathrm{d}{w}_{0}\left(t\right)\hfill \\ \hfill r\left(t\right)& ={m}_{0}+k{\int }_{-\infty }^{t}{\mathrm{e}}^{-\alpha \left(t-{t}^{\prime }\right)}\mathrm{d}w\left(t\right)+\frac{{k}_{0}}{\alpha -{\alpha }_{0}}{\int }_{-\infty }^{t}\left({\mathrm{e}}^{-{\alpha }_{0}\left(t-{t}^{\prime }\right)}-{\mathrm{e}}^{-\alpha \left(t-{t}^{\prime }\right)}\right)\mathrm{d}{w}_{0}\left(t\right).\hfill \end{align} \tag{ 18 }$$

We can easily see that this extended process shows the same average, $\mathbb{E}\left[r\left(t\right)\right]={m}_{0}$, than the simpler OU version but with greater variance (see equations (4) and (5), respectively)

$$\begin{equation}\mathrm{Var}\left[r\left(t\right)\right]={k}^{2}/2\alpha +{k}_{0}^{2}/2{\alpha }_{0}.\end{equation} \tag{ 19 }$$

The autocorrelation function now reads [9] (see equation (6)):

$$\begin{align}\hfill K\left(t-{t}^{\prime }\right)& =\mathbb{E}\left[\left(r\left(t\right)-{m}_{0}\right)\left(r\left({t}^{\prime }\right)-{m}_{0}\right)\right]\hfill \\ \hfill & =\left(\frac{{k}^{2}}{2\alpha }-\frac{{k}_{0}^{2}\alpha }{2\left({\alpha }^{2}-{\alpha }_{0}^{2}\right)}\right){\mathrm{e}}^{-\alpha \vert t-{t}^{\prime }\vert }+\frac{{k}_{0}^{2}{\alpha }^{2}}{2\left({\alpha }^{2}-{\alpha }_{0}^{2}\right){\alpha }_{0}}{\mathrm{e}}^{-{\alpha }_{0}\vert t-{t}^{\prime }\vert },\hfill \end{align}$$

where the first term with an exponential decay with 1/α would dominate for short time difference |t − t'|. In the opposite situation, for longer time difference, the second exponential decay expressed in terms of 1/α₀ would dominate. The extended process now have five parameters (α, k, α₀, k₀, and m₀), while basic OU process had only three (α, k, and m).

We can finally look at the effects on the asymptotic discount. It will shown in a future publication that the process has

$$\begin{equation}{r}_{\infty }^{\mathrm{e}\mathrm{x}\mathrm{t}}={m}_{0}-\frac{1}{2}\left(\frac{{k}^{2}}{2{\alpha }^{2}}+\frac{{k}_{0}^{2}}{2{\alpha }_{0}^{2}}\right).\end{equation} \tag{ 20 }$$

The result brings rates which are even lower than the one provided by the maximum likelihood estimation procedure in the simple Ornstein–Uhlenbeck process. As a simple exercise we can estimate a combination of k₀ and α₀ with the historical variance of the whole process (see equation (19)). To estimate α₀ is not that simple since our historical data sets are too short and its estimation becomes too noisy. However, jointly with the values already obtained for Ornstein–Uhlenbeck maximum likelihood estimation for m (now equivalent to m₀), k, and α, it is possible to observe the effects for different values of α₀:

$$\begin{equation*}{r}_{\infty }^{\mathrm{e}\mathrm{x}\mathrm{t}}={m}_{0}-\frac{1}{2}\left[\frac{{k}^{2}}{2{\alpha }^{2}}+\left(\mathrm{Var}\left(r\left(t\right)\right)-\frac{{k}^{2}}{2\alpha }\right)\frac{1}{{\alpha }_{0}}\right].\end{equation*}$$

In this case we can see that the long-run rate r_∞ for the United States is practically zero when α₀ = α/20.

Our empirical analysis proves that real interest rates are often negative—roughly a quarter of the time—which implies that one must use a discount model that is compatible with this property. For this purpose we have proposed the Ornstein–Uhlenbeck model which has the additional advantage that it can be solved analytically in a relatively simple way. This model facilitates the understanding of why the long-run discount rate can be so low. A first reason is that real interest rates are themselves typically low. As we have showed the average over all countries surveyed is negative, and even the average over stable countries (those with a positive long-run rate, r_∞ > 0) is 2.8%. A second reason is that the fluctuating part on the right hand side of equation (13), which depends both on the noise intensity k and the persistence term 1/α, typically lowers rates for the stable countries by about 7%. In some cases, such as Spain, the effect is much more dramatic: even though the mean short term rate has the high value of m = 6.7%, the long-term discounting rate is r_∞ = −36% which would imply a great increasing discount. The estimation is being done with a maximum likelihood procedure that includes an error analysis that demonstrates the robustness of the results obtained.

Our analysis here makes several simplifications such as ignoring non-stationarity. We have here partially address this issue by extending the Ornstein–Uhlenbeck which allows for slower fluctuations in the normal level and resulting even in lower long-run discount rates⁹. Correlations between the environment and the economy have also being ignored but, in any case, despite the variety of results, the long-run discount rate is always smaller than Nordhaus estimates by other methods as we have exemplified with the German case [21]. We have also not considered the market price of risk [42, 46], in other words, we have assumed that markets are risk neutral and the average in equation (1) defining the discount function, is evaluated using the empirical probability measure without any risk adjustment [47, 48]. These issues are under present investigation and some results are expected soon [49].

In any case the methods that we have introduced here provide a foundation on which to incorporate more realistic assumptions. We do not mean to imply that it is realistic to actually use the increasing discounting functions that occur for countries with less stable interest rate processes. There is some validity to treating hyper-inflation as an aberration—when it occurs government bonds are widely abandoned in favor of more stable carriers of wealth such as land and gold, and as a result under such circumstances the difference between nominal interest and inflation may underestimate the actual real rate of interest.

Nonetheless, real interest rates are closely related to economic growth, and economic downturns are a reality. The great depression lasted for 15 years, and the fall of Rome triggered a depression in western Europe that lasted almost a thousand years. In light of our results here, arguments that we should wait to act on global warming because future economic growth will easily solve the problem should be viewed with extreme skepticism. Our analysis clearly supports Stern over Nordhaus. When we plan for the future we should always bear in mind that sustained economic downturns may visit us again, as they had in the past.

Effective responses to this multifaceted problem have been slow to develop, in large part because many experts have not only underestimated its impact, but also overlooked the underlying institutional structure, organizational power and financial roots [50, 51]. A growing body of sophisticated research is currently emerging with a large set of multidisciplinary strategies that wants to exploit socioeconomic tipping points (as in any complex dynamical system) to magnify the impact of each political intervention [52] and also integrate science-policy perspectives with public awareness, citizen-led research and citizen science practices (see for instance [53, 54]). In all cases the final purpose is to better respond to global challenges such as climate action in a near future, sooner rather than later.

We thank to an anonymous referee for the several comments that have allowed to improve the final version of the paper. We also would like to thank National Science Foundation grant 0624351 (JG) and the Institute for New Economic Thinking (J D F). This work was supported by MINECO (Spain) FIS2013-47532-C3-2P (J M, M M and J P), FIS2016-78904-C3-2P (J M, M M and J P); by Generalitat de Catalunya (Spain) through Complexity Lab Barcelona (contract no. 2017 SGR 1064, J M, M M and J P)