On the social cost of carbon and discounting in the DICE model

1.   Introduction

To quantify the damages from anthropogenic emissions of heat-trapping greenhouse gases, economists model the dynamics of climate–socio-economic interactions using integrated assessment models (IAMs) that incorporate mathematical models of phenomena from geophysical science, economics, and social sciences. A central role for IAMs is to estimate the social cost of carbon (SCC), defined as the dollar value of the economic damage caused by a one metric ton increase in carbon dioxide (CO₂) emissions to the atmosphere. The SCC can then form the basis of a Pigouvian tax on greenhouse gas emissions, guide carbon prices in emissions trading schemes, inform cost-benefit analyses driving climate policy, or be used by investors looking to consider potential emissions costs in investment decisions. This wide range of uses highlights the importance of methods for estimating the SCC.

A critical and much debated issue in computing the SCC is the choice of discount rate. CO₂ is a long-lived greenhouse gas (GHG), with CO₂ emissions having a substantial impact on climate for thousands of years; e.g., 20–40% of a CO₂ perturbation remains in the atmosphere even after equilibration with the ocean on a timescale of 2–20 centuries ^[1,2]. Comparing costs and benefits over such long time horizons necessitates a careful choice of discount rates and methodologies in order to provide a fair comparison of costs paid now—via emissions reductions—with economic benefits from curbed global warming in the far distant future.

Conceptually, the most common approach to computing the SCC using IAMs is to compute baseline emissions and economic consumption. These computations are then repeated with an increment of CO₂ emissions in a chosen year, yielding reduced economic consumption into the future as the result of a climate damages function in the IAM. The results of these baseline and perturbed scenarios are then used to compute the SCC as the net present value (NPV) of the "damage stream" (i.e., pathway of lost consumption) (see, for example, ^[3,4]). This NPV approach is employed directly by two of the most widely used IAMs: PAGE (Policy Analysis of the Greenhouse Effect ^[5]) and FUND (Climate Framework for Uncertainty, Negotiation, and Distribution ^[6]). While descriptively this is what is intended with DICE (Dynamic Integrated model of Climate and Economy; see Figure 1 and associated discussion by Nordhaus W ^[4]), in practice the publicly available implementation of DICE uses a different method, based on shadow prices or Lagrange multipliers obtained as a result of solving an optimal control problem (^[7]). Hence, a critical question, answered herein, is whether or not these different computational approaches are in fact equivalent.

The classical Ramsey formula for social discount rate is given by $r = \rho + \theta g$, where $\rho$ is the pure rate of time preference (or "impatience factor"), $\theta$ is the elasticity of marginal utility of consumption, and $g$ is the growth rate of per-capita consumption ^[8]. While the calculation of the SCC as the NPV of the damage stream due to CO₂ pulse implies use of a discount rate, there is widespread disagreement in the literature regarding the appropriate discount rate to use in this context. For example, in some representative papers, the SCC is described as being computed as the NPV of a damage stream variously discounted at (a) the classical Ramsey formula (i.e., assuming a constant growth rate of consumption), (b) a time-varying Ramsey formula (reflecting the growth rate of time-varying consumption), (c) the pure rate of time preference, and (d) a fixed rate without any explicit reference to the Ramsey formula; see the studies by Marten AL, et al. ^{[3,4,5,6,9,10,11,12,13,14,15,16]}, and also the study by Tol RSJ^[17], which proposes exercises based on both (b) and (d). The report Economides 2018 provides a chapter on various discounting approaches in the economics of climate change. It is worth noting that in some cases it is unclear precisely how the discounting is being performed as it is described in words rather than by providing an explicit mathematical expression.

As briefly mentioned above, rather than computing the NPV of the damages stream, a conceptually distinct approach to computing the SCC relies on maximizing a social welfare function. As explicitly defined in the study of Newbold S C, et al. ^[19], "The [SCC] in a particular year is the decrease in aggregate consumption in that year that would change the current expected value of social welfare by the same amount as a one unit increase in carbon emissions in that year." Expressed formally, the SCC at time $t$ is defined as

where a simple manipulation shows that this is the ratio of the marginal prices of consumption and emissions. In fact, this is the computation contained in the code implementing the DICE model and optimal control problem as presented in the study by D. Nordhaus W, et al. ^[7]. Here, $W$ denotes social welfare, $E$ denotes emissions, and $C$ denotes consumption, whose definitions are made precise in Section 2.

The main contribution of this paper is to discuss and compare the three possible computational methods of the SCC indicated in the discussion above. In particular, we discuss a direct discretization of (1.1), an approach using the ratio of marginals obtained from solving an optimal control problem, and an NPV approach based on experimental CO₂ emission pulses. For the last approach, there are several options for how to perform the required discounting. The three approaches are compared theoretically, in terms of their accuracy and their computational complexity, and numerically, based on simulations of the DICE model. In addition, the three approaches and their analysis lead us to the following question: Which (if any) of the above options for computing the SCC are consistent with each other? Importantly, neither FUND nor PAGE involve an optimal control problem, meaning the first two methods mentioned here are not relevant. Hence, answering this consistency question is critical for comparing results between IAMs based on optimal control methodologies with those based on scenario methodologies.

