Valuing the Future and Discounting in Random Environments: A Review

Abstract

We address the process of discounting in random environments, which allows valuation of the future in economic terms. We review several approaches to the problem regarding different well-established stochastic market dynamics in the continuous-time context and include the Feynman–Kac approach. We also review the relation between bond-pricing theory and discounting and introduce both the market price of risk and the risk neutral measure from an intuitive point of view devoid of excessive formalism. We provide the discount for each economic model and discuss their key results. We finally present a summary of our previous empirical studies for several countries on the long-run discount problem.

1. Introduction

The introduction around three decades ago of the view and methods of statistical physics into economics and finance signaled the appearance of a new interdisciplinary aspect of physics, which is sometimes called “econophysics” [1,2,3]. The fact that financial prices are random with sudden and uncontrollable ups and downs has been long-known; however, the first step towards a systematic mathematical analysis of price randomness was taken by Bachelier in 1900, who proposed a model for the market dynamics in which the prices follow ordinary Brownian motion [4].

However, Bachelier’s model is not completely satisfactory because, in such a representation, prices can be either positive or negative, contradicting one of the most fundamental tenets of economics, the “principle of limited liability”, which affirms that prices cannot attain negative values. This limitation in Bachelier’s model was remedied more than six decades later by Osborne [5] by assuming the geometric Brownian motion where prices are described by the exponential of the ordinary Brownian motion and, hence, they can never attain negative values.

Let us denote by $S\left(t\right)$ a speculative price (or an economic index) at time t. In the continuous time framework, the geometric Brownian motion assumes that

$$S\left(t\right)=S_0{e}^{x\left(t\right)},$$

(1)

where $S_0=S\left(t_0\right)$ is the price at some initial time $t_0$ and $x\left(t\right)$, the so-called return, is described by the ordinary Brownian motion, that is to say, by the stochastic differential equation

$$dx\left(t\right)=mdt+\sigma dW\left(t\right),$$

(2)

where $W\left(t\right)$ is the standard Wiener process with zero mean and unit variance. Note that the return is assumed to be a diffusion process with constant drift m and diffusion coefficient $\sigma$. In this model, the return is a Gaussian process with the mean and variance given, respectively, by m and ${\sigma }^2$. The price is, hence, a log-normal process, and the geometric Brownian motion is also called the log-normal model.

Despite the log-normal model being used in countless financial applications, it has certain limitations [6], which have given rise to several generalizations. One of them assumes that the return is a more complex diffusion process obeying a stochastic equation of the form

$$dx\left(t\right)=f(x,t)dt+g(x,t)dW\left(t\right),$$

(3)

which is interpreted in the sense of Itô. In this case, returns, and hence prices, are driven by an external “force” and by multiplicative noise, which, in the most general case, depend explicitly on the time and on the return level. Function $f(x,t)$ drives prices, and function $g(x,t)$ modulates the intensity of the fluctuations around the deterministic motion set by $f(x,t)$. In any case, and regardless of the values taken by $x\left(t\right)$, the prices given by Equation (1) are always nonnegative, thus, keeping the principle of limited liability.

Another significant shortcoming of the geometric Brownian motion model is the absence of both “fat tails” and skewness in the distribution of log-prices (i.e., returns). Indeed, empirical distributions of log-prices not only show fat tails—meaning that extreme losses and profits have a higher probability than those of the log-normal model—but also an asymmetric shape in the sense that losses are usually more probable than profits [6]. In order to address these and other problems, intense research has been prompted both in mathematics and physics, which, among others, may involve the use of the Lévy process as driving noise (instead of the Wiener process). One of the most popular alternative financial models—proposed by Mandelbrot [7] and Fama [8] in the early 1960s—is provided by substituting the Wiener process by the Lévy process, which can take into account the appearance of fat tails in the probability distribution of prices, a widely accepted empirical fact [6]. A major inconvenience of the non-Poissonian Lévy jump processes is, however, their lack of finite moments apart from (at most) the first one, which does not seem to be case in empirical data [6]. or models in which the variance ${\sigma }^2$ (or the noise intensity g) is a random process, such as in the so-called “stochastic volatility models” [9,10,11].

In economics and finance, one of the most consequential developments is that of “discounting”, which essentially attempts to answer the crucial question of what the price will be in the future. In other words, discounting weighs the future relative to the present. Traditionally, the weighting procedure has been performed through a decreasing exponential. Thus, under a constant interest rate r, continuously compounded, a dollar invested today at time $t=0$ yields $e^{rt}$ dollars at time t [12]; hence, one dollar in any future time t is worth $e^{-rt}$ today. This statement is true under constant and fixed rates; however, in real life, rates are random, and this uncertainty makes it completely unrealistic to represent rates by constant quantities or even by deterministic functions of time, and, as a consequence, random models for rates must be addressed.

The problem of discounting is widely known in finance where it has been thoroughly studied closely related to bond pricing particularly over short periods of time [13]. Discounting, particularly in the long run, is of importance not only in the context of finance but to many other aspects of the global economy. For instance, we may consider the long-term environmental planning, which is certainly sensitive in relation to climate action. Thus, in an oversimplified way, an environmental problem, which costs X to fix at a time t is worth an investment of $e^{-rt}X$ today.

This analysis assumes that the interest rate remains constant between today and the distant future t. The rate r becomes a key magnitude to decide whether it becomes more beneficial to take action today with a significant investment or whether the discount gives negligible value to today’s investment. The choice of discount rate is perhaps the biggest factor influencing the debate on the urgency of the response to global warming as it relates today’s investments with potential climate-related losses in future [14]. No wonder that, in recent years, obtaining a long-run discount rate valid for decades ahead has been the object of intense controversy.

Thus, Nicholas Stern, in an influential report commissioned by the UK government, advocated for a long-run discount of 1.4% [15], which on a 100-year horizon implies a present value of 25% (meaning that the future is worth 25% as much as the present). On the other hand, Willian Nordhaus proposed a discount rate of 4% [16] (implying a present value of 2%) and even the higher value of 6% [17], which implies a present value of 0.3%. Stern has been widely criticized for using such a low rate [18,19,20], and the question is far from being settled.

Economists present a variety of reasons for discounting, particularly for environmental problems in the long run. These reasons include, among others, ethical considerations [21,22], impatience, economic growth [23] and arguments based on the maximization of utility functions that are mostly chosen for mathematical convenience [24], all of them ingrained in a phenomenological expression called the Ramsey formula [25], which constitutes the standard approach to discounting in the economics literature (particularly in the long-run) [14].

From an empirical point of view, any practical economist involved in environmental debates might consider the average of historical interest rates, which occurred in the last 200 hundred years to estimate the forward discount rate (which is 2.7% in the less unstable countries [26]) or take the average of Wall Street forward looking models, with price bonds of maturity as long as 30 years. Unfortunately, due to historical fluctuations of real interest rates, the appropriate rate is considerably below such averages [26].

In econophysics, the problem of discounting, despite its relevance, is virtually unknown. The main purpose of this paper is to offer a survey of the problem devoid of excessive formalism and abstraction as well as to review some of our recent work on the problem [26,27,28,29,30]

2. The Process of Discounting—Fundamentals

Let us denote by $M=M\left(t\right)$ a given quantity of wealth at time t. In economics, the increment of $M\left(t\right)$ is assumed to be proportional to the quantity itself and the duration of the variation. For a continuous and instantaneous infinitesimal variation, this can be written as

$$dM\left(t\right)\propto M\left(t\right)dt.$$

(4)

This starting phenomenological law is built on the empirical observation that the larger $M\left(t\right)$, the greater its variation at a given time along with the simpler assumption that such a variation is linear in $M\left(t\right)$ and not, for instance, quadratic.

2.1. Definitions and General Setting

We define the interest rate as the relative time derivative

$$r\left(t\right)\equiv \frac{1}{M\left(t\right)}\frac{dM\left(t\right)}{dt}=\frac{dlnM\left(t\right)}{dt},$$

(5)

i.e., the rate is the time derivative of the logarithm of wealth. Let us incidentally note that the linearity shown in Equation (4) is equivalent to assume that $r\left(t\right)$ is independent of $M\left(t\right)$.

In the simplest case, the law (4) represents a direct proportionality, that is to say, $r\left(t\right)$ is constant and and from (5), we see that

$$dM\left(t\right)=rM\left(t\right)dt,$$

(6)

where r is the constant interest rate that has units of 1/(time) (wealth is assumed to be dimensionless). By direct integration, we have the usual exponential law [12]

$$M\left(t\right)=e^{r(t-t_0)}M\left(t_0\right),$$

(7)

which connects wealth at some time $t_0$, for instance, today (which, without loss of generality, we can take equal to zero) with wealth at some future time $t>t_0$.

The growth law (6) appears in many branches of natural and social sciences. Thus, in radioactivity, if $N\left(t\right)$ represents the number of active nuclei at time t, the usual hypothesis is that this number decreases as

$$dN\left(t\right)=-\lambda N\left(t\right)dt.$$

where $\lambda >0$ is the decay constant. Similar considerations apply to many other situations, as in chemical reactions or population dynamics, to name a few.

As we mentioned above, discounting is the procedure of linking wealth at different times. This is done through the discount function defined as

$$\delta \left(t\right)\equiv \frac{M\left(0\right)}{M\left(t\right)},$$

(8)

where $M\left(0\right)$ is today’s wealth. In the case of constant rates, we see from Equation (7) that this function is given by the decreasing exponential

$$\delta \left(t\right)=e^{-rt}.$$

(9)

The assumption of constant rates is actually unrealistic. A first generalization would be to assume that rates are known functions of time $r=r\left(t\right)$. In such a case, the growth law (6) would be given by

$$dM\left(t\right)=r\left(t\right)M\left(t\right)dt,$$

(10)

which, after integrating, yields

$$\delta \left(t\right)=exp\left(-{\int }_0^{t}r\left(t^{\prime }\right)d{t}^{\prime }\right).$$

(11)

However, the assumption of rates being given by constants or by deterministic functions of time is unjustified, in particular over long periods of time. Financial interest rates are typically described as random, as the many models for stochastic interest rates appearing in the literature [13,31,32]. In other words, $r\left(t\right)$ is a random function of time, and in consequence the discount function $\delta \left(t\right)$ given by Equation (11) is a stochastic process. The effective discount function is then defined as the average of $\delta \left(t\right)$,

$$D\left(t\right)=\mathbb{E}\left[exp\left(-{\int }_0^{t}r\left(t^{\prime }\right)d{t}^{\prime }\right)\right],$$

(12)

taken over all possible realizations of $r\left(t\right)$. The function $r\left(t\right)$ can, in principle, be any random process. However, the most natural and simplifying assumption is that rates are Markovian processes with continuous paths—that is, they are diffusion processes [13]. This approach was proposed after the success of taking an identical approach to model stock prices with log-normal process in 1959 and contrasting with empirical data [5]. The first interest rate model was proposed by O. Vasicek in 1977 [33] and during the same decade when stochastic differential equations became crucial to obtain European option prices. Therefore, rates are solutions to stochastic differential equations of the form

$$dr=f\left(r\right)dt+g\left(r\right)dW\left(t\right),$$

(13)

where $W\left(t\right)$ is the Wiener process and the stochastic differential equation is interpreted in the sense of Itô. We assume that drift $f\left(r\right)$ and noise term $g\left(r\right)$ do not depend explicitly on time, that is to say, the time dependence is only implicit through $r=r\left(t\right)$, which means that the interest rate process is time homogeneous and may be stationary [34]. This is certainly an idealization because real markets do not appear to be time-homogeneous, at least over long periods of time [35].

On the other hand, explicit expressions for $f\left(r\right)$ and $g\left(r\right)$ should be proposed based on characteristics obtained by actual data, which is observed to have a reversion toward a mean value, and it is thus claimed to attain a stationary regime in contrast with, for instance, stock market price evolution where no stationary behavior is observed. We return to the comparison between models and empirical data in Section 4.

In order to obtain an operative expression for the effective discount function (12), we define the additional random process

$$x\left(t\right)={\int }_0^{t}r\left(t^{\prime }\right)d{t}^{\prime }.$$

(14)

The interpretation of $x\left(t\right)$ is apparent after substituting Equation (5) into Equation (14) and integrating. We find

$$x\left(t\right)=ln\frac{M\left(t\right)}{M\left(0\right)}\Rightarrow M\left(t\right)=M\left(0\right)e^{x\left(t\right)},$$