Our main mathematical result, which covers the third computation method, is summarized in Theorem 1. The result demonstrates that in order to obtain consistency between the NPV calculation and the SCC values calculated via optimal control approaches, it is necessary to use a time-varying discount rate in the NPV calculation. In particular, the time-varying component is inversely proportional to economic growth, leading to a declining discount rate ^[20,21,22]. Specifically, returning to the Ramsey formula $r = \rho + \theta g$, mentioned above, the quantity $g$ is time-varying and given by the growth rate of consumption. The quantities $\rho$ and $\theta$ are "exogenous" to the model and need to be selected by the modeler on the basis of data and/or economic expertise. Note that it has been argued that the use of declining rates is justified in the presence of significant uncertainty in future economic growth (see ^[23]). Here, we see that the use of a declining discount rate is required for consistency with outputs from calculations based on shadow prices obtained when solving optimal control problems.

This paper is organized as follows. In Section 2 we succinctly present the DICE model as described in the study by Kellett C M, et al. ^[24], with model parameters provided in Appendix E. In Section 3, we investigate three methods for computation of the SCC using the DICE model. Section 3.1 presents the calculation based on a discretization of the difference quotient of the two quantities in (1.1). In Section 3.2 we demonstrate how to rewrite the DICE optimal control problem so that emissions and consumption are decision variables, which allows us to obtain the marginals or dual variables directly as an output of solving the optimal control problem. This has the computational benefit of removing the welfare function as an intermediary in the SCC calculation. In Section 3.3 we present the main mathematical result of the paper (Theorem 1), which shows an equivalence between calculating the SCC as a ratio of the dual variables versus the NPV when the discounting approach is computed via the time-varying Ramsey formula accounting for a time-varying growth rate of consumption. Section 4 presents some numerical results and discussion, and we conclude in Section 5. Mathematical necessities and the proof of our main result are collected in Appendices.

Notation: The set $\mathbb{N}_{0} = \{0, 1, 2, \ldots\}$ denotes the natural numbers including zero, while $\mathbb{N} = \{1, 2, \ldots\}$ denotes the natural numbers without zero. The real numbers are denoted as $\mathbb{R}$, and $\mathbb{R}_{\geq 0} = \{x\in \mathbb{R}| \ x\geq 0\}$ denotes the real numbers bigger or equal to zero. For $n\in \mathbb{N}$, the set $\mathbb{R}^n = \mathbb{R}\times \cdots \times \mathbb{R}$ denotes the $n$-fold Cartesian product of the real numbers, and $\mathit{\boldsymbol{x}}\in \mathbb{R}^n$ denotes an $n$-dimensional vector with $x(i)\in \mathbb{R}$, $i\in\{1, \ldots, n\}$. To simplify expressions and to define convergence rates we use the big $O$-notation, i.e., two functions $f_1, f_2:\mathbb{R}\rightarrow \mathbb{R}$, $f_2(x) > 0$ for all $x > 0$, satisfy $f_1(x) = O(f_2(x))$ if and only if there are $C > 0$ and $\delta > 0$ such that $|f_1(x)| \le Cf_2(x)$ for all $0 < x < \delta$. For a continuously differentiable function $f:\mathbb{R}^n \times \mathbb{R}^m \rightarrow \mathbb{R}$, $(x, y) \mapsto f(x, y)$, the gradient is denoted by $\nabla f(x, y)$. The gradient and the second derivative of $f$ with respect to a subset of variables, that is with respect to $x$, for example, are denoted by $\nabla_x f(x, y)$ and $\nabla_{xx}^2f(x, y)$, respectively.

2.   The DICE model

For the computation of the social cost of carbon, we consider the Dynamic Integrated model of Climate and Economy (DICE). For a comprehensive introduction of the model and its parameters we refer to the studies ^[4,24]. Here, we only give a brief introduction of the dynamic equations and their interpretation. The presentation follows the notation introduced in ^[24]. We focus on the parameters defining the R-DICE2016 model, but the results derived in this paper also hold for the parameters defining the DICE-2013R model. Following the R-DICE2016 setting we denote time in years by $t \in \mathbb{R}_{\geq 0}$ so that

with initial time $t_0 = 2015$, a sample period $\Delta = 5$, and $i\in \mathbb{N}_{0}$. The DICE model is described through the dynamics

modeling the carbon cycle (CAR), the climate (CLI), the capital (CAP), the population (POP), the total factor productivity (TFP), and the emissions intensity (EI). In (CAR), $M = [M_{\text{AT}}, M_{\text{UP}}, M_{\text{LO}}]^\top$ is a three dimensional vector describing the average mass of carbon in the atmosphere, in the upper ocean, and in the deep or lower ocean. The climate dynamics (CLI) capture the evolution of the temperature $T = [T_{\text{AT}}, T_{\text{LO}}]^\top$ in the atmosphere ($T_{\text{AT}}$) and in the lower ocean ($T_{\text{LO}}$). The one dimensional state $K$ represents the economic output in the economic dynamics (CAP). For the world-wide population $L$, the total factor productivity $A$, and the emissions intensity $\sigma$, the DICE model uses the Hassell model ^[25], a logistic-function model, and a monotonically decreasing function, respectively, to capture the dynamics.

The remaining time dependent quantities are described through the equations

and denoted as the mitigation effort $\theta_1$, the effect of greenhouse gases $F_{\rm{ EX}}$, and the emissions due to land use $E_{\rm{ Land}}$.

Control inputs are the savings rate, $s(i) \in [0, 1]$ for all $i \in \mathbb{N}_0$, and the mitigation rate, $\mu(i) \in [0, 1]$ for all $i \in \mathbb{N}_0$. An additional parameter, $E_{\rm{ Pulse}}(i) \in \mathbb{R}$, for all $i \in \mathbb{N}_0$, is used for the calculation of the partial derivatives needed contained in the definition of the SCC, cf. equation (3.2), below. The additional parameter $E_{\rm{ Pulse}}(i)$ has an impact on the total emissions

at time $i\in \mathbb{N}$, which we can identify in (CAR). The parameter is not present in the baseline DICE model where $E_{\rm{ Pulse}}(i) = 0$ for all $i\in \mathbb{N}_0$. We collect the states in the variable $X$, i.e.,

and use $X(0) = X_0$ to denote a fixed initial condition.

The control inputs are used to maximize the social welfare depending on the consumption and the utility at time $i\in \mathbb{N}_0$ as

Here, $C_{\rm{ Pulse}}(i) \in \mathbb{R}$, for all $i \in \mathbb{N}_0$, plays the same role as $E_{\rm{ Pulse}}(i)$ in Equation (CAR), manipulating the consumption in a specific year. Then, for a fixed horizon, $N\in \mathbb{N}$, the social welfare is defined as

and the optimal control problem (2.9) which maximizes the social welfare and generates optimal pathways is defined as

The optimal control problem (2.9) constitutes a nonlinear program (NLP). We use the notation $W^*$ to denote the optimal value with respect to $C_{\rm{Pulse}}(j) = E_{\rm{Pulse}}(j) = 0$, for all $j\in \{0, \ldots, N-1\}$, i.e., $W^*$ denotes the nominal or baseline solution of (2.9). In the case that $C_{\rm{Pulse}}(j) \neq 0$ or $E_{\rm{Pulse}}(j) \neq 0$ for $j\in \{0, \ldots, N-1\}$ is considered in (2.9), the optimal solution is denoted by $W^*[C_{\rm{Pulse}}(j)]$ or $W^*[E_{\rm{Pulse}}(j)]$, respectively. Similarly, optimal pathways are denoted using $\cdot^*$. In this context, $C^*(\cdot)$, $C^*[C_{\rm{Pulse}}(j)](\cdot)$, and $C^*[E_{\rm{Pulse}}(j)](\cdot)$ denote the optimal consumption pathways obtained through (2.9) with $C_{\rm{Pulse}}(j) = E_{\rm{Pulse}}(j) = 0$, $C_{\rm{Pulse}}(j)\neq 0$, and $E_{\rm{Pulse}}(j)\neq 0$, respectively (in all cases we set $C_{\rm{Pulse}}(i) = E_{\rm{Pulse}}(i) = 0$ for all $i\ne j$). Additionally, $N\in \mathbb{N}$ denotes a fixed horizon and $\mathit{\boldsymbol{s}}, \mathit{\boldsymbol{\mu}} \in \mathbb{R}^{N}$ define a short-hand notation for

For completeness, the parameters used in R-DICE2016, and in the simulations presented in this paper, are reported in Appendix E.

3.   Calculation of the social cost of carbon

In this section we take as a formal definition of the SCC at time $j\in \{0, \ldots, N-1\}$, a scaled version of the change of the optimal consumption $C^*(j)$ with respect to a change of the optimal emissions $E^*(j)$. Accordingly, the SCC in US dollars per kilogram of emitted CO₂ is defined as the change $-\frac{\partial C^*(j)}{\partial E^*(j)}$. However, instead of kilograms of emitted CO₂, it is common to express the SCC in terms of tonnes of CO₂, i.e., we define

(see [24, Sec. 2.8]). Since $E^*(j)$ is itself an outcome of the optimal control problem (2.9) rather than an independent variable on which $C^*(j)$ would depend, it is not immediately clear how this partial derivative has to be interpreted. Its meaning is as follows: We consider $E_{\rm{Pulse}}(j)\in\mathbb{R}\setminus\{0\}$ and $E_{\rm{Pulse}}(i): = 0$ for all $i\ne j$. Then, we solve the optimal control problem (2.9) in order to obtain $C^*[E_{\rm{Pulse}}(j)](\cdot)$. The partial derivative in (3.1) is then defined as

In the same way, partial derivatives with respect to other quantities occurring in the optimization can be defined.

Under appropriate differentiability assumptions on $E^*$, $C^*$ and the optimal value function $W^*$, Equation (3.1) can be rewritten as

used as the common definition in [4, Sec. 2.1], for example. Here, "[t]he numerator is the marginal impact of emissions at time $j$ on welfare, while the denominator is the marginal welfare value of a unit of aggregate consumption in period $j$" [4, Section 2.1].

Since the optimal value function $W^*$ and the optimal total emissions $E^*(j)$ and optimal consumptions $C^*(j)$ are only implicitly known as the optimal solutions of the optimal control problems (2.9), a direct calculation of the social cost of carbon through (3.3) is not possible. Nevertheless, there are several efficient ways to approximate (3.3).

3.1.   Direct discretization based on the difference quotient

Using (3.2), we can approximate the quantity in (3.3) at a fixed time $j\in\{0, \ldots, N-1\}$ by performing three experiments to compute approximations

To this end three optimal control problems can be solved.

1. A baseline solution of (2.9) that generates a welfare value $W^*$ corresponding to $E_{\rm{ Pulse}}(j) = 0$ and $C_{\rm{ Pulse}}(j) = 0$ for all $j = 0, \ldots, N-1$.