which can be taken as an alternative definition of the accumulated return $x\left(t\right)$.

Substituting Equation (14) into Equation (12), we see that the effective discount function can be written as

$$D\left(t\right)=\mathbb{E}\left[e^{-x\left(t\right)}\right],$$

which implies that, in terms of the probability density function (PDF) $p(x,r,t|r_0)$ of the bidimensional diffusion process $\left(x\right(t),r(t\left)\right)$, we can write

$$D\left(t\right|r_0)={\int }_{-\infty }^{\infty }dr{\int }_{-\infty }^{\infty }{e}^{-x}p(x,r,t|r_0)dx,$$

(15)

where we have included the dependence on the initial rate, $r_0=r\left(0\right)$, in the discount function $D\left(t\right|r_0)$.

From Equations (13) and (14), we see that the bidimensional process $\left(x\right(t),r(t\left)\right)$ is defined by the following pair of stochastic differential equations

$$\begin{array}{ccc}\hfill dx& =& rdt,\hfill \\ \hfill dr& =& f\left(r\right)dt+g\left(r\right)dW\left(t\right).\hfill \end{array}$$

(16)

Therefore, the joint density obeys the (forward) Fokker–Planck equation (FPE) [34]

$$\frac{\partial p}{\partial t}=-r\frac{\partial p}{\partial x}-\frac{\partial }{\partial r}\left[f\left(r\right)p\right]+\frac{1}{2}\frac{{\partial }^2}{\partial r^2}\left[g^2\left(r\right)p\right],$$

(17)

with the initial condition

$$p(x,r,0|r_0)=\delta \left(x\right)\delta (r-r_0).$$

(18)

After solving the initial-value problem (17)—(18) and obtaining the joint PDF $p(x,r,t|r_0)$, the discount function follows from Equation (15). There are, however, two different approaches for achieving it. One of them, which is standard in the financial literature, is based on the backward Fokker–Planck equation, and this is called the Feynman–Kac approach [13]. A second procedure is based on Fourier analysis [27]. We will explain both approaches next.

2.2. The Feynman–Kac Approach

Using this method, one obtains a partial differential for the discount function $D\left(t\right|r_0)$, which is based on the backward Fokker–Planck equation for the joint density $p(x,r,t|r_0)$. In what follows, we assume that $t_0\ne 0$ and denote $x_0=x\left(t_0\right)$. By definition, $x_0=0$ (cf. Equation (14)). However, we temporally keep $x_0\ne 0$ and set $x_0=0$ at the end of the calculation when needed.

The backward FPE for the PDF $p(x,r,t|x_0,r_0,t_0)$ that corresponds to the bidimensional process (16) is [34]

$$\frac{\partial p}{\partial t_0}=-r_0\frac{\partial p}{\partial x_0}-f\left(r_0\right)\frac{\partial p}{\partial r_0}-\frac{1}{2}{g}^2\left(r_0\right)\frac{{\partial }^{2}p}{\partial r_0^2},$$

(19)

with the final condition as $t_0\to t$,

$$p(x,r,t|x_0,r_0,t)=\delta (x-x_0)\delta (r-r_0).$$

(20)

Let us observe that the problem (19)—(20) is invariant under translations of both time and $x_0$. We thus define the new variables

$$t^{\prime }=t-t_0,x^{\prime }=x-x_0,$$

(21)

so that

$$\frac{\partial p}{\partial t_0}=-\frac{\partial p}{\partial t^{\prime }},\frac{\partial p}{\partial x_0}=-\frac{\partial p}{\partial x^{\prime }},$$

and Equation (19) reads

$$\frac{\partial p}{\partial t^{\prime }}=-r_0\frac{\partial p}{\partial x^{\prime }}+f\left(r_0\right)\frac{\partial p}{\partial r_0}+\frac{1}{2}{g}^2\left(r_0\right)\frac{{\partial }^{2}p}{\partial r_0^2}.$$

(22)

Under this change of variables, we also have

$$p=p(x,r,t|x_0,r_0,t_0)=p(x,r,t|x-x^{\prime },r_0,t-t^{\prime })=p(x^{\prime },r,t^{\prime }|r_0),$$

where the last equality comes from the invariance under time and x translations, that is,

$$p(x,r,t|x_0,r_0,t_0)=p(x-x_0,r,t-t_0|r_0).$$

Consequently, the final condition (20) becomes the initial condition

$$p(x^{\prime },r,t^{\prime }=0|r_0)=\delta \left(x^{\prime }\right)\delta (r-r_0).$$

(23)

Having set the backward FPE in the form given by Equation (22), we next obtain the equation satisfied by the effective discount $D\left(t\right|r_0)$. To this end, we multiply Equation (22) by $e^{-x^{\prime }}$ and integrate over $x^{\prime }$ and r, we have

$$\begin{array}{ccc}\hfill \frac{\partial }{\partial t^{\prime }}{\int }_{-\infty }^{\infty }dr{\int }_{-\infty }^{\infty }{e}^{-x^{\prime }}pd{x}^{\prime }& =& -r_0{\int }_{-\infty }^{\infty }dr{\int }_{-\infty }^{\infty }{e}^{-x^{\prime }}\frac{\partial p}{\partial x^{\prime }}d{x}^{\prime }\hfill \\ & & +\left[f\left(r_0\right)\frac{\partial }{\partial r_0}+\frac{1}{2}{g}^2\left(r_0\right)\frac{{\partial }^2}{\partial r_0^2}\right]{\int }_{-\infty }^{\infty }dr{\int }_{-\infty }^{\infty }{e}^{-x^{\prime }}pd{x}^{\prime }.\hfill \end{array}$$

(24)

From Equation (15), we see that

$${\int }_{-\infty }^{\infty }dr{\int }_{-\infty }^{\infty }{e}^{-x^{\prime }}p(x^{\prime },r,t^{\prime }|r_0)d{x}^{\prime }=D\left(t^{\prime }\right|r_0).$$

(25)

On the other hand, integrating by parts, the first integral on the right hand side of Equation (24) and using (25), we have

$${\int }_{-\infty }^{\infty }dr{\int }_{-\infty }^{\infty }{e}^{-x^{\prime }}\frac{\partial p}{\partial x^{\prime }}d{x}^{\prime }={\int }_{-\infty }^{\infty }dr{\int }_{-\infty }^{\infty }{e}^{-x^{\prime }}p(x^{\prime },r,t^{\prime }|r_0)d{x}^{\prime }=D\left(t^{\prime }\right|r_0),$$

(26)

where we have considered the boundary condition (otherwise implicit in the definition of D given in Equation (15))

$$\underset{x^{\prime }\to \pm \infty }{lim}\left[e^{-x^{\prime }}p(x^{\prime },r,t^{\prime }|r_0)\right]=0.$$

Substituting Equations (25) and (26) into Equation (24) and setting $t_0=0$, which implies $t^{\prime }=t$ [cf. Equation (21)], we finally obtain

$$\frac{\partial D}{\partial t}=-r_{0}D+f\left(r_0\right)\frac{\partial D}{\partial r_0}+\frac{1}{2}{g}^2\left(r_0\right)\frac{{\partial }^{2}D}{\partial r_0^2},$$

(27)

with the initial condition (cf. Equations (23) and (25))

$$D\left(0\right|r_0)=1.$$

(28)

The method for obtaining the discount function $D\left(t\right|r_0)$ by solving the initial-value problem (27)—(28) is called the Feynman–Kac approach, and Equation is (27) the Feynman–Kac equation. In some applications (see, for instance, Section 3), it is convenient to consider $t_0\ne 0$ so that $t^{\prime }=t-t_0\ne t$. In these cases, it is appropriate to denote $D=D\left(t\right|r_0,t_0)$ and the Feynman–Kac Equation (27) reads

$$\frac{\partial D}{\partial t_0}=r_{0}D-f\left(r_0\right)\frac{\partial D}{\partial r_0}-\frac{1}{2}{g}^2\left(r_0\right)\frac{{\partial }^{2}D}{\partial r_0^2},$$

(29)

with the final condition $D\left(t\right|r_0,t)=1$.

2.3. The Fourier Transform Approach

An alternative method for obtaining the discount function is based on the joint characteristic function—that is, on the Fourier transform of the joint density,

$$\tilde{p}({\omega }_1,{\omega }_2,t|r_0)={\int }_{-\infty }^{\infty }{e}^{-i{\omega }_{2}r}dr{\int }_{-\infty }^{\infty }{e}^{-i{\omega }_{1}x}p(x,r,t|r_0)dx.$$

(30)

One of the chief advantages of working with the characteristic function is that obtaining the effective discount is straightforward. Indeed, comparison of Equation (15),

$$D\left(t\right|r_0)={\int }_{-\infty }^{\infty }dr{\int }_{-\infty }^{\infty }{e}^{-x}p(x,r,t|r_0)dx,$$

with Equation (30) shows that

$$D\left(t\right|r_0)=\tilde{p}\left({\omega }_1=-i,{\omega }_2=0,t|r_0\right).$$

(31)

Therefore, in order to obtain the discount function, we only need to know the joint characteristic function of the bidimensional process $(x,r)$. The procedure is quite advantageous in linear cases. In a forthcoming section, we will apply this approach to some standard models of interest rates.

2.4. Adding Risk Aversion

As we will see in the next section, the process of discounting just described is very closely related to an important problem in finance called bond pricing. In the context of bond pricing, there can be two kinds of investors. For one hand, if investors are risk neutral, then bond prices can be modeled based on the data generating measure p, which is the solution of the Fokker–Planck Equation (17) with the initial condition (18). This is sometimes called the Local Expectation Hypothesis (LEH) [36,37]. Nonetheless, a more general assumption is that investors are sensitive to risk.

In such a case, bonds are somewhat more accurately priced using an artificial density $p^{*}$ usually called a risk-neutral (or risk-correcting) probability measure. Both magnitudes, the data generating measure p and the risk-neutral measure $p^{*}$, are related through a quantity that is denoted by $q(r,t)$ and called market price of risk, which, as described in the next section, is the extra return per unit of risk that investors demand to bear risk. This additional return is thus determined by a function $q=q(r,t)$ that, in its most general form, may depend on the rate r and current time t, although the most usual assumption is that $q=q\left(r\right)$ only depends on the rate [33]. Following a standard procedure for bond pricing [33,38], which we will present in Section 3, one takes risk into account by replacing the drift $f\left(r\right)$ by $f^{*}\left(r\right)$,

$$f\left(r\right)\to f^{*}\left(r\right),$$

where

$$f^{*}\left(r\right)=f\left(r\right)+g\left(r\right)q\left(r\right),$$

(32)

and $q\left(r\right)\ge 0$ is the market price of risk. The form of $q\left(r\right)$ is, in principle, unknown and has to be conjectured. The simplest and most common assumption is that $q\left(r\right)=q$ is constant, in such a case, the value of q may be more easily estimated from empirical data. Now, the risk-neutral measure $p^{*}(x,r,t|r_0)$ is given by the Fokker–Planck Equation (17) with $f\left(r\right)$ replaced by $f^{*}\left(r\right)$; that is,

$$\frac{\partial p^{*}}{\partial t}=-r\frac{\partial p^{*}}{\partial x}-\frac{\partial }{\partial r}\left[\left[f\left(r\right)+g\left(r\right)q\left(r\right)\right]p^{*}\right]+\frac{1}{2}\frac{{\partial }^2}{\partial r^2}\left[g^2\left(r\right)p^{*}\right],$$

(33)

with the initial condition given by

$$p^{*}(x,r,0|r_0)=\delta \left(x\right)\delta (r-r_0).$$

(34)

In an analogous way, the discount function adjusted for risk will now be given by the Feynman–Kac Equation (27) with $f\left(r\right)$ replaced by $f^{*}\left(r\right)$. Or, using the Fourier method, the discount function will be given in terms of the risk-neutral characteristic function, ${\tilde{p}}^{*}({\omega }_1,{\omega }_2,t|r_0)$, by (cf. Equation (31))

$$D\left(t\right|r_0)={\tilde{p}}^{*}\left({\omega }_1=-i,{\omega }_2=0,t|r_0\right).$$

(35)

3. Pricing Bonds—The Term Structure of Interest Rates