2. A second solution of (2.9) with a small non-zero perturbation on emissions at time $j\in \{0, \ldots, N-1\}$; i.e., $E_{\rm{ Pulse}}(j) \neq 0$. This yields a second (optimal) welfare value $W^*[E_{\rm{ Pulse}}(j)]$.

3. A third solution of (2.9) with a small non-zero perturbation on consumption at time $j\in\{0, \ldots, N-1\}$, i.e., $C_{\rm{ Pulse}}(j) \neq0$, that yields a third (optimal) welfare value $W^*[C_{\rm{ Pulse}}(j)]$.

Then, the SCC in (3.3) can be approximated by

Since (3.4) is a first order approximation of the derivative of $W^*$ with respect to $E^*(j)$ and $C^*(j)$, respectively, the right-hand side of (3.5) converges with order $O(C_{\rm{ Pulse}}(j)+E_{\rm{ Pulse}}(j))$ to (3.3) for $C_{\rm{ Pulse}}(j)\rightarrow 0$ and $E_{\rm{ Pulse}}(j) \rightarrow 0$ and under the assumption that the optimal control problem (2.9) can be solved with arbitrary precision. To calculate ${\rm S_{CC}}(j)$ for all $j\in \{0, \ldots, N-1\}$ using the approximation (3.5), $2N+1$ optimal control problems need to be solved.

Remark 1. The order of convergence follows from Taylor expansions of $W^*$ under the assumption that $W^*$ is twice continuously differentiable with respect to $E$ and $C$. In particular, the derivation of the error term $O(C_{\rm{ Pulse}}(j)+E_{\rm{ Pulse}}(j))$ relies on

and is given in detail in Appendix D.

3.2.   SCC calculation based on dual variables

Most modern NLP solvers, including IPOPT ^[26] and CONOPT ^[27], for example, return the optimal Lagrange multipliers, i.e., the dual variables of all explicitly enforced equality and inequality constraints in addition to the optimal solution $W^*$, $\mathit{\boldsymbol{s}}^*$, $\mathit{\boldsymbol{\mu}}^*$, are calculated. This can be used to compute the SCC directly via dual variables. All that is necessary to do is formulating (2.9) in a way such that $C(i)$ and $E(i)$ are decision variables of the NLP. Using (2.5) allows a reformulation of (CAR) in terms of

Instead of (2.9), we can thus consider the modified optimal control problem

where

are defined similar to (2.10). The optimal control problems (2.9) and (3.8) are identical in the sense that the optimal solutions $W^*$, $\mathit{\boldsymbol{s}}^*$, and $\mathit{\boldsymbol{\mu}}^*$ coincide. However, in (3.8) we avoid the substitution of (2.6) and (2.5).

The main difference between (2.9) and (3.8) is that in the latter case, consumption and emissions are treated as pseudo decision variables defined via explicit equality constraints, while in the former case one substitutes $C(i)$ and $E(i)$ in the objective and the dynamics by the expressions given in (2.6) and (2.5). The implementation of (3.8) in an NLP solver such as IPOPT yields values for the optimal Lagrange multipliers of all constraints.

To make this more precise, note that we can rewrite (3.8) in the form of the parameter dependent optimization problem

in the optimization variables $\mathit{\boldsymbol{x}} = \{\mathit{\boldsymbol{s}}, \mathit{\boldsymbol{\mu}}, \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{E}}\}$. Here,

denotes the objective function in the variables $\mathit{\boldsymbol{s}}, \mathit{\boldsymbol{\mu}}, \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{E}}$, and

capture (2.5) and (2.6) for $i = 0, \ldots, N-1$. The pulses $E_{\rm{ Pulse}}$ and $C_{\rm{ Pulse}}$ are included through the parameter ${{p\mathit{\boldsymbol{}}}} = [\mathit{\boldsymbol{p}}_E^T \ \ \mathit{\boldsymbol{p}}_C^T]^T$,

in the right-hand side, and the optimal value $f^*(\cdot)$ depends on the selection of $\mathit{\boldsymbol{p}}$. The equality and inequality constraints $\mathit{\boldsymbol{x}} \mapsto k(\mathit{\boldsymbol{x}})$ and $\mathit{\boldsymbol{x}} \mapsto g(\mathit{\boldsymbol{x}})$ cover the remaining constraints of (3.8), i.e., $X(0) = X_0$, (3.7), (CLI)–(EI), (2.2)–(2.4), $s(i), \mu(i) \in [0, 1]$ for $i\in \{0, \ldots, N-1\}$ are encoded through appropriate equality and inequality constraints in the decision variable ${{x\mathit{\boldsymbol{}}}}$.

The Lagrangian corresponding to (3.9) with Lagrange multipliers $(\mathit{\boldsymbol{\lambda}}_E, \mathit{\boldsymbol{\lambda}}_C, \mathit{\boldsymbol{\delta}}, \mathit{\boldsymbol{\kappa}})$ is defined as

If $\mathit{\boldsymbol{x}}^*$ and $(\mathit{\boldsymbol{\lambda}}_E^*, \mathit{\boldsymbol{\lambda}}_C^*, \mathit{\boldsymbol{\delta}}^*, \mathit{\boldsymbol{\kappa}}^*)$ are an optimal solution-Lagrange multiplier pair (see [28, Sec. 5.1.3], for example), then the optimal value function satisfies

Moreover, under appropriate conditions (second order sufficiency conditions [28, Prop. 3.3.2], for example), the partial derivatives of the optimal value function in a neighborhood around $\mathit{\boldsymbol{p}} = 0$ satisfy

for all $i\in\{0, \ldots, N-1\}$, according to ^[28].

While in the optimal value function, necessarily, the terms containing Lagrange multipliers vanish, and only the objective function is important according to (3.11), in $\nabla_{\mathit{\boldsymbol{p}}} f^*(\mathit{\boldsymbol{p}})$ in (3.12), it turns out that the dependency of $(\mathit{\boldsymbol{x}}^*, \mathit{\boldsymbol{\lambda}}_E^*, \mathit{\boldsymbol{\lambda}}_C^*, \mathit{\boldsymbol{\delta}}^*, \mathit{\boldsymbol{\kappa}}^*)$ on $\mathit{\boldsymbol{p}}$ can be ignored when computing the gradient through the Lagrangian (3.10). Thus, the Lagrange multipliers $\mathit{\boldsymbol{\lambda}}_E$, $\mathit{\boldsymbol{\lambda}}_C$ are a measure of the sensitivity of the objective function with respect to changes in the right-hand side of the equality constraints. In economics, this sensitivity, i.e., the change in the objective function under a marginal change of a constraint, is also referred to as the shadow price (see ^[29,30], for example) and it provides an interpretation of how valuable a specific commodity is.

These observations finally allow us to return to (3.3). If we recall the connection between (3.8) and (3.9) we observe that

is satisfied for all $j\in\{0, \ldots, N-1\}$. Thus, according to (3.3), the SCC is given by

Using (3.14), $\mathrm{S_{CC}}(j)$ can be computed for all $j\in\{0, \ldots, N-1\}$ by solving only a single optimization problem (3.8). In addition, unlike (3.5), Equation (3.14) is not an approximation of (3.3) but an exact expression. However, the pulse experiments might be more intuitive and from an interpretation point of view easier to understand than the connection to the Lagrange multipliers.

3.3.   SCC calculation based on the discounted damages stream

In this section we present a third way to compute the SCC. This is based on the net present value, which, as described in the introduction, has been used in much of the climate economics research but not in a consistent manner. Here, we will show that a time-varying Ramsey discount factor, accounting for a time-varying growth rate of consumption, is required for consistency between calculating the SCC as a discounted damages stream and calculating the SCC as a ratio of marginal prices. The result relies on the assumption that (2.9) satisfies second order sufficiency conditions. The precise definition and a corresponding result are reported in Theorem 2 in Appendix A.

Theorem 1. Consider the optimal control problem (2.9) and assume that (2.9) satisfies second order sufficiency conditions. Then, the SCC at time $j\in \{0, \ldots, N-1\}$ is given by

where the Ramsey discount factor between periods $j$ and $i\in \{0, \ldots, N-1\}$ is defined as

and the growth rate of per capita consumption up to period $i$ is defined as

A proof of Theorem 1 is given in Appendix C. Based on Theorem 1, the SCC can again be computed using a finite difference approximation. While $\delta_{j, i}$ in (3.16) can be computed explicitly, the partial derivative $\frac{\partial C^*(i)}{\partial E^*(j)}$ in (3.15) needs to be approximated. We can thus perform two experiments to approximate the SCC at time $j\in \{0, \ldots, N-1\}$ using Equation (3.15):

1. A baseline solution of (2.9) that generates a welfare value $W^*$ corresponding to $E_{\rm{ Pulse}}(i) = 0$ and $C_{\rm{ Pulse}}(i) = 0$ for all $i = 0, \ldots, N-1$ and additionally returns the nominal consumption pathway $C^*(i)$ for $i = 0, \ldots, N-1$.

2. A second solution of (2.9) with a small non-zero perturbation on emissions at time $j\in\{0, \ldots, N-1\}$, i.e., $E_{\rm{ Pulse}}(j) \neq 0$. This yields a second (optimal) welfare value $W^*[E_{\rm{Pulse}}(j)]$ and additionally a second (optimal) consumption pathway $C^*[E_{\rm{Pulse}}(j)](i)$ for $i = 0, \ldots, N-1$.

With the resulting pathways, we can calculate the discount rate $\delta_{j, i}$ and can approximate the partial derivatives, similar to (3.4), by $(C^*[E_{\rm{Pulse}}(j)](i) - C^*(i))/(E_{\rm{Pulse}}(j))$. The social cost of carbon can then be approximated by

Remark 2. Using the Taylor expansion

i.e.,

it follows that (3.18) converges with order $O(E_{\rm{ Pulse}}(j))$ to $\mathrm{S_{CC}}(j)$, similar to the first order convergence of (3.5).

Similar to the method discussed in Section 3.1, the SCC calculation discussed in this section relies on pulse experiments and on first order approximations. In particular, while the method in Section 3.1 relies on a first order approximation of $\frac{\partial W^*}{\partial E^*(j)}$ and $\frac{\partial W^*}{\partial C^*(j)}$, the third method discussed here uses a first order approximation of $\frac{\partial C^*}{\partial E^*(j)}$, respectively. However, as discussed in Remarks 1 and 2, the order of convergence of both methods is the same. The value of the formula in Theorem 1 is that the expression (3.15) (or (3.18), respectively, ) has a familiar interpretation as the net present value (NPV) of the stream of damages due to an emissions pulse at time $j$ and is conventionally expressed in units of ＄/$\mathrm{tCO}_2$. Note, however, that in computing the NPV in (3.18), the discount factor $\delta_{j, i}$ is time-varying due to its dependence on the growth rate $g(\cdot)$ of per capita consumption; cf. (3.16). This marks a subtle, but potentially interesting point of departure between (3.18) and a large body of literature which applies the well-known, continuous-time Ramsey equation for the discount rate $r = \rho + \alpha \cdot g$ for constant $g$, to the discrete-time setting (see ^[31], for example).