Pricing bonds is a traditional objective in finance and intimately related to the problem of discounting. It constitutes a vast subject with countless studies, many of them rather abstract, which have appeared in the mathematical finance literature over the last decades. We present a short and intentionally simple, yet rigorous, introduction to the subject devoid as much as possible of technicalities and mathematical subtleties and refer the interested reader to more specialized works for further information [13].

A bond is a financial instrument that one purchases now and that provides a payment in the future. From a more technical point of view, we say that a (discount) bond is a default-free claim on a specified sum of money to be delivered at a given future date called the maturity time. Such claims are bought and issued by investors. Let us denote by $B(t_0,t)$ the price at time $t_0$ of a discount bond maturing at time $t\ge t_0$, with unit maturity value,

$$B(t,t)=1.$$

Let us incidentally note that, if the final maturity price is not 1 (say, $B(t,t)=\beta$) then the price of the bond at $t_0$ would be $\beta B(t_0,t)$.

Bonds are classified according to the time interval to maturity $\tau$ defined as

$$\tau =t-t_0.$$

Thus, if $\tau =10$ years, we talk about a 10-year bond that is traded at $t_0$ (for instance, today) with price $B(t_0,t_0+10)$ and that, after 10 years, has unit value. Similarly for a 3-year bond, 3-month bond, etc.

The central question is to know the backward evolution of the bond price, from unit maturity to the initial purchasing price $B(t_0,t)$. Note that the problem is virtually identical to the problem of discounting discussed in Sect. II, with the sole difference that, in discounting, we look for the forward evolution from a known initial value to an unknown final value, while, in bond pricing, the situation is reversed, since we know the final value but not the initial one.

In order to proceed further, we define the instantaneous rate of return $b(t_0,t)$ (also called forward rate) as the relative time variation of the bond price (compare with Equation (5))

$$b(t_0,t)\equiv \frac{1}{B(t_0,t)}\frac{\partial B(t_0,t)}{\partial t_0}=\frac{\partial lnB(t_0,t)}{\partial t_0}.$$

(36)

The knowledge of the forward rate $b(t_0,t)$ allows us to relate the initial price $B(t_0,t)$ and the maturing price $B(t,t)=1$. Indeed, the integration of the above equation directly leads to

$$B(t_0,t)=exp\left[-{\int }_{t_0}^{t}b(t_0^{\prime },t)d{t}_0^{\prime }\right].$$

(37)

The close analogy between bond pricing and discounting is now apparent. Indeed, the comparison of Equation (37) with Equation (11) shows that $B(t_0,t)$ is the equivalent of the discount function $\delta \left(t\right)$ and that the forward rate $b(t_0,t)$ is the equivalent of the discount rate $r\left(t\right)$. However, in what follows, we will use the notation $r\left(t\right)$ not for the forward rate $b(t_0,t)$ but for the so-called spot rate (also called nominal rate), which we define in Equation (39).

Another quantity of interest is the yield to maturity $y(t_0,\tau )$ defined by

$$y(t_0,\tau )\equiv -\frac{1}{\tau }lnB(t_0,t_0+\tau )\Rightarrow B(t_0,t_0+\tau )=e^{-\tau y(t_0,\tau )}.$$

(38)

From (37), we see that

$$y(t_0,\tau )=\frac{1}{\tau }{\int }_{t_0}^{t_0+\tau }b(t_0^{\prime },t)d{t}_0^{\prime },$$

that is to say, the yield is the time average of the forward rate over the maturity period $\tau$.

A final quantity is needed, the spot or nominal rate, which is defined as the limit of the yield when the maturity tends to 0. In dealing with bonds, one sometimes uses, for the nominal rate, the notation $n\left(t_0\right)$ instead of $r\left(t_0\right)$—the later reserved for real interest rates, which can be negative due to inflation (see Section 3).

$$r\left(t_0\right)\equiv \underset{\tau \to 0}{lim}y(t_0,\tau )=\underset{\tau \to 0}{lim}\left[\frac{1}{\tau }{\int }_{t_0}^{t_0+\tau }b(t_0^{\prime },t)d{t}_0^{\prime }\right].$$

(39)

Solving the indeterminacy by expanding the integral in powers of $\tau$, we see that the spot rate is given in terms of the forward rate by

$$r\left(t_0\right)=b(t_0,t_0).$$

(40)

In other words, the spot rate is the instantaneous forward rate.

Let us finally note that a loan of amount M subscribed at time $t_0$ with an interest rate $r\left(t_0\right)$ (the spot rate) will, at time $t_0+d{t}_0$, increase in value to $M+dM$, where

$$dM=r\left(t_0\right)Md{t}_0.$$

(41)

Indeed, at any time $t_0$, the value of the spot rate $r\left(t_0\right)$ is the instantaneous increase of the loan value, that is, $r\left(t_0\right)=dlnM\left(t_0\right)/d{t}_0$ (compare with Equation (36)). All of this clearly heightens the close similarities with discounting mentioned above.

However, subsequent values of the spot rate are not necessarily certain. We will see next, the consequences of this fact on the time evolution of the bond price $B(t_0,t)$.

3.1. Dynamics of the Bond Price

Suppose the spot rate $r\left(t_0\right)$ is not deterministic but random. In such a case, and analogously to discounting, the usual assumption is that $r_0=r\left(t_0\right)$ is a Markovian random process with continuous trajectories; that is, a diffusion process obeying a stochastic differential equation of the form

$$d{r}_0=f\left(r_0\right)d{t}_0+g\left(r_0\right)dW\left(t_0\right),$$

(42)

where $W\left(t_0\right)$ is the standard Wiener process. We assumed that the drift and the noise intensity are independent of time, thus, the time dependence of these coefficients is implicit through $r_0=r\left(t_0\right)$. We know that this implies invariance under time translations, and we can set $t_0=0$ when needed without loss of generality.

We will now obtain the time evolution of the bond price $B(t_0,t)$ from the purchasing time $t_0$ to maturity t, and to this end, we follow Oldrich Vasicek [33]. Let us first observe that the most natural hypothesis consists in assuming that the bond price B is a function of the initial spot rate $r_0=r\left(t_0\right)$ and write

$$B=B\left[t_0,t|r\left(t_0\right)\right].$$

(43)

In this way, $B(t_0,t|r_0)$ represents the price of a bond issued at time $t_0$ and maturing at time t, given that the initial interest rate is $r_0=r\left(t_0\right)$. The infinitesimal variation of the bond price is then defined by

$$dB=B[t_0+d{t}_0,t|r(t_0+d{t}_0)]-B\left[t_0,t|r\left(t_0\right)\right].$$

We expand in Taylor series up to second order

$$\begin{array}{ccc}\hfill B[t_0+d{t}_0,t\left|r\right(t_0& +& d{t}_0\left)\right]=B\left[t_0,t|r\left(t_0\right)\right]+\frac{\partial B}{\partial t_0}d{t}_0+\frac{\partial B}{\partial r_0}d{r}_0\hfill \\ & +& \frac{1}{2}\left[\frac{{\partial }^{2}B}{\partial t_0^2}d{t}_0^2+\frac{{\partial }^{2}B}{\partial r_0^2}d{r}_0^2+2\frac{{\partial }^{2}B}{\partial t_0\partial r_0}d{t}_{0}d{r}_0\right]+\cdots .\hfill \end{array}$$

Substituting for Equation (42) and taking into account that $dW\left(t_0\right)=O\left(d{t}_0^{1/2}\right)$ [34,39], we write

$$dB=\left[\frac{\partial B}{\partial t_0}+f\left(r_0\right)\frac{\partial B}{\partial r_0}\right]d{t}_0+g\left(r_0\right)\frac{\partial B}{\partial r_0}dW\left(t_0\right)+\frac{1}{2}{g}^2\left(r_0\right)\frac{{\partial }^{2}B}{\partial r_0^2}{\left[dW\left(t_0\right)\right]}^2+O\left(d{t}_0^{3/2}\right).$$

However, ${\left[dW\left(t_0\right)\right]}^2=d{t}_0$ (in mean square sense) [34,39], and, up to the first order in $d{t}_0$, we obtain

$$dB=\left[\frac{\partial B}{\partial t_0}+f\left(r_0\right)\frac{\partial B}{\partial r_0}+\frac{1}{2}{g}^2\left(r_0\right)\frac{{\partial }^{2}B}{\partial r_0^2}\right]d{t}_0+g\left(r_0\right)\frac{\partial B}{\partial r_0}dW\left(t_0\right).$$

(44)

Defining

$$\mu (t_0,t|r_0)\equiv \frac{1}{B}\left[\frac{\partial B}{\partial t_0}+f\left(r_0\right)\frac{\partial B}{\partial r_0}+\frac{1}{2}{g}^2\left(r_0\right)\frac{{\partial }^{2}B}{\partial r_0^2}\right],$$

(45)

and

$$\sigma (t_0,t|r_0)\equiv -\frac{1}{B}g\left(r_0\right)\frac{\partial B}{\partial r_0},$$

(46)

we see from (44) that the bond price satisfies the stochastic differential equation

$$\frac{dB}{B}=\mu (t_0,t|r_0)d{t}_0-\sigma (t_0,t|r_0)dW\left(t_0\right),$$

(47)

showing that the bond price is also a diffusion process.

Averaging Equation (47) and recalling that $\mathbb{E}\left[dW\left(t_0\right)\right]=0$, we see that

$$\mu (t_0,t|r_0)=\mathbb{E}\left[\frac{1}{B}\frac{dB}{d{t}_0}\right],$$

which proves that $\mu (t_0,t|r_0)$ is the average of the instantaneous rate of return [cf. Equation (36)] at time $t_0$ on a bond with maturing date t, given that the current spot rate is $r_0$. In an analogous way, one can easily show that ${\sigma }^2(t_0,t|r_0)$ is the variance [33].

We therefore see from the above development that the bond price is a random quantity but with a final fixed price (the maturity price). The question is: what is the (initial) price that an investor has to buy (or sell) a bond at time $t_0$ maturing at time $t=t_0+\tau$ with the current spot rate $r_0$? One possible answer would be proceeding as in discounting to take the average over all possible realizations of the bond price. However, this procedure implies that the expected rate of return of a bond is invariant under risk variation—that is, under changes of the variance ${\sigma }^2(t_0,t|r_0)$—a fact that investors always have in mind.

We explain next a procedure resulting in a deterministic bond price, which takes into account the risk aversion of investors (in practice this is only true to some extend because the mathematical procedure assumes that the market is driven by Gaussian white noise—that is, the Wiener process, which is an idealized noise presenting, among other shortcomings, no fat tails, a key characteristic of real markets [6]).

3.2. The Market Price of Risk

Consider an investor who, at time $t_0$, sells an amount $M_1$ of a bond maturing at time $t_1$ and, at the same time, buys an amount $M_2$ of another bond with a different maturing date $t_2$. The total worth of the portfolio, thus, constructed is $M=M_2-M_1$. Note that each amount $M_i$ ($i=1,2$) is a multiple of the bond price $B(t_0,t_i|r_0)$ ($i=1,2$) and, hence, they also obey the stochastic differential Equation (47). That is,

$$\frac{d{M}_i}{M_i}=\mu (t_0,t_i|r_0)d{t}_0-\sigma (t_0,t_i|r_0)dW\left(t_0\right).$$

As a consequence, the infinitesimal variation $dM=d{M}_2-d{M}_1$ of the worth of the portfolio changes over time according to

$$\begin{array}{ccc}\hfill dM& =& \left[\mu (t_0,t_2|r_0)M_2-\mu (t_0,t_1|r_0)M_1\right]d{t}_0\hfill \\ & & -\left[\sigma (t_0,t_2|r_0)M_2-\sigma (t_0,t_1|r_0)M_1\right]dW\left(t_0\right).\hfill \end{array}$$

(48)

Suppose we choose the amounts $M_1$ and $M_2$ such that

$$M_1=\frac{M}{{\sigma }_1-{\sigma }_2}{\sigma }_2,M_2=\frac{M}{{\sigma }_1-{\sigma }_2}{\sigma }_1,$$

(49)

where $M=M_2-M_1$ and ${\sigma }_i=\sigma (t_0,t_i|r_0)$ ($i=1,2$). Hence $M_1$ is proportional to ${\sigma }_2$, while $M_2$ is proportional to ${\sigma }_1$. With this choice, we have