4.   Numerical results

In this section we discuss three aspects in the calculation of the SCC using the DICE model. For the calculations we use the standard setting of the DICE-2016R model (see [24, Table 2] for the parameter selection) if not specified differently. For the numerical calculations we use the implementation ^[32]. While we use a horizon of $N = 100$, consistent with the parameters in DICE-2016R, we only focus on the first time steps in the visualizations.

Before we focus on the calculation of the SCC, we recall one of its illustrative definitions and visualizations used in the literature and highlight the anticipatory response of solutions of the optimal control problem (2.9) when emission or consumption impulses are added at different time steps. In the study by Nordhaus W ^[4] the SCC is illustrated as "the economic cost caused by an additional ton of carbon dioxide emissions" which is usually illustrated by adding a pulse to the first time step in the future [4, Figure 1]. In particular, with respect to [4, Figure 1], if time step 1 indicates the initial conditions, an additional ton of carbon dioxide emissions is added at time step 2. Under these assumptions, baseline emission and consumption pathways are compared with perturbed emission and consumption pathways. The baseline emissions and consumption pathways using the DICE-2016R model are visualized in Figure 1.

By adding the additional tonne at the earliest possible time in the future, one might get the misleading impression that a perturbation $E_{ \rm{Pulse}}(j)$, $j\in \{0, \ldots, N-1\}$, only affects emissions and consumption at time steps $k\in \{j, \ldots, N-1\}$. This is however not the case, as shown in Figure 2, where a pulse is added at time $j = 1$, $j = 10$, and $j = 20$, respectively. In Figure 2 the differences between the baseline solution, i.e., $E_{ \rm{Pulse}}(j) = C_{ \rm{Pulse}}(j) = 0$ for all $j\in \{0, \ldots, N-1\}$, and perturbed solutions with $C_{ \rm{Pulse}}(j) = 1$ for $j = 1$, $j = 10$, and $j = 20$ are shown. A pulse at time $j\in \{0, \ldots, N-1\}$ has an impact on the whole pathways $C[E_{ \rm{Pulse}}(j)](\cdot)$ and $E[E_{ \rm{Pulse}}(j)](\cdot)$ due to the optimization involved in their computation. This also becomes clear from the definition of the partial derivative (3.2) used in the formula for the SCC: For computing this derivative, first the pulse is added, and then the optimization is performed, implying that the optimizer is aware of the pulse at time $j$ when performing the optimization for the times $i = 0, \ldots, j-1$. In other words, if there is a change of emissions in the future, which importantly by assumption is known, then the change is anticipated and encoded in the optimal control problem (2.9). This explains the non-negligible differences in $C[C_{ \rm{Pulse}}(20)](i)$ and $C(i)$ and in $E[C_{ \rm{Pulse}}(i)](20)$ and $E(i)$, respectively, for $i < 20$, for example.

As a next step, we show that all three methods discussed in Section 3 essentially lead to the same SCC value. Figure 3 compares the three methods with each other, where for the pulse experiments the values $C_{ \rm{Pulse}}(j) = 1$ and $C_{ \rm{Pulse}}(j) = 1$, $j\in \{0, \ldots, N-1\}$, were used for the finite difference approximations (which corresponds to a one tonne impulse in the implementation). On the scale on the left of Figure 3, the SCC values are indistinguishable. However, the difference between the Lagrange multiplier method and the finite difference approximation (3.5) and the formula (3.18) derived in Theorem 1 appears to grow nonlinearly with $j\in \{0, \ldots, N-1\}$ as visualized in Figure 3 on the right.

In Figure 4 the errors in the SCC calculation with respect to $E_{ \rm{Pulse}}$ and $C_{ \rm{Pulse}}$ are shown. In particular, the impulses $E_{ \rm{Pulse}}(j)$ and $C_{ \rm{Pulse}}(j)$ are varied from 0.01 to 5 for different values of $j\in \{0, \ldots, N-1\}$. As expected from Remarks 1 and 2, the error grows linearly with $E_{ \rm{Pulse}}$ and $C_{ \rm{Pulse}}$, respectively.

Remark 3. While in theory the SCCs computed with the three methods coincide, the precise numerical values we obtain are sensitive with respect to $E_{ \rm{Pulse}}$ and $C_{ \rm{Pulse}}$, the chosen optimization algorithms and the parameter selection of the optimization algorithm. Figure 5 shows the same results as Figure 4 with the only difference that the horizon $N = 100$ has been replaced with $N = 150$. In this case the parameter selection used in the IPOPT solver does not lead to reliable results for all impulse perturbations.

As a final visualization we show the change of the SCC with respect to the discount rate in Figure 6, i.e., $\frac{d}{d \rho}\mathrm{S_{CC}}(j)$ is calculated and shown in Figure 6. Here, $\frac{d}{d \rho}\mathrm{S_{CC}}(j)$ is approximated through a finite difference approximation

where $\mathrm{S_{CC}}[\rho + \Delta \rho](j)$ and $\mathrm{S_{CC}}[\rho](j)$ are calculated using the Lagrange multiplier method presented in Section 3.2.

5.   Conclusions