$${\sigma }_2{M}_2-{\sigma }_1{M}_1={\sigma }_2\frac{{\sigma }_{1}M}{{\sigma }_1-{\sigma }_2}-{\sigma }_1\frac{{\sigma }_{2}M}{{\sigma }_1-{\sigma }_2}=0,$$

and the random term in Equation (48) vanishes. This renders the portfolio composed of such amounts of the two bonds instantaneously riskless:

$$dM=\frac{M}{{\sigma }_1-{\sigma }_2}({\mu }_2{\sigma }_1-{\mu }_1{\sigma }_2)d{t}_0,$$

(50)

where ${\mu }_i=\mu (t_0,t_i|r_0)$. The rate of return $r_M$ of this portfolio is

$$r_M\equiv \frac{1}{M}\frac{dM}{d{t}_0}=\frac{{\mu }_2{\sigma }_1-{\mu }_1{\sigma }_2}{{\sigma }_1-{\sigma }_2}.$$

In order to avoid arbitrage opportunities—that is, making profits without taking any risk—the rate $r_M$ must be equal to the spot rate $r_0$. If not, the portfolio can be purchased by taking funds borrowed at the spot rate, or otherwise sold and the profits lent out to accomplish a riskless arbitrage [33]. Therefore (compare also Equation (41) with Equation (50))

$$r_0=\frac{{\mu }_2{\sigma }_1-{\mu }_1{\sigma }_2}{{\sigma }_1-{\sigma }_2}.$$

Rearranging terms, we find $({\mu }_1-r_0)/{\sigma }_1=({\mu }_2-r_0)/{\sigma }_2$, so that

$$\frac{\mu (t_0,t_1|r_0)-r_0}{\sigma (t_0,t_1|r_0)}=\frac{\mu (t_0,t_2|r_0)-r_0}{\sigma (t_0,t_2|r_0)}.$$

This equation is valid for arbitrary maturities $t_1,t_2,\dots$, it then follows that the ratio $[\mu (t_0,t|r_0)-r_0]/\sigma (t_0,t|r_0)$must be independent of the maturity timet.

Let us denote by $q\left(t_0\right|r_0)$ the common value of such a ratio for a bond of any maturity date, given that the current spot rate (at time $t_0$) is $r_0$,

$$q\left(t_0\right|r_0)\equiv \frac{\mu (t_0,t|r_0)-r_0}{\sigma (t_0,t|r_0)},(t\ge t_0).$$

(51)

The quantity $q\left(t_0\right|r_0)$ is called the market price of risk, as it gives the variation of the expected rate of return on a bond (specified by the risk premium$\mu -r_0$) per an additional unit risk (specified by the standard deviation $\sigma$). The market price of risk $q\left(t_0\right|r_0)$ is the so-called Sharpe ratio [40] of the excess return $\mu -r_0$.

Note that, if $q=0$, the spot rate $r_0=r\left(t_0\right)$ and the average rate of return $\mu$ coincide.

$$\mu (t_0,t|r_0)=r\left(t_0\right)$$

$(t=t_0+\tau )$ meaning that the expected instantaneous rates of return on bonds are the same for all maturities.

3.3. The Term Structure Equation and the Risk-Neutral Measure

The introduction of the market price of risk implies a non-random bond price $B=B(t_0,t|r_0)$, which, in turn, allows a deterministic equation for B. In effect, rewriting Equation (51) as

$$\mu (t_0,t|r_0)-r_0=\sigma (t_0,t|r_0)q\left(t_0\right|r_0),$$

and substituting $\mu$ and $\sigma$ for their definitions given in Equations (45) and (46), we have

$$\frac{1}{B}\left[\frac{\partial B}{\partial t_0}+f\left(r_0\right)\frac{\partial B}{\partial r_0}+\frac{1}{2}{g}^2\left(r_0\right)\frac{{\partial }^{2}B}{\partial r_0^2}\right]-r_0=-q\left(t_0\right|r_0)\frac{1}{B}g\left(r_0\right)\frac{\partial B}{\partial r_0},$$

which, after rearranging terms, yields

$$\frac{\partial B}{\partial t_0}=r_{0}B-[f\left(r_0\right)+g\left(r_0\right)q\left(t_0\right|r_0)]\frac{\partial B}{\partial r_0}-\frac{1}{2}{g}^2\left(r_0\right)\frac{\partial B}{\partial r_0^2}.$$

(52)

This equation, called the term structure equation, is a partial differential equation for $B(t_0,t|r_0)$, that is obtained once we know the random character of the spot rate process $r\left(t\right)$ (through its drift f and noise intensity g) and once the market price of risk $q\left(t_0\right|r_0)$ is specified. Bond prices are thus obtained after solving the deterministic Equation (52) with the final condition:

$$B(t,t|r_0)=1.$$

(53)

Let us observe that the term structure Equation (52) for the bond price B is identical to the Feynman–Kac Equation (27) for the discount function D as long as we make the following change of drift

$$f\left(r_0\right)⟶f\left(r_0\right)+g\left(r_0\right)q\left(t_0\right|r_0).$$

(54)

On the other hand, as we have seen in Section 2, the solution of the Feynman–Kac Equation (27) for the discount function $D\left(t\right|r_0)$ is written as the average (cf. Equation (15))

$$D\left(t\right|r_0,t_0)={\int }_{-\infty }^{\infty }dr{\int }_{-\infty }^{\infty }{e}^{-x}p(x,r,t|r_0,t_0)dx,$$

where $p(x,r,t|r_0,t_0)$ is the probability density function of the bidimensional diffusion process $\left(x\right(t),r(t\left)\right)$ defined by Equation (16),

$$dx=rdt,dr=f\left(r\right)dt+g\left(r\right)dW\left(t\right).$$

Now the analogy between the term structure Equation (52) and the Feynman–Kac Equation (29) suggests that we can write the bond price $B(t_0,t|r_0)$ as an average over the different realizations of the spot rate $r\left(t_0\right)$. However, this averaging procedure is taken using a modified PDF called the risk-free measure. Thus, it can be proven in a more rigorous way that [32,33]

$$B(t_0,t|r_0)={\int }_{-\infty }^{\infty }dr{\int }_{-\infty }^{\infty }{e}^{-x}{p}^{*}(x,r,t|r_0,t_0)dx,$$

(55)

where $p^{*}(x,r,t|r_0,t_0)$ is the risk-free measure that is the PDF of the bidimensional process $(x\left(t_0\right),r\left(t_0\right))$ defined by the following pair of stochastic differential equations that include the market price of risk (see Equation (54)):

$$\begin{array}{ccc}\hfill dx& =& rdt,\hfill \\ \hfill dr& =& \left[f\right(r)+g(r\left)q\right(t\left|r\right)]dt+g(r\left)dW\right(t).\hfill \end{array}$$

(56)

That is, $p^{*}$ is the solution to the FPE

$$\frac{\partial p^{*}}{\partial t}=-r\frac{\partial p^{*}}{\partial x}-\frac{\partial }{\partial r}\left[\left[f\left(r\right)+g\left(r\right)q\left(t\right|r)\right]p^{*}\right]+\frac{1}{2}\frac{{\partial }^2}{\partial r^2}\left[g^2\left(r\right)p^{*}\right],$$

(57)

with the initial condition

$$p^{*}(x,r,t_0|r_0,t_0)=\delta \left(x\right)\delta (r-r_0).$$

(58)

Since, as we have shown in Section 2.2, the Feynman–Kac approach to discounting is equivalent to the Fourier method described in Section 2.3, we can apply the latter to directly obtain the bond price knowing only the risk neutral PDF, without having to solve the Feynman–Kac Equation (52) with condition (23). Indeed, the characteristic function of the risk neutral density $p^{*}$ is the joint Fourier transform

$${\tilde{p}}^{*}({\omega }_1,{\omega }_2,t|r_0,t_0)={\int }_{-\infty }^{\infty }{e}^{-i{\omega }_{2}r}dr{\int }_{-\infty }^{\infty }{e}^{-i{\omega }_{1}x}{p}^{*}(x,r,t|r_0,t_0)dx,$$

which, after comparing with Equation (55), yields

$$B(t_0,t|r_0)={\tilde{p}}^{*}\left({\omega }_1=-i,{\omega }_2=0,t|r_0,t_0\right).$$

(59)

Finally, once we know the bond price, the yield to maturity $y(t_0,\tau |r_0)$ (also called the term structure of interest rates) is readily evaluated from Equation (38):

$$y(t_0,\tau |r_0)=-\frac{1}{\tau }lnB(t_0,t_0+\tau |r_0).$$

(60)

The graphic representations of $y(t_0,\tau |r_0)$ as a function of $t_0$ and for different values of the maturity interval $\tau$ are called yield curves and are of prime importance for practitioners.

4. Standard Models

Throughout the above development, it is clear that, in order to proceed further in the discounting process (as well as in pricing bonds), we need to identify the specific diffusion process chosen for modeling rates. Such a choice is mostly based on the analysis of empirical data [26,29]. Clearly, in any proposed market model, there are idealizations, otherwise a complete treatment of the problem would be very problematic, not to say impossible, not only analytically but computationally as well.

In addition to assuming a diffusive behavior for the market, the first of such idealizations is supposing that the market is stationary, i.e., the structural conditions of the market do not change over time. However, and in particular after the 1980s, market circumstances have largely changed due to a great increase of transaction volumes along with transparency and, to a lesser extend, changes in investor perspectives. In this review, we only address stationary models, although we are working on new models dealing with some non-stationary features of the market [29,35], and we refer the interested reader to these works for further information.

On the other hand, there is a property of the market that appears to be well founded on empirical grounds. This is the property of mean reversion meaning that prices tend to return to some fundamental value, called the normal level, which is typically identified as the long-time (i.e., stationary) mean value. The simplest method of introducing this feature in the diffusion market model is to assume a linear drift of the form $f\left(r\right)=-\alpha (r-m)$, where $\alpha >0$ is the strength of the reversion to the normal level, identified by m.

In such a case, the drift acts like a linear restoring force driving $r\left(t\right)$ towards m as time increases. Despite that the introduction of mean reversion might create some arbitrage opportunities, the property of mean reversion is widely accepted in the literature [6], and we previously discussed this issue in the context of option pricing when considering the Ornstein–Uhlenbeck model [41].

4.1. Bonds and Real Rates

Before proceeding with the introduction of some standard models for the market evolution, we briefly explain the link between bonds and (real) interest rates.

Financial economists have developed a large number of models of interest rate processes to enable them to price bonds and other cash flows. In these models, interest rates are described by positive random processes since financial interest rates rarely take negative values. Although the models could be, in principle, extended to arbitrary horizons, they have only been studied carefully over time horizons of up to 30 years, since bonds are seldom issued for periods longer than this.

On the other hand, environmental economists are interested in the real behavior of the economic growth over longer horizons, in contrast to financial economists who are typically more interested in nominal rates over shorter periods of time. The behavior of real and nominal rates usually differ as, due to inflation, real rates can take on negative values. In this way, real rates $r\left(t\right)$ are generally defined by the so-called Fisher procedure:

$$r\left(t\right)=n\left(t\right)-i\left(t\right),$$

(61)

where $i\left(t\right)$ is the inflation rate that is usually generated from consumer price indexes as we will explain in the next section. The quantity $n\left(t\right)$ represents nominal rates, which are typically constructed out of government bonds and are usually positive (even though, in recent years, nominal rates have taken slightly negative values). In order to explain the close relationship between nominal rates and bonds, let us first recall that nominal rates were called spot rates in the previous section on bond pricing where we used the notation $r\left(t\right)$ instead of $n\left(t\right)$. We thus define the spot (i.e., nominal) rate as (see Equations (39) and (40))

$$n\left(t\right)\equiv \underset{\tau \to 0}{lim}\left[\frac{1}{\tau }{\int }_t^{t+\tau }b(t^{\prime },t)d{t}^{\prime }\right]=b(t,t),$$

(62)

where $b(t^{\prime },t)$ is the forward rate for bonds defined in Equation (36) (see also Equation (37)), that is

$$b(t^{\prime },t)=\frac{\partial lnB(t^{\prime },t)}{\partial t^{\prime }}\Rightarrow B(t^{\prime },t)=exp\left[-{\int }_{t^{\prime }}^{t}b(t^{″},t)d{t}^{″}\right]$$