In this paper, we have investigated and compared three different methods for the computation of the SCC in integrated assessment models. While our discussion is focused on the DICE model, it is equally relevant for other models including PAGE and FUND, to obtain a characterization of the SCC that is consistent among the literature. While the first two methods, which rely on perturbations of a baseline model and on dual variables of an appropriately defined optimal control problem, respectively, are commonly known and used, the third method relies on a novel formula stated in Theorem 1. The novel formula allows a direct analysis of the impact of the discount factor in the calculation of the NPV and the SCC through Eq (3.16), which subsequently provides insights in SCC values stemming from different integrated assessment models. The consistency of the three methods is analyzed analytically through first order approximations of Taylor polynomials and through dual variables and shadow price interpretations. The theoretical findings are additionally confirmed through numerical simulations, while it is also highlighted that an appropriate parameter selection in the optimal control formulation and in the underlying solver is important. Future work will focus on the reliability of the three methods to compute the SCC with respect to the horizon $N$ and the parameters of the DICE model as discussed in Remark 3.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare no conflict of interest.

Appendix

A.   Optimality conditions for constrained and parametric optimization

In this section we recapitulate optimality conditions for constrained optimization problems from Bertsekas DP' study ^[28], which are needed for the proof of the main Theorem 1 in Appendix C. In particular, for the convenience of the reader, we rewrite results for minimization problems in terms of maximization problems discussed in this paper. We consider optimization problems of the form

where $x\in \mathbb{R}^n$, $f:\mathbb{R}^n \rightarrow \mathbb{R}$, $h:\mathbb{R}^n \rightarrow \mathbb{R}^m$, and $g: \mathbb{R}^n \rightarrow \mathbb{R}^p$ for $n, m, p \in \mathbb{N}$. With the Lagrange multipliers $\lambda \in \mathbb{R}^m$ and $\mu \in \mathbb{R}^p_{\geq 0}$, the Lagrangian corresponding to (A.1) is defined as

Additionally, for $x\in \mathbb{R}^n$ we define

as the corresponding active set with respect to the inequality constraints $g$.

Theorem 2 (Second order sufficiency conditions, [28, Prop. 3.3.2]). Consider the optimization problem (A.1) and assume that $f$, $h$, and $g$ are twice continuously differentiable. Let $x^* \in \mathbb{R}^n$, $\lambda^* \in \mathbb{R}^m$, and $\mu^* \in \mathbb{R}^m$ satisfy

for all $y\neq 0$ such that

Finally, assume that

Then, $x^*$ is a strict local maximum of (A.1) and we say that (A.1) satisfies the second order sufficient condition at $x^*$.

Additionally we will need the following result for parameter dependent optimization problems where the parameter only appears in the objective function $f:\mathbb{R}^n \times \mathbb{R} \rightarrow \mathbb{R}$.

Theorem 3 ([33, Corollary 3.4.2]). Consider the parameter dependent optimization problem

with Lagrangian

Lagrange multipliers $\lambda$, $\mu$, and parameter $p\in \mathbb{R}$ in the objective function. Assume that the second order sufficiency conditions of Theorem 2 are satisfied. Then, in a neighborhood around $p = 0$, it holds that

B.   Proof of Theorem 1

To prove Theorem 1 we will make use of results in ^[33] and the results summarized in Appendix A. Before we begin with the actual proof, we note that for

the optimization problems (2.9) can be written in the general form

where the equality constraints are eliminated, and the box constraints are captured in the function $g$. Here, the parameter $p$ represents the pulses $E_{ \rm{Pulse}}(j)$ and $C_{ \rm{Pulse}}(j)$, or equivalently, $p$ represents a change in $C(j)$ and $E(j)$, for $j\in \{0, \ldots, N-1\}$. Since we only consider the perturbation with respect to $C$ or $E$ and only with respect to one $j\in \{0, \ldots, N-1\}$ at a time, we assume that $p$ is one dimensional. With these definitions we can use Theorem 3 (without the function $h$) in the following. Additionally we note that according to the definition of the functions involved in the optimal control problem (2.9), the functions $f$ and $g$ are twice continuously differentiable.

Proof of Theorem 1: We begin the proof by rewriting the partial derivatives

for $j\in \{0, \ldots, N-1\}$. To this end we define

which is a scaling of $C(i)$ through the known parameter $L(i)$, $i\in \{0, \ldots, N-1\}$. Thus, the utility $U$ defined in (2.7) can be written as a function of $c_i$,

For the derivative we obtain

and with respect to the original variables, it holds that

With Theorem 3, and in particular with (A.4) for $p = E_{\rm{Pulse}}(j)$, $j\in \{0, \ldots, N-1\}$, it holds that

Using the notation

to shorten the expressions and the equation derived in Equation (B.5), the right-hand side can be written in the form

Similarly, the application of Equation (A.3) in Theorem 3, where $p$ represents a change in $C(j)$, $j\in \{0, \ldots, N-1\}$, leads to the expression

As a next step, we derive Taylor expansions of the optimal value function $W^*$. In particular, for parameters $p_i$ representing a change in $C^*(j)$, (i.e., a pulse $C_{ \rm{Pulse}}(i)$) for $i\in\{0, \ldots, N-1\}$, we derive Taylor expansions in a neighborhood of $p_i = 0$. Using the definition of $c_i$ in Equation (B.2), the consumption (2.6) can be equivalently written as

Moreover, with the definition $c_{\rm{ Pulse}}(i) = \tfrac{10^3}{L(i)}C_{\rm{ Pulse}}(i)$ and $c_i = \tilde{c}_i + c_{\rm{ Pulse}}(i)$, Equation (B.8) can be written as