where $B(t^{\prime },t)$ is the price at time $t^{\prime }$ of a (government) bond maturing at time $t\ge t^{\prime }$. Let us recall the definition of the yield to maturity $y(t,\tau )$ given in Equation (38),

$$y(t,\tau )\equiv -\frac{1}{\tau }lnB(t,t+\tau )\Rightarrow B(t,t+\tau )=e^{-\tau y(t,\tau )}$$

($\tau \ge 0$), so that

$$y(t,\tau )=\frac{1}{\tau }{\int }_t^{t+\tau }b(t^{\prime },t)d{t}^{\prime }$$

and comparing with Equation (62), we see that, in terms of the yield nominal rates, $n\left(t\right)$ can be defined as

$$n\left(t\right)=\underset{\tau \to 0}{lim}y(t,\tau ).$$

(63)

Thus, for empirical analysis, the yield can be used as an estimator of the nominal rates:

$$n\left(t\right)\sim y(t,\tau ),$$

(64)

and the accuracy of such an estimator increases as $\tau \to 0$. We will return to this discussion in the next section.

In this way, taking nominal rates corrected by inflation as a proxy of economic growth, we recently demonstrated [26,29,30] through a detailed empirical study on many countries that real interest rates are negative around 25% of the time (see next section). To understand how discounting depends on the random process used to characterize interest rates, we focused on three different models and obtained exact analytical expressions for the discount function [27]. The three models describe to, varying degrees, a number of relevant characteristics observed in the data, while being simple enough to allow for complete analytical treatment. The main results are summarized in

Table 1. Key statistical features for three standard models: the Vasicek (Ornstein–Uhlenbeck), the Cox–Ingersoll–Ross (Feller) and the log-normal models. The average and variance are provided in terms of the model parameters to better compare the asymptotic behavior of $D\left(t\right)$. The asymptotic discount is provided by showing an exponential decay with a long-run rate of discount $r_{\infty }$ for the Vasicek and the Cox–Ingersoll–Ross models and also in the log-normal case for a specific combination of parameters ($k^2/2<\alpha$, mild fluctuations). The parameter $\delta$ is defined in Equation (111).

Model	$\mathbb{E}\left[\mathit{r}\right(\mathit{t}\left)\right]$	Var$\left[\mathit{r}\right(\mathit{t}\left)\right]$	$\mathit{D}(\mathit{t}\to \mathbf{\infty })$		${\mathit{r}}_{\mathbf{\infty }}$
Vasicek	m	$k^2/\alpha$	$exp(-r_{\infty }t)$		$m-k^2/2{\alpha }^2$
Feller	m	$m{k}^2/\left(2\alpha \right)$	$exp(-r_{\infty }t)$		$\frac{2m}{1+\sqrt{1+2k^2/\alpha }}$
Log-normal	$r_0{e}^{\alpha t}$	$r_0^2{e}^{2\alpha t}[e^{k^{2}t}-1]$	constant	($k^2/2>\alpha$)	$--$
			$exp(-r_{\infty }t)$	($k^2/2<\alpha$)	$(\alpha -k^2/2)/\delta$
			$t^{-1/2}$	($k^2/2=\alpha$)	$--$

.

The first model is based on the Ornstein–Uhlenbeck (OU) process—also called the Vasicek model in the financial literature [13]—which allows for negative rates and is, therefore, suitable for pricing environmental problems. The model has a stationary probability distribution and exhibits reversion to the mean, which means that the process tends to return to its average stationary value. We will review this model below.

The second and third models that we considered are given by the Feller and log-normal processes, respectively. For these processes, the rates cannot be negative. The Feller process—also known as the Cox–Ingersoll–Ross (CIR) model [42]—has reversion to the mean and stationary probability distribution.

This is one of the most popular models in finance [13], and we recently reviewed the main properties of the Feller process in previous works [27,43]. A third model, also implying positive rates, is the log-normal process (occasionally called the Dotham model in the financial literature [44]). The model does not have reversion to the mean nor a stationary distribution. Despite these shortcomings, the log-normal process has also been used in the financial literature mainly because it is positive and allows for analytical treatment [13]. We refer the interested reader to our previous work [27] for details on this model.

As remarked in the introduction, we are primarily interested in valuing the far future for environmental problems rather than the short time discount of finance, the latter implying positive interest rates, while the former involves positive as well as negative rates. For this reason, we next review in more detail the Vasicek model allowing for both positive and negative rates than the CIR and log-normal models of which we only present a sketched review.

4.2. The Vasicek (Ornstein–Uhlenbeck) Model

In this model, the rates are described by the Ornstein–Uhlenbeck process [33], which is a diffusion model with linear drift and constant noise intensity:

$$dr\left(t\right)=-\alpha \left[r\right(t)-m]+kdW\left(t\right),$$

(65)

where $r\left(t\right)$ is the rate and $W\left(t\right)$ is the Wiener process. The parameter m (the normal level) is the mean value to which rates revert, $k>0$ is the amplitude of fluctuations, and $\alpha >0$ is the strength of the reversion to the mean. These parameters have to be estimated from empirical data.

In this case, the Fokker–Planck equation for the joint density $p(x,r,t|r_0)$ of the bidimensional process $\left(x\right(t),r(t\left)\right)$, given by Equation (17), reads

$$\frac{\partial p}{\partial t}=-r\frac{\partial p}{\partial x}+\alpha \frac{\partial }{\partial r}\left[(r-m)p\right]+\frac{1}{2}{k}^2\frac{{\partial }^{2}p}{\partial r^2},$$

(66)

with the initial condition

$$p(x,r,0|r_0)=\delta \left(x\right)\delta (r-r_0).$$

(67)

The joint Fourier transform of these equations results in a simpler initial-value problem for the joint characteristic function $\tilde{p}({\omega }_1,{\omega }_2,t|r_0)$, which can be readily solved to yield the Gaussian density [27]

$$\tilde{p}({\omega }_1,{\omega }_2,t|r_0)=exp\left\{-A({\omega }_1,t){\omega }_2^2-B({\omega }_1,t|r_0){\omega }_2-C({\omega }_1,t|r_0)\right\},$$

(68)

where $A({\omega }_1,t)$, $B({\omega }_1,t)$ and $C({\omega }_1,t)$ are given by [27]

$$A({\omega }_1,t)=\frac{k^2}{4\alpha }\left(1-e^{-2\alpha t}\right),$$

(69)

$$B({\omega }_1,t|r_0)=i{r}_0{e}^{-\alpha t}+\frac{k^2{\omega }_1}{2{\alpha }^2}\left(1-2e^{-\alpha t}+e^{-2\alpha t}\right)+im\left(1-e^{-\alpha t}\right),$$

(70)

and

$$\begin{array}{ccc}\hfill C({\omega }_1,t|r_0)& =& i{\omega }_1{r}_0\frac{1}{\alpha }\left(1-e^{-\alpha t}\right)+\frac{k^2{\omega }_1^2}{2{\alpha }^3}\left[\alpha t-2\left(1-e^{-\alpha t}\right)+\frac{1}{2}\left(1-e^{-2\alpha t}\right)\right]\hfill \\ & & +im{\omega }_1\left[t-\frac{1}{\alpha }\left(1-e^{-\alpha t}\right)\right].\hfill \end{array}$$

(71)

The characteristic function of the rate $r\left(t\right)$ is obtained by setting ${\omega }_2=0$ in Equation (68), which also results in the Gaussian density

$$\tilde{p}({\omega }_2,t|r_0)=exp\left\{-\frac{k^2}{4\alpha }\left(1-e^{-2\alpha t}\right){\omega }_2^2-i\left[r_0{e}^{-\alpha t}+m\left(1-e^{-\alpha t}\right)\right]{\omega }_2\right\},$$

(72)

and in the stationary state $(t\to \infty$), we have

$${\tilde{p}}_{st}\left({\omega }_2\right)=e^{-(k^2/4\alpha ){\omega }_2^2-im{\omega }_2}⟹p_{st}\left(r\right)={\left(\frac{\alpha }{\pi k^2}\right)}^{1/2}{e}^{-\alpha {(r-m)}^2/k^2},$$

(73)

which proves that the normal level m is the stationary mean value,

$$m=\mathbb{E}\left[r\right(t\left)\right].$$

(74)

It can also be shown that the correlation function of the process, defined as the average

$$C\left(\tau \right)=\mathbb{E}\left[r(t+\tau )r\left(t\right)\right]-\left[\mathbb{E}\right[r\left(t\right)]{]}^2,$$

($\tau \ge 0$) in the stationary state reads [27]

$$C\left(\tau \right)=(k^2/2\alpha )e^{-\alpha \tau },$$

(75)

which means that ${\alpha }^{-1}$ is the correlation time, ${\tau }_c$, of the rate. Indeed,

$${\tau }_c\equiv \frac{1}{C\left(0\right)}{\int }_0^{\infty }C\left(\tau \right)d\tau ={\alpha }^{-1}.$$

Let us observe that the volatility, ${\sigma }^2=C\left(0\right)$, is independent of the normal level and given by

$${\sigma }^2=k^2/2\alpha .$$

(76)

The discount function $D\left(t\right|r_0)$ is also obtained from Equations (68)–(71) although, in this case, after setting ${\omega }_1=-i$ and ${\omega }_2=0$ (cf. Equation (31)). We have

$$\begin{array}{ccc}\hfill lnD\left(t\right)& =& -\frac{r_0}{\alpha }\left(1-e^{-\alpha t}\right)\hfill \\ & +& \frac{k^2}{2{\alpha }^3}\left[\alpha t-2\left(1-e^{-\alpha t}\right)+\frac{1}{2}\left(1-e^{-2\alpha t}\right)\right]-m\left[t-\frac{1}{\alpha }\left(1-e^{-\alpha t}\right)\right],\hfill \end{array}$$

which, after rearranging terms, can be written as

$$lnD\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-e^{-\alpha t}\right)\right]\left(1-e^{-\alpha t}\right),$$

(77)

where $r_0=r\left(0\right)$ is the initial rate. Note that, as $t\to \infty$ (in fact when $t\gg {\alpha }^{-1}$, i.e., for times much greater than the correlation time ${\alpha }^{-1}$) Equation (77) shows at once that the discount function of the Vasicek model has the typical exponential decay

$$D\left(t\right)\simeq e^{-r_{\infty }t},$$

(78)

where

$$r_{\infty }=m-k^2/2{\alpha }^2,$$

(79)

is the long-run discount rate. Let us note that the long-run rate can be defined as the limit

$$r_{\infty }=-\underset{t\to \infty }{lim}\frac{lnD\left(t\right)}{t},$$

(80)

as long as the limit exists. Let us also note the important fact that $r_{\infty }$ is smaller than the mean value of the return given by the normal level m. This reduction is quantified by the “noise-to-signal” ratio $k/\alpha$, which means that either a long persistence (recall that this is equivalent to long correlation time, i.e., $\alpha$ small) or an increase of the noise fluctuations (i.e., k large) reduce the long-run discount rate as compared with the average rate m.

Finally, we easily see from Equation (77) that, as $t\to 0$, the discount function approximates to $D\left(t\right)\simeq e^{-r_{0}t}$, which would correspond to a fixed interest rate without random fluctuations or deterministic changes.

Risk Aversion

As mentioned above, risk aversion is taken into account by introducing the market price of risk $q\left(r\right)$ and changing the drift according to Equation (32). For the Vasicek model, in which $f\left(r\right)=-\alpha (r-m)$ and $g\left(r\right)=k$, we have

$$f^{*}\left(r\right)=-\alpha (r-m)+kq\left(r\right),$$

(81)

and assuming $q\left(r\right)=q$ to be a constant independent of r, we write

$$f^{*}\left(r\right)=-\alpha (r-m^{*}),$$

(82)

where

$$m^{*}=m+\frac{qk}{\alpha }.$$

(83)

Since the modified drift $f^{*}\left(r\right)$ has the same form that $f\left(r\right)$, we conclude that the adjusted-for-risk discount function will be given by Equation (77) after the replacement $m\to m^{*}$. In particular, the adjusted long-run discount now reads (cf. Equation (79))

$$r_{\infty }^{*}=m+\frac{qk}{\alpha }-\frac{k^2}{2{\alpha }^2}.$$

(84)

We thus see that the long-run discount depends on the historical rate m; however, this is shifted by two terms. The first term 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 $\alpha$). We rewrite Equation (84) as

$$r_{\infty }^{*}=m+\frac{k}{\alpha }\left(q-\frac{k}{2\alpha }\right).$$

(85)

This shows that the overall shift in the long-run discount rate will be positive or negative depending on the size of the market price of risk and on the noise-to-signal ratio between the volatility parameter and the reversion rate.

It is not surprising that the market price of risk raises the long term rate; however, it is not so obvious that uncertainty and persistence can lower it. Indeed, for any given mean interest rate m, by varying k and $\alpha$, the long-run discount rate $r_{\infty }$ can take on any value less than m, including negative values, while, at the same time, the standard deviation $\sigma$ can also be made to take on any arbitrary positive value.

A negative long-run rate is due to the amplification of negative real interest rates $r\left(t\right)$. 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 Vasicek model, and even taking into account risk aversion, the long-run discounting rate can be much lower than the mean and, indeed, can correspond to low interest rates that are rarely observed.

4.3. The Cox–Ingersoll–Ross (Feller) Model

In the financial literature, one of the most accepted models for interest rates is the Cox–Ingersoll–Ross (CIR) model [42] where rates follow the Feller process described by drift and noise intensity given, respectively, by [45]

$$f\left(r\right)=-\alpha (r-m),g\left(r\right)=k\sqrt{r}.$$

(86)

The Feller model is thus a diffusion process described by the stochastic differential equation

$$dr\left(t\right)=-\alpha \left[r\left(t\right)-m\right]dt+k\sqrt{r(t})dW\left(t\right),$$

(87)

where $W\left(t\right)$ is the standard Wiener process, and, as in the OU process, $m>0$ represents the mean stationary rate (the normal level), and ${\alpha }^{-1}$ is the correlation time [27]. Let us note that, in one-dimensional diffusions, the diffusion coefficient is given by the square of the noise intensity, and we thus see that the Feller process has a linear diffusion vanishing at the origin. This turns the origin into a singular boundary, which results in significant properties for the process [43].

As in the Vasicek model, the linear drift results in a restoring force, which, in the absence of noise, makes the process decay toward the normal level m. On the other hand, the state-dependent noise intensity $k\sqrt{r}$ for large values of r magnifies the effect of noise, while when r goes to zero, this effect vanishes. Therefore, as the process approaches the origin, the drift drags r towards m. Hence, since $m>0$, starting at some positive value $r_0>0$ the process cannot attain negative values with the overall result that the Feller process always remains positive.

Previous works [27,43] that we reviewed rather thoroughly presented the properties of the Feller process, and we refer the reader to these works for more detailed information. The process is not Gaussian, and the stationary PDF as $t\to \infty$ is the Gamma distribution [27]

$$p_{st}\left(r\right)=\frac{{(2\alpha /k^2)}^{\theta }}{\Gamma \left(\theta \right)}{r}^{\theta -1}{e}^{-(2\alpha /k^2)r},$$

(88)

where

$$\theta =\frac{2\alpha m}{k^2}$$

(89)

is a positive and dimensionless constant that combines all the parameters of the model into a single expression. As mentioned above, a major characteristic of the Feller process is that $r\left(t\right)$ cannot attain negative values, which makes the model a convenient tool for pricing bonds, which are never negative [13].

In the Feller model, the joint density of the discounting process $\left(x\right(t),r(t\left)\right)$ defined in Equation (16) obeys the Fokker–Planck equation (FPE) (cf. Equations (17) and (18))

$$\frac{\partial p}{\partial t}=-r\frac{\partial p}{\partial x}+\alpha \frac{\partial }{\partial r}\left[(r-m)p\right]+\frac{k^2}{2}\frac{{\partial }^2}{\partial r^2}\left(rp\right),$$

(90)

with the initial condition

$$p(x,r,0|r_0)=\delta \left(x\right)\delta (r-r_0).$$

(91)

The joint Fourier transform, Equation (30), turns Equations (90)–(91) into a more manageable problem:

$$\frac{\partial \tilde{p}}{\partial t}=\left({\omega }_1-\alpha {\omega }_2-i\frac{k^2}{2}{\omega }_2^2\right)\frac{\partial \tilde{p}}{\partial {\tilde{\omega }}_2}-i\alpha m{\omega }_2\tilde{p},$$

(92)

$$\tilde{p}({\omega }_1,{\omega }_2,0|r_0)=e^{-i{\omega }_2{r}_0}.$$

(93)

Equation (92) is a linear partial differential equation of first order whose solution can be obtained by the method of characteristics, and we refer the interested reader to our work [27] for a detailed information. Once we know the solution $\tilde{p}({\omega }_1,{\omega }_2,t|r_0)$, the discount function is then obtained through Equation (31) with the result [27]

$$D\left(t\right)={\left[\frac{2\lambda e^{-(\lambda -\alpha )t/2}}{(\lambda +\alpha )+(\lambda -\alpha )e^{-\lambda t}}\right]}^{\theta }exp\left\{-\frac{2(1-e^{-\lambda t})r_0}{(\lambda +\alpha )+(\lambda -\alpha )e^{-\lambda t}}\right\},$$

(94)

where $\theta$ is defined in Equation (89) and

$$\lambda =\sqrt{{\alpha }^2+2k^2}.$$

(95)

Notice that $\lambda >\alpha$ and the time scale represented by ${\lambda }^{-1}$ is smaller than the correlation time ${\alpha }^{-1}$.

In this case, the long-run discount rate, defined by the limit (cf. Equation (80))

$$r_{\infty }=-\underset{t\to \infty }{lim}\frac{lnD\left(t\right)}{t},$$

is directly obtained from Equation (94) with the result

$$r_{\infty }=\frac{1}{2}(\lambda -\alpha )\theta ,$$

(96)

and, as in the Vasicek model, the effective discount reduces to the expected exponential decay

$$D\left(t\right)\simeq e^{-r_{\infty }t}(t\to \infty ).$$

(97)

Substituting into Equation (96) the expressions for $\theta$ and $\lambda$ given in Equations (89) and (95), we write

$$r_{\infty }=\frac{2m}{1+\sqrt{1+2k^2/{\alpha }^2}},$$

(98)

which clearly shows that the long-run discount rate is always smaller than the stationary average rate:

$$r_{\infty }<m.$$

Figure 1 shows the discount function $D\left(t\right)$ along with the quantity $-lnD\left(t\right)/t$ (cf. Equation (80)) and compare them with the Vasicek model with equivalent parameters.

Risk Aversion

For the Feller process, the adjusted drift for risk defined in Equation (54) reads

$$f^{*}\left(r\right)=-\alpha (r-m)+kq\left(r\right)\sqrt{r},$$

(99)

where $q\left(r\right)$ is the market price of risk as discussed in the previous section. For any function $q\left(r\right)$ (including a constant market price of risk q), this adjusted drift leads to an unsolvable Fokker–Planck equation with no analytical expression for the adjusted discount and the long-run discount rate. It is, nonetheless, possible to obtain analytical expressions for these quantities if the market price of risk has the following functional form

$$q\left(r\right)=q\sqrt{r},$$

(100)

where $q\ge 0$ is a positive quantity. In such a case, we may write

$$f^{*}\left(r\right)=-{\alpha }^{*}(r-m^{*}),$$

(101)

where

$${\alpha }^{*}=\alpha -kq,m^{*}=\frac{\alpha m}{\alpha -kq}.$$

(102)

The adjusted drift has the same form as $f\left(r\right)$. Therefore, the adjusted discount function will be given Equation (94) with the replacements $\alpha \to {\alpha }^{*}$ and $m\to m^{*}$, and the long-run discount is (cf. Equation (98))

$$r_{\infty }^{*}=\frac{2m^{*}}{1+\sqrt{1+2k^2/{{\alpha }^{*}}^2}}.$$

(103)

From the definitions of ${\alpha }^{*}$ and $m^{*}$, we easily see that ${\alpha }^{*}\le \alpha$ and ${\alpha }^{*}{m}^{*}=\alpha m$. Hence, writing $r^{*}$ as

$$r_{\infty }^{*}=\frac{2{\alpha }^{*}{m}^{*}}{{\alpha }^{*}+\sqrt{{{\alpha }^{*}}^2+2k^2}}\ge \frac{2\alpha m}{\alpha +\sqrt{{\alpha }^2+2k^2}}=r_{\infty },$$

so that $r_{\infty }^{*}\ge r_{\infty }$, and, if the market price of risk has the form given in Equation (100), then, in the CIR model, risk always increases the long-run discount rate regardless of the noise intensity and persistence.

4.4. The Log-Normal Model

In this model, rates are described by the the geometric Brownian motion (log-normal process), and the model is determined by the stochastic differential equation

$$\frac{dr}{r}=\alpha dt+kdW\left(t\right),$$

(104)

where r is the interest rate, $\alpha$ and k are constant parameters. $\alpha$ may be positive or negative, whereas k is always positive, and $W\left(t\right)$ is the standard Wiener process. Equation (104) can be integrated at once yielding

$$r\left(t\right)=r_{0}exp\left\{\left(\alpha -\frac{k^2}{2}\right)t+kW\left(t\right)\right\},$$

(105)

showing that $r\left(t\right)$ is never negative ($r_0>0$). Therefore, the log-normal model is more suited for modeling nominal interest rates and bonds in finance than for the long-run real rates of environmental economics. Contrary to the OU and Feller processes, the log-normal process does not show reversion to the mean. Indeed, as t increases, we see from Equation (105) that the rate either diverges when $\alpha >0$ or goes to zero if $\alpha <0$. In an equivalent way, one can also show from Equation (105) that the mean and variance of the process are [27]

$$\langle r\left(t\right)\rangle =r_0{e}^{\alpha t},\mathrm{Var}\left[r\left(t\right)\right]=r_0^2{e}^{2\alpha t}\left(e^{k^{2}t}-1\right).$$

The discount associated with the log-normal process model was studied in 1978 by L. U. Dothan [44], and, in finance, it is sometimes refereed to as the Dothan model. As it allows for analytical treatment, it is one of the models used in the literature [13]. For this model, the FPE for the joint density of the discounting processes $\left(x\right(t),r(t\left)\right)$ is given by (cf. Equation (17))

$$\frac{\partial p}{\partial t}=-r\frac{\partial p}{\partial x}-\alpha \frac{\partial }{\partial r}\left(rp\right)+\frac{1}{2}{k}^2\frac{{\partial }^2}{\partial r^2}\left(r^{2}p\right),$$

(106)

with the usual initial condition given by Equation (18). The Fourier transform of this expression leads to the following equation for the characteristic function $\tilde{p}({\omega }_1,{\omega }_2,t|r_0)$

$$\frac{\partial \tilde{p}}{\partial t}=\left({\omega }_1+\alpha {\omega }_2\right)\frac{\partial \tilde{p}}{\partial {\tilde{\omega }}_2}+\frac{1}{2}{k}^2{\omega }_2^2\frac{{\partial }^2\tilde{p}}{\partial {\tilde{\omega }}_2^2}$$

(107)

and the initial condition (91). Equation (106) is a partial differential equation of second order, which cannot be solved by the method of characteristics, and we refer the interested reader to our work [27] for more information on how to solve Equation (107) using the time–Laplace transform. Hence—and contrary to Vasicek and CIR models where it is possible to obtain exact expressions for the discount function $D\left(t\right)$—for the log-normal case, we can only achieve the exact expression of its Laplace transform,

$$\widehat{D}\left(s\right)={\int }_0^{\infty }{e}^{-st}D\left(t\right)dt.$$

The resulting formula—written as an integral of special functions, the Kummer function—is rather intricate, and we will not write it here (see [27] for more information). However, from the exact expression for $\widehat{D}\left(s\right)$, we can obtain asymptotic expressions as $t\to \infty$ of the discount function $D\left(t\right)$ in real time. This is done using the so-called Tauberian theorems, which relate the small s behavior of $\widehat{D}\left(s\right)$ with the long-time behavior of $D\left(t\right)$ [46,47]. The final result is the following asymptotic expression for the discount function $D\left(t\right)$ in the long run as $t\to \infty$ [27]