and the utility function can be written as

Thus, instead of the parameter $C_{\rm{Pulse}}(j)$ in the consumption (2.6), we can use a scaled parameter directly in the utility function (2.7), and the optimal value function satisfies

with $p_i = c_{\rm{ Pulse}}(i)$. According to Theorem 3 (using Equation (A.5)), the optimal value function is twice continuously differentiable with respect to the parameter $p_i$, and using Theorem 3, Equation (A.4), the Taylor approximation is given by

where

for $\xi \in [0, p_i]$. For $i, j\in \{0, \ldots, N-1\}$, we take specific perturbations $p_i\neq 0$, $p_j \neq 0$ with the property $W^*[p_i] = W^*[p_j]$ so that the Taylor approximation satisfies

A $p_j$ satisfying (B.12) can be found for each sufficiently small $p_i$ because the terms in front of $p_i$ and $p_j$ in (B.12) are not zero due to the assumption that (2.9) satisfies the second order sufficiency conditions. In fact, $\gamma_j L(j) (\tilde{c}^*_j[0])^{-\alpha} = 0$ implies that $s^*(i) = 1$ or $\theta_1(i) \mu^*(i)^{\theta_2} = 1$ (see Equation (B.9)). If $s^*(i) = 1$, then $\mu^*(i)\in [0, 1]$ can be selected arbitrarily, and if $\theta_1(i) \mu^*(i)^{\theta_2} = 1$, then $s^*(j) \in [0, 1]$ can be selected arbitrarily, i.e., the optimal solution of (2.9) is not characterized through a strict local maximum, which contradicts Theorem 2. Under these assumptions stated above, $p_j$ is a function of $p_i$, and $\lim_{p_i\rightarrow 0} p_j(p_i) = 0$ (since $W^*$ is continuously differentiable). Rearranging terms yields

and allows us to investigate the limit $p_i \rightarrow 0$, $p_i\neq 0$, (i.e., $p_j(p_i)\rightarrow 0$). We observe that

and since the remainder in Equation (B.11) satisfies $R_i(p_i) = O(p_i^2)$, we obtain

To be able to neglect the second term in the right-hand side of Equation (B.13), it is left to show that $\frac{1}{p_i}R_j(p_j(p_i))$ converges to zero for $p_i \rightarrow 0$. With the definition of the remainder $R_j(p_j) = \frac{1}{2} \frac{d^2}{d p_j^2} W^*_j[\xi] p_j(p_i)^2$ in Equation (B.11), the limit satisfies

Consequently, using Equation (B.4) together with $\tilde{c}_i = c_i$ in the case $p_i = 0$, it holds that

To return to the original variables $C(i)$, we recall the coordinate transformation (B.2), and use the notation $P_i = C_{\rm{Pulse}}(j) = \frac{L(i)}{10^3} p_i.$

With the above computations we are finally in the position to derive the partial derivative of $C^*(i)$ with respect to $C^*(j)$ for $i, j\in \{0, \ldots, N-1\}$. In particular, using the chain rule and the definition of the derivative, it holds that

Using the relation between $p_i$ and $P_i$, this expression can be equivalently written as

Finally, with the result derived in Equation (B.14), it holds that

where the second term is equal to one since $W^*[p_i] = W^*[p_j(p_i)]$ according the assumption made at the beginning of the derivation of (B.14) in (B.12).

To complete the proof we combine the results so far in the definition of $\mathrm{S_{CC}}(j)$. With Equations (B.6) and (B.7), the SCC can be expressed as

and which becomes

through the identity (B.15). Finally, reformulating individual terms and recalling the definition of $\delta_{ji}$ in Theorem 1 completes the proof:

C.   Order of convergence

As a last result in this paper, we show that (3.5) converges with order $O(C_{\rm{ Pulse}}(j)+E_{\rm{ Pulse}}(j))$ to ${\rm S_{CC}}(j)$. Using the notation introduced in (3.4) and (3.5), the SCC satisfies

for unknown functions $\Phi_E$, $\Phi_C$, and $\Phi$. However, according to the Taylor expansions in (3.6), it holds that $\Phi_E(E_{\rm{ Pulse}}(j)) = O(E_{\rm{ Pulse}}(j))$ and $\Phi_C(C_{\rm{ Pulse}}(j)) = O(C_{\rm{ Pulse}}(j))$.

We investigate properties of the function $\Phi$. To simplify the notation, we drop the index in the derivations, i.e., we investigate the equation

or equivalently

We further simplify the notation and define $\Delta W_E^* = W^*[E] - W^*$ and $\Delta W_C^* = W^*[C] - W^*$. Then, multiplying by the denominator leads to the equation

After canceling and rearranging terms, we are left with the expression:

Since $W^*$ is differentiable, i.e., the limits

exist, it holds that $\Delta W_E^* = O(E_{\rm{ Pulse}}(j))$ and $\Delta W_C^* = O(C_{\rm{ Pulse}}(j))$. Using these properties in the last equation simplifies the expression to

or equivalently

In particular, $\Phi$ satisfies $\Phi(E, C) = \mathcal{O}(C+E)$, proving the assertion.

D.   Parameters of the DICE model

In this section we report the parameters used to define the R-DICE2016 model. The parameters in Table 1 are taken from [24, Appendix A], which additionally contains the parameters for R-DICE2013. Moreover, the initial conditions used to define $X_0$ in (3.9) are defined as