$$D\left(t\right)\sim \left\{\begin{array}{cc}\mathrm{constant}\hfill & k^2/2>\alpha ,\hfill \\ e^{-r_{\infty }t}\hfill & k^2/2<\alpha ,\hfill \\ t^{-1/2}\hfill & k^2/2=\alpha .\hfill \end{array}\right.$$

(108)

The asymptotic form of the discount function thus depends on the values taken by the ratio $\alpha /k^2$ between the strength of the constant deterministic drift $\alpha$ and the amplitude of fluctuations given by $k^2/2$ (which can be considered the “signal-to-noise ratio” of this model).

(i) The case $k^2/2>\alpha$ corresponds to strong fluctuations, where the noise intensity $k^2/2$ is greater than the drift parameter $\alpha$. In this case, the discount tends to a constant value (for the actual value of this constant, see [27]).

(ii) The case $k^2/2<\alpha$ corresponds to mild fluctuations for which the deterministic drift is stronger than noise. In such a case, the discount function has the expected exponential decay

$$D\left(t\right)\sim e^{-r_{\infty }t},$$

(109)

with a long-run rate of discount given by [27]

$$r_{\infty }=\frac{1}{\delta }\left(\alpha -\frac{k^2}{2}\right),$$

(110)

where $\delta >1$ is a positive numerical factor that only depends on the ratio $2\alpha /k^2$ and reads

$$\delta =\psi \left(2\alpha /k^2\right)+\frac{1}{2\alpha /k^2-1},$$

(111)

where $\psi (·)$ is the digamma function.

Let us write Equation (109) in a more characteristic form. Indeed, from Equation (105), we see that

$$\mathbb{E}\left[ln\frac{r\left(t\right)}{r_0}\right]=\left(\alpha -\frac{k^2}{2}\right)t,$$

and, with the help of Equation (109), we write Equation (109) as

$$D\left(t\right)\sim exp\left\{-\frac{1}{\delta }\mathbb{E}\left[ln\frac{r\left(t\right)}{r_0}\right]\right\},$$

(112)

($t\to \infty$ and $k^2/2<\alpha$). Note that the average $\mathbb{E}[lnr\left(t\right)/r_0]$ is what a practitioner would take as an estimate of the discount rate up to time t within the log-normal model. Since $\delta >1$, the analytical result (112) shows that the actual long-run rate of the model is a fraction of the average rate. We indicated elsewhere that the long-run discount rate is at most 73% of the average rate [26].

In this way, when $2\alpha /k^2>1$, the log-normal model follows a similar pattern to that of the OU and Feller models: In all of them, the long-run rate is smaller than the average rate. This general statement is a direct consequence of Jensen’s inequality, which states that the average of a convex function is greater than or equal to the function of the average; that is, $\mathbb{E}\left[f\right(X\left)\right]\ge f\left(\mathbb{E}\right[X\left]\right)$. Assuming f to be the decreasing exponential and X as the cumulative process $x\left(t\right)$ defined in Equation (14), it follows immediately that the long-run rate $r_{\infty }$ must be always less than or equal to the average rate. Nonetheless, our procedure quantifies the difference among averages [27].

(iii) The critical case $\alpha =k^2/2$, in which deterministic motion and fluctuations are balanced, leads to the hyperbolic discount function as obtained by Farmer and Geanakoplos [48,49]. The hyperbolic $D\left(t\right)$ is substantially greater than any exponential decaying function, showing that there is no long-run rate of interest in this case. In fact, the long-run rate of interest is 0; however, that does not convey as precise information as saying that $D\left(t\right)$ is approximately $k/\sqrt{t}$ for all large t. Since the sum (i.e., the integral) of all these $D\left(t\right)$ is infinite, such $D\left(t\right)$ assigns infinite value to any permanent positive flow of consumption: the infinite future is infinitely valuable.

Risk Aversion

Let us very briefly comment on the inclusion of risk aversion in the Dothan model. For the log-normal process $f\left(r\right)=\alpha r$ and $g\left(r\right)=kr$ and

$$f^{*}\left(r\right)=[\alpha +kq\left(r\right)]r.$$

Assuming a constant market price of risk, $q\left(r\right)=q\ge 0$, we have

$$f^{*}\left(r\right)={\alpha }^{*}r,{\alpha }^{*}=\alpha +q.$$

Again, $f^{*}\left(r\right)$ has the same form than $f\left(r\right)$, and all previous results will apply after making the replacement $\alpha \to \alpha +q$.

5. Some Empirical Results

In order to choose an appropriate model for rates that would allow us to obtain realistic long-run discount functions, we performed a rather complete empirical study on interest rates combined with inflation. Our study follows the line partly initiated by Newell and Pizer [50] (see also [51]). To our knowledge, there are few empirical studies on real rates with some exceptions. We remark here the recent and excellent survey by Giglio et al. on the housing market in London and Singapore [52,53,54], which allowed for a rather realistic estimation of long-run discount rates.

Our first concern was knowing how the discount process depended on the underlying random process that characterizes interest rates. To this end, we collected data for the nominal interest rates and inflation of fourteen countries over time spans ranging from 87 to 318 years [26]. The countries in our sample are Argentina (ARG, 1864–1960), Australia (AUS, 1861–2012), Chile (CHL, 1925–2012), Germany (DEU, 1820–2012), Denmark (DNK, 1821–2012), Spain (ESP, 1821–2012), United Kingdom (GBR, 1694–2012), Italy (ITA, 1861–2012), Japan (JPN, 1921–2012), Netherlands (NLD, 1813–2012), Sweden (SWE, 1868–2012), the United States (USA, 1820–2012) and South Africa (ZAF, 1920–2012). The data are summarized in

Table 2. List of the countries analyzed. CPI stand for Consumer Price Index. Data has different specificities, particularly in terms of empty records as has been reported elsewhere [26,29]. *, We have taken the discount (ID) rate since the government bond yield data was not available.

Country	CPI	Bond Yield	From	To	Records
Italy	CPITAM	IGITA10	12/31/1861	09/30/2012	565
		annual from 12/31/1861	quarterly
		quarterly from 12/31/1919
Chile	CPCHLM	IDCHLM	03/31/1925	09/30/2012	312
		quarterly	quarterly
Canada	CPCANM	IGCAN10	12/31/1913	09/30/2012	357
		quarterly	quarterly
Germany	CPDEUM	IGDEU10	12/31/1820	09/30/2012	729
		annual from 12/31/1820	quarterly
		quarterly from 12/31/1869
Spain	CPESPM	IGESP10	12/31/1821	09/30/2012	709
		annual from 12/31/1821	quarterly
		quarterly from 12/31/1920
Argentina	CPARGM	IGARGM	12/31/1864	03/31/1960	342
		annual from 12/31/1864	quarterly
		quarterly from 12/31/1932
Netherlands	CPNLDM	IGNLD10D	12/31/1813	12/31/2012	189
		annual	annual
Japan	CPJPNM	IGJPN10D	12/31/1921	12/31/2012	325
		quarterly	quarterly
Australia	CPAUSM	IGAUS10	12/31/1861	09/30/2012	564
		annual from 12/31/1861	quarterly
		quarterly 12/31/1991
Denmark	CPDNKM	IGDNK10	12/31/1821	09/30/2012	725
		annual from 12/31/1821	quarterly
		quarterly from 12/31/1914
South Africa	CPZAFM	IGZAF10	12/31/1920	09/30/2012	329
		quarterly	quarterly
Sweden	CPSWEM	IGSWE10	12/31/1868	09/30/2012	135
		annual	annual
United Kingdom	CPGBRM	IDGBRD ${}^{*}$	12/31/1694	12/31/2012	309
		annual	annual
United States	CPUSAM	TRUSG10M	12/31/1820	10/30/2012	183
		annual	annual

.

Since all but two of our nominal interest rate processes are for 10-year government bonds, which pay out over a 10-year period, we smoothed out inflation rates with a 10-year moving average and subtracted the annualized inflation index from the annualized nominal rate to compute the real interest rate, as explained in the previous section by means of the Fisher’s procedure (cf. Equation (61)),

$$r\left(t\right)=n\left(t\right)-i\left(t\right),$$

where $n\left(t\right)$ is the nominal rate and $i\left(t\right)$ is the inflation rate. The particular case of the United States is plotted in Figure 2.

In our empirical analysis, the nominal rates are determined by IG rates constructed from the 10-year Government Bond Yield $y(t,\tau )$ with $\tau$= 10 years. Thus, looking at Equations (63) and (64), we estimate the nominal rates by

$$n\left(t\right)\sim y(t,\tau =10 \mathrm{years}).$$

Let us recall that denoting, by $B(t,t+\tau )$, the government bond issued at time t and maturing at time $t+\tau$ with unit maturity, $B(t,t)=1$, the yield $y(t,\tau )$ is defined as (cf. Equation (38))

$$y(t,\tau )\equiv -\frac{1}{\tau }lnB(t,t+\tau )⟹B(t,t+\tau )=e^{-\tau y\left(t\right|\tau )}.$$

One can argue that $\tau =10$ years is not a short period of time in order to consider $y(t,\tau =10 \mathrm{years})$ a very accurate estimator of $n\left(t\right)$ (cf. Equations (63) and (64)). Although this may be true, we must bear in mind that 10 year bonds are the shortest bonds available for most of the countries analyzed.

The inflation rate is estimated through the Consumer Price Index (CPI) as

$$i\left(t\right)\sim \frac{1}{\tau }ln[I(t+\tau )/I\left(t\right)],$$

where $I\left(t\right)$ is the aggregated inflation up to time t, and $\tau =$10 years. The relation between $I\left(t\right)$ and the Consumer Price Index (CPI) is

$$I(t+\tau )=I\left(t\right)\prod _{j=0}^{\tau -1}\left[1+C(t+j)\right],$$

where $C\left(t\right)$ is the time series of the empirical CPI. The instantaneous rate of inflation $i\left(t\right)$ is, therefore, estimated by the quantity $i(t+\tau )$, which is written in terms of the CPI reads

$$i\left(t\right)\sim i(t+\tau )=\frac{1}{\tau }\sum _{j=0}^{\tau -1}ln\left[1+C(t+j)\right].$$

A remarkable characteristic observed for all countries is that real interest rates frequently become negative as the real interest rates are mostly dominated by inflation $i\left(t\right)>0$ (see Figure 2). In some cases, as we can see in

Table 3. The OU (Vasicek) model parameter estimation in yearly units using stationary averages. “Neg RI” provides the time percentage and the number of years with negative real interest rates. The columns $\widehat{m}$, $\widehat{k}$ (in %) and $\widehat{\alpha }$ are estimates from the country time series; ${\widehat{r}}_{\infty }$ (in %) is evaluated from Equation (79). The Min and Max columns give reasons regarding the level of robustness of the estimation as they provide the minimum and the maximum values of the parameter estimation for four data blocks of equal length. The parameter $\alpha$ is estimated by fitting the empirical correlation function to an exponential (cf. Equation (75)) after using the whole data block. Countries in boldface are those considered historically more stable with positive long-run rates ${\widehat{r}}_{\infty }>0$.

Country	Neg RI	$\widehat{\mathit{m}}$	Min	Max	$\widehat{\mathit{k}}$	Min	Max	$\widehat{\mathit{\alpha }}$	${\widehat{\mathit{r}}}_{\mathbf{\infty }}$
Italy	$28\%\left(40y\right)$	$-0.3$	$-9.1$	$5.6$	$6.9$	$0.8$	$10.1$	$0.22$	$-5.4$
Chile	$56\%\left(43y\right)$	$-6.8$	$-20.2$	$12.0$	$25.2$	$5.6$	$44.1$	$0.40$	$-26$
Canada	$22\%\left(20y\right)$	$\mathbf{2.9}$	$0.1$	6	$\mathbf{2.3}$	$1.1$	$2.0$	$0.26$	$\mathbf{2.5}$
Germany	$14\%\left(25y\right)$	$-10.7$	$-51.0$	$4.0$	$33.9$	$0.9$	$61.4$	$0.20$	$-160$
Spain	$25\%\left(45y\right)$	$5.7$	$-0.5$	$13.5$	$2.9$	$1.2$	$3.6$	$0.06$	$-6.4$
Argentina	$20\%\left(17y\right)$	$\mathbf{2.4}$	$-2.9$	$6.8$	$\mathbf{6.2}$	$2.8$	$6.7$	$0.39$	$\mathbf{1.1}$
Netherlands	$17\%\left(33y\right)$	$\mathbf{3.2}$	$0.8$	$5.4$	$\mathbf{1.6}$	$0.8$	$2.2$	$0.14$	$\mathbf{2.4}$
Japan	$33\%\left(26y\right)$	$-2.2$	$-7.8$	$4.0$	$9.7$	$1.1$	$13.2$	$0.24$	$-10$
Australia	$23\%\left(33y\right)$	$\mathbf{2.6}$	$-0.7$	$4.9$	$\mathbf{2.3}$	$0.7$	$2.8$	$0.19$	$\mathbf{1.9}$
Denmark	$18\%\left(33y\right)$	$\mathbf{3.2}$	$1.5$	$4.3$	$\mathbf{2.3}$	$1.1$	$2.9$	$0.23$	$\mathbf{2.7}$
South Africa	$43\%\left(36y\right)$	$\mathbf{1.8}$	$-2.2$	$5.5$	$\mathbf{2.5}$	$1.2$	$2.0$	$0.21$	$\mathbf{1.1}$
Sweden	$28\%\left(38y\right)$	$\mathbf{2.3}$	$-0.3$	$3.9$	$\mathbf{2.5}$	$0.6$	$3.4$	$0.25$	$\mathbf{1.9}$
United Kingdom	$14\%\left(45y\right)$	$\mathbf{3.3}$	$1.4$	$4.3$	$\mathbf{1.9}$	$1.0$	$2.4$	$0.19$	$\mathbf{2.8}$
United States	$31\%\left(36y\right)$	$\mathbf{2.6}$	$1.0$	$4.0$	$\mathbf{1.8}$	$1.2$	$2.1$	$0.18$	$\mathbf{2.1}$
Stable countries	$23\%\left(33y\right)$	$\mathbf{2.7}$	$-0.14$	$5.0$	$\mathbf{2.6}$	$1.04$	$2.94$	$0.23$	$\mathbf{2.1}$
Unstable counntries	$31\%\left(36y\right)$	$-\mathbf{2.9}$	$17.7$	$1.8$	$\mathbf{16}$	$1.9$	$26.5$	$0.22$	$-\mathbf{42}$

(see also Figure 2), negative real rates show high frequency and long periods of time, and, on average, real interest rates are negative one quarter of the time.

This makes the Feller and log-normal models—as well as any other model assuming positive interest rates [13]—less interesting or at least less appropriate to model real interest $r\left(t\right)=n\left(t\right)-i\left(t\right)$ instead of solely considering nominal rates $n\left(t\right)$. It is, however, necessary to remember the fact that nominal rates can indeed become negative as has recently been observed in Western economies. We, therefore, confined the empirical work to the OU (Vasicek) model and then assumed the Local Expectation Hypothesis [36,37,38], according to which, we live in a risk neutral world, and the market price of risk is zero. Let us recall, as explained in Section 3, that the market price of risk $q=q(r,t)$ may be any function of the rate and time. There is, hence, no unique expression for it. Thus, in Section 4, we presented several expressions of the long-run rate, which include risk in different forms for all market models analyzed. The usual assumption in the literature [33,38] is that the price of risk is a constant that is independent of time and the value of the rate but without any empirical justification. This is a sensitive issue since data is quite scarce, particularly in environmental applications, for obtaining a credible estimation of q. Moreover, to our knowledge, in environmental problems, the estimations of the long-run rates do not take into account, nor even mention, the market price of risk [14,15,16,50,54]. In any case, we do not lessen the importance of taking into account some kind of risk in estimating log-run rates; however, unfortunately, with the data available to us, we cannot make any reliable estimation of q. For this reason, we have not taken into account the market price of risk, assuming risk-neutral investors and following the Local Expectation Hypothesis. In any case, the question is under consideration).

We can estimate the parameters m, $\alpha$ and k of the Vasicek model to each of the data series. There are several possible procedures. One of the possible methods is to deal with stationary averages. The parameter m can be estimated through the stationary mean value of the rate (cf. Equation (74))

$$m=\mathbb{E}\left[r\right(t\left)\right].$$

Parameters $\alpha$ and k can be estimated via the correlation function of the Ornstein–Uhlenbeck process. Thus, from Equation (75), we have

$$C(t-t^{\prime })=\frac{k^2}{2\alpha }{e}^{-\alpha |t-t^{\prime }|}.$$

The empirical correlation can then be fitted by an exponential, which in turn allows us to estimate $\alpha$ (measured in 1/year units) for each country. The parameter k is obtained from the empirical standard deviation ${\sigma }^2={\mathbb{E}\left[\right|r\left(t\right)-m|}^2]$, and for the Vasicek model, it is given by

$$k=\sigma \sqrt{2\alpha }.$$

The resulting parameters are shown in Table 3. The minimum and maximum values for each country allows us to show that parameters may indeed fluctuate over different periods of time.

Finally, the long-run discount rate can be evaluated from Equation (79),

$$r_{\infty }=m-k^2/2{\alpha }^2.$$

For this calculation, we neglected the market price of risk as mentioned above.

The countries studied can be divided into two groups. Nine countries have long-run positive rates (boldface in Table 3). The average historical rate for these nine countries is $\overline{m}=2.7$% while the average long-run rate is ${\overline{r}}_{\infty }=2.1$%, which, on average, is 29% lower than $\overline{m}$. Five countries with less stable behavior have long-run negative rates and an exponentially increasing discount.

Four cases of this group have a negative average rate m due to at least one period of runaway inflation; the exception is Spain, which has a (highly positive) mean real interest rate but still has a long-run negative rate. Convergence in this case to the long-run rate happens within 30 years and typically within less than a decade. This contrasts with other treatments, which assume that short term rates are always (or nearly always) positive and predict that the decrease in the discounting rate happens over a much longer timescale, which can be measured in hundreds of years [50,51,55,56,57,58].

Alternatively, we can estimate parameters using the well-established maximum likelihood procedure. For the Vasicek model, the maximum likelihood estimation is extensively documented in the financial mathematics literature (see for instance [13]). The approach differs from the previous one as it focuses attention on two consecutive steps of our time series (generally consecutive years) and takes the conditional probability to perform the estimation.

Table 4. Maximum likelihood estimation of the long-run interest rate for the Vasicek model. $\widehat{m}$ estimates of the mean real interest rate in 1/years (in %). $\widehat{\alpha }$ estimates the characteristic reversion time in 1/years. The squared root of $\widehat{{\mathit{k}}^2}$ is given in terms of $1/{\left(\mathrm{year}\right)}^3$ (multiplied by ${10}^4$ to be comparable with the results in Table 3). These estimators are accompanied by the square root of the variance of each estimator. ${\widehat{r}}_{\infty }$ estimates the long-run real interest rate with 1/year (in %). Negative values of ${\widehat{r}}_{\infty }$ imply that the discount function is asymptotically increasing. The standard error is obtained through error propagation. The last two rows show the average over all countries with the more stable countries ($r_{\infty }>0$) and the less stable countries ($r_{\infty }<0$). The error provided corresponds to the standard deviation of the ${\widehat{r}}_{\infty }$ for the different countries.

Country	$\widehat{\mathit{m}}$	${\mathit{\sigma }}_{\widehat{\mathit{m}}}$	$\widehat{\mathit{\alpha }}$	${\mathit{\sigma }}_{\widehat{\mathit{\alpha }}}$	$\widehat{{\mathit{k}}^2}$	${\mathit{\sigma }}_{\widehat{{\mathit{k}}^2}}$	${\widehat{\mathit{r}}}_{\mathbf{\infty }}$	${\mathit{\sigma }}_{{\widehat{\mathit{r}}}_{\mathbf{\infty }}}$
Italy	1.97	15.95	0.0056	0.0089	0.1146	0.068	−177.8	19.2
Chile	−5.79	31.46	0.0201	0.0227	31.07	2.49	−391.7	44.2
Canada	2.66	3.91	0.0142	0.0178	0.275	0.021	−4.15	3.94
Germany	−9.45	66.95	0.0071	0.0089	41.72	2.19	−4094	228
Spain	6.71	6.92	0.0167	0.0137	2.371	0.126	−35.78	7.28
Argentina	3.15	7.09	0.0228	0.0231	2.240	0.171	−18.31	7.27
Netherlands	5.99	0.78	0.1648	0.0550	1.797	0.243	5.66	0.78
Japan	5.02	24.68	0.0053	0.0114	1.396	0.109	−243.1	31.4
Australia	3.97	4.50	0.0089	0.0112	0.223	0.013	−10.29	4.58
South Africa	2.69	4.72	0.0154	0.0193	0.435	0.034	−6.49	4.77
Sweden	2.79	1.66	0.0676	0.0317	1.692	0.206	0.95	1.67
Denmark	4.10	2.59	0.0161	0.0133	0.315	0.017	−1.97	2.61
United Kingdom	3.42	0.62	0.1635	0.0326	3.137	0.253	2.83	0.62
United States	3.19	1.23	0.0603	0.0257	1.003	0.105	1.81	1.24
Stable countries	3.85	1.07	0.1140	0.0362	1.907	0.202	2.81	1.08
Unstable countries	1.50	16.86	0.0132	0.0150	8.120	0.523	−498.4	35.3

shows that the most inaccurate estimator is $\widehat{\alpha }$, an unsurprising fact since the estimation of $\alpha$ is known to be quite difficult for the Vasicek model [59]. The last two columns in Table 4 include the long-run interest rate estimator ${\widehat{r}}_{\infty }$ and its error calculated through error propagation.

Only four countries (the Netherlands, Sweden, the United Kingdom and the United States) show a positive long-run rate, $r_{\infty }>0$. This estimation procedure leads to more negative $r_{\infty }$. This feature can be attributed to the fact that, in most of the countries, estimating $\alpha$ via the maximum likelihood brings smaller values, which in turn leads to more drastic corrections to the long-run rate as $r_{\infty }$ is inversely proportional to $\alpha$ (remember that $r_{\infty }=m-k^2/2{\alpha }^2$). This effect is particularly relevant in most turbulent countries during last century (e.g., Germany) thus signaling a more intense lack of stationarity in empirical data. The averaged $r_{\infty }$ over all countries estimated via maximum likelihood is also sensitively smaller.

However, if we focus the attention on stable countries (with $r_{\infty }>0$) both estimation procedures bring quite similar results (see, for instance, the United States case in Table 3 and Table 4, 2.1% versus 1.8%). As in the previous estimation procedure, we also neglected the effects of risk aversion and the market price of risk.

The Vasicek model is therefore the only one among the three classic models allowing for negative rates, and for this reason, both the Feller and the log-normal models have been excluded from our analysis. However, for the Cox–Ingersoll–Ross (Feller) model, it is possible to redefine the model by shifting the process $y=r-r_{min}$ where $r_{min}<0$.

The estimation through the maximum likelihood procedure and its error analysis is then possible [59], and Figure 1 includes the shifted Cox–Ingersoll–Ross discount and compares it with the equivalent result assuming the Vasicek model. We demonstrated in Ref. [29] how to redefine the Feller process and how maximum likelihood estimation could be possible.

6. Discussion

We reviewed one of the most important aspects of economics and finance, i.e., the problem of discount, which weights the future relative to the present. The problem is clearly very relevant in finance over relatively short time spans; however, it is even more crucial for long-run planning in addressing environmental problems on how to act now with measures to mitigate the effects of climate change.

To our knowledge, this is a rather unknown issue to the econophysics community, and this review is particularly intended for this group. We thus addressed the problem with a simple approach and yet with a high level of rigor and generality. In this way, we also developed the traditional method used in mathematical finance to address the problem, i.e., the Feynman–Kac approach. In addition, we reviewed the bond pricing theory and its close similarity with discounting and presented a short introduction to the term structure of interest rates along with the market price of risk.

We obtained quantitative results on the problem by studying, in some detail, three standard models for the dynamical evolution of rates. These models are based on the Ornstein–Uhlenbeck process (the Vasicek model), thus, allowing for both positive and negative rates and also on the Feller and log-normal processes for positive rates. We presented the exact results for the discount function and asymptotic expressions as $t\to \infty$ leading to the long-run discount rate, and we discussed the modifications of these expressions when the market price of risk is taken into account.

An important conclusion is that, for all models, the long-run discount rate is always less than the long-time average rate. This is a conclusion that necessarily has to have consequences in any long-run economic planning. Finally, we reviewed our recent empirical study on 14 different countries, which obtained numerical values for the parameters that appear in the Vasicek model. We demonstrated two different estimation procedures and briefly discussed their differences and similarities.