Behind the scenes solving ODEs of a model of optimal climate policy
by Stephen DeCanio, Alan Sanstad, Charles Manski
by Stephen DeCanio, Alan Sanstad, Charles Manski
Our paper (Charles F. Manski, Alan H. Sanstad, and Stephen J. DeCanio, “Addressing partial identification in climate modeling and policy analysis”, Proceedings of the National Academy of Sciences, Vol 118, No. 15 (2021), https://doi.org/10.1073/pnas.2022886118), examines the problem of choosing carbon emissions abatement policies to minimize the total cost of abatement and climate damages, in a simple climate/economic integrated assessment model. Because there are multiple, co-existing computational models of the physical climate system that are considered equally plausible by climate scientists, we used the Minimax Regret technique from decision theory to show how a near-optimal abatement policy can be chosen in the face of this model uncertainty without applying poorly-justified weighting schemes that produce model averages.
In the course of employing Mathematica for both symbolic and numerical computations described in the paper, we encountered several issues, specifically involving the solution of a system of ordinary differential equations, that may be of interest to other Mathematica users. These involve (a) the fact that in this case human intervention was required to get DSolve to return all the solutions, and (b) the algebraic expressions generated by Mathematica’s simplification algorithm(s) can be difficult to parse. Numerical experimentation using Mathematica assisted us in identifying and then resolving these issues.
The optimization problem and its solution are described on pages 6-7 of the main text of the paper and pages 3-4 of the online Supporting Information that accompanies the paper. The total cost function was the sum of abatement costs and climate damages, with these costs specified as quadratic functions and parameterized with values drawn from the literature.
The first-order conditions for minimization of this cost function give two relatively simple differential equations, one for abatement A[t] and the other for the cumulative buildup of greenhouse gases S[t]. We used Mathematica’s symbolic manipulation capabilities to check “by hand” solutions to the differential equations (these are shown on pages 5-7 of the Supporting Information), and to carry out the numerical implementation of the solutions. The issues discussed in this blog post are covered in Sections 3.1-3.4 of the Supporting Information.
The model gives rise to a system of two relatively simple, first-order linear ordinary differential equations with constant coefficients:
In the course of employing Mathematica for both symbolic and numerical computations described in the paper, we encountered several issues, specifically involving the solution of a system of ordinary differential equations, that may be of interest to other Mathematica users. These involve (a) the fact that in this case human intervention was required to get DSolve to return all the solutions, and (b) the algebraic expressions generated by Mathematica’s simplification algorithm(s) can be difficult to parse. Numerical experimentation using Mathematica assisted us in identifying and then resolving these issues.
The optimization problem and its solution are described on pages 6-7 of the main text of the paper and pages 3-4 of the online Supporting Information that accompanies the paper. The total cost function was the sum of abatement costs and climate damages, with these costs specified as quadratic functions and parameterized with values drawn from the literature.
The first-order conditions for minimization of this cost function give two relatively simple differential equations, one for abatement A[t] and the other for the cumulative buildup of greenhouse gases S[t]. We used Mathematica’s symbolic manipulation capabilities to check “by hand” solutions to the differential equations (these are shown on pages 5-7 of the Supporting Information), and to carry out the numerical implementation of the solutions. The issues discussed in this blog post are covered in Sections 3.1-3.4 of the Supporting Information.
The model gives rise to a system of two relatively simple, first-order linear ordinary differential equations with constant coefficients:
In[]:=
A'[t]=pA-(β)αS[t];
2
m
and
In[]:=
S'[t]=q[t]-A[t];
where p is the discount rate, β and m are parameters of the climate damage function, α is a parameter of the abatement cost function, and q[t] is a forcing function that represents an exogenous baseline emissions scenario (absent any policy) that is commonly used in the literature for climate simulations. These two equations are given as equation SI.3.3 of the Supporting Information, with the notational change that the published paper uses “E” for cumulative greenhouse gas accumulations instead of “S,” because Mathematica reserves “E” for the base of the natural logarithms. In addition the published paper uses “B” instead of “q” to denote the forcing function and δ rather than p for the discount rate.
The parameters θ and φ are found by non-linear fitting (NonlinearModelFit) of an appropriately-shaped function to the exogenous no-policy emissions scenario (see equation SI.2.1 and Figure S1 of the Supporting Information.), specifically
where p is the discount rate, β and m are parameters of the climate damage function, α is a parameter of the abatement cost function, and q[t] is a forcing function that represents an exogenous baseline emissions scenario (absent any policy) that is commonly used in the literature for climate simulations. These two equations are given as equation SI.3.3 of the Supporting Information, with the notational change that the published paper uses “E” for cumulative greenhouse gas accumulations instead of “S,” because Mathematica reserves “E” for the base of the natural logarithms. In addition the published paper uses “B” instead of “q” to denote the forcing function and δ rather than p for the discount rate.
The parameters θ and φ are found by non-linear fitting (NonlinearModelFit) of an appropriately-shaped function to the exogenous no-policy emissions scenario (see equation SI.2.1 and Figure S1 of the Supporting Information.), specifically
In[]:=
q[t]=(θ*t)+q0*;
θφ
-θ(t-φ)
The matrix of this system has two eigenvalues, λ1 and λ2, given by
In[]:=
λ1=Eigensystemp,-(β)α,{-1,0}[[1]][[1]]
2
m
Out[]=
pα-α+4β
α
2
p
2
m
2α
In[]:=
λ2=Eigensystemp,-(β)α,{-1,0}[[1]][[2]]
2
m
Out[]=
pα+α+4β
α
2
p
2
m
2α
It is evident that λ1 is negative and λ2 is positive, and that λ1 + λ2 = p.
In[]:=
Clear[λ1,λ2]
After uncoupling the two equations using a standard diagonalization procedure, applying DSolve, reversing the diagonalization to recover the original variables, and finding constants of integration consistent with the initial conditions, we found the following solutions for A[t] and S[t] expressed in terms of the eigenvalues and the parameters of the forcing function:
In[]:=
A=-λ1S0(λ1-λ2)-θ+q0(θ+λ1)λ2+q0(θ+λ1)+θ(1+t(θ+λ1))λ2+λ1θ+q0(θ+λ2)-λ2q0(θ+λ2)+θ(1+t(θ+λ2))((λ1-λ2));
tλ1
2
(θ+λ1)
2
(θ+λ2)
θφ
2
(θ+λ2)
-t(θ+λ1)
θφ
2
(θ+λ2)
2
(θ+λ1)
θφ
-tθ
2
(θ+λ1)
θφ
2
(θ+λ1)
2
(θ+λ2)
and
In[]:=
S=S0(λ1-λ2)-θ+q0(θ+λ1)λ2+q0(θ+λ1)+θ(1+t(θ+λ1))λ2+λ1θ+q0(θ+λ2)-λ1q0(θ+λ2)+θ(1+t(θ+λ2))((λ1-λ2));
tλ1
2
(θ+λ1)
2
(θ+λ2)
θφ
2
(θ+λ2)
-t(θ+λ1)
θφ
2
(θ+λ2)
2
(θ+λ1)
θφ
-tθ
2
(θ+λ1)
θφ
2
(θ+λ1)
2
(θ+λ2)
Substituting the the basic parameters back into the expressions for the eigenvalues, we obtained the following expressions as the general solutions for A[t] and S[t]:
In[]:=
Ageneral=pα-α+4β-pα+4β2tθ+q0+tθα+4β+2tp-α+4β-16S0α+4β-pα+α+4β-q0α+4β+2αθ+q0α(p+2θ)+pα-α+4βq0α+4β+2αθ+q0α(p+2θ)+t-p+α+4β-2θpα+α+4βq0pα-α+4β+2αθ+θ-tα+4β+α(2+pt+2tθ)32α+4ββ-θ(p+θ);
α
2
p
2
m
-tθ
α
+2
p
2
m
2
p-αα+4β+2θ
3/2
α
2
p
2
m
3/2
α
θφ
α
θ+pα
q0+θφ
θφ
2
p
2
m
α
θ+1
2
2
p
2
m
α
3/2
α
2
p
2
m
2
(β-αθ(p+θ))
2
m
α
2
p
2
m
2
pα+α+4β+2αθ
α
2
p
2
m
α
2
p
2
m
θφ
α
2
p
2
m
2
pα-α+4β+2αθ
α
2
p
2
m
α
2
p
2
m
θφ
1
2
2
p
2
m
α
α
2
p
2
m
2
pα+α+4β+2αθ
α
2
p
2
m
α
2
p
2
m
θφ
α
2
p
2
m
9/2
α
2
p
2
m
2
2
m
α
and
In[]:=
Sgeneral=--pα+4β2tθ+q0+tθα+4β+2tp-α+4β-16S0α+4β-pα+α+4β-q0α+4β+2αθ+q0α(p+2θ)+pα-α+4βq0α+4β+2αθ+q0α(p+2θ)+t-p+α+4β-2θpα+α+4βq0pα-α+4β+2αθ+θ-tα+4β+α(2+pt+2tθ)16α+4ββ-θ(p+θ);
-tθ
α
+2
p
2
m
2
αα+4β-(p+2θ)
2
p
2
m
3/2
α
θφ
α
θ+pα
q0+θφ
θφ
2
p
2
m
α
θ+1
2
2
p
2
m
α
3/2
α
2
p
2
m
2
(β-αθ(p+θ))
2
m
α
2
p
2
m
2
pα+α+4β+2αθ
α
2
p
2
m
α
2
p
2
m
θφ
α
2
p
2
m
2
pα-α+4β+2αθ
α
2
p
2
m
α
2
p
2
m
θφ
1
2
2
p
2
m
α
α
2
p
2
m
2
pα+α+4β+2αθ
α
2
p
2
m
α
2
p
2
m
θφ
α
2
p
2
m
7/2
α
2
p
2
m
2
2
m
α
The next step in our paper was to calculate the cost of the optimal abatement path. This cost, as described on pages 8-9 of the paper, is given by:
In[]:=
Integral1=(1/2)Integrate[*((α*)+(β**)),{t,0,∞},Assumptions{θ>0,φ>0,p>0,α>0,β>0,m>0}]
-pt
2
Ageneral
2
m
2
Sgeneral
Out[]=
This is obviously quite a complicated expression (over 1000 terms). Can it be simplified?
In[]:=
FullSimplify[%,Assumptions{θ>0,φ>0,p>0,α>0,β>0,m>0}]
Out[]=
Even after simplification, this is an extremely complicated expression. Economists like to get a sense of what their solutions actually mean, but there is no way to tell what is going on with this expression just from looking at it. So, we conducted some numerical experiments to explore the shape of the cost function for different combinations of parameters. For example, here is a 3-D plot of the total cost function varying the parameters α and β while holding the other parameters constant.
In[]:=
shapeIntegral1=Table[Integral1/.{p0.05,m0.002286328,θ0.012312473880977292,φ339.5651624412206,S0470,q08},{α,0.000078,0.0002,0.000001},{β,0.014,0.022,0.0001}];
In[]:=
αindex0=Table[i,{i,0.000078,0.0002,0.000001}];
In[]:=
βindex0=Table[i,{i,0.014,0.022,0.0001}];
In[]:=
ListPlot3D[{αindex0,βindex0,shapeIntegral1},AxesLabel{"βindex","αindex"},LabelStyleDirective[Blue,Large]]
Out[]=
It can be seen that the shape is smooth and regular, except that there appear to be some anomalies cutting through the middle of the plot. Examination of the denominator term of the simplified expression for Integral1 reveals a term that, for some possible parameter combinations, can be zero. (Everything else in the denominator is positive and non-zero.) A little algebra shows that β-αθ(p+θ) = 0 iff θ + λ1 = 0. Looking back at the solutions for A and S expressed in terms of the eigenvalues, we see that θ + λ1 appears in the denominators of those expressions as well. Hence, for this particular subspace of the parameter space, the solutions obtained by DSolve break down. It should be noted that this same issue arises if the differential equations are solved by hand using the textbook method of applying an integrating factor. This can be seen by following the derivation on pages 5-6 of the Supporting Information published online with the main text of our paper. The anomalies are the result of numerical instability arising when the parameters are close to the β-αθ(p+θ) = 0 locus. The fact that this instability produces some anomalies that are negative --- which is theoretically impossible in this model --- suggests that the numerators of Ageneral and Sgeneral may also be zero for these parameter values.. L’Hospital’s Rule can evaluate indeterminate expressions where the numerator and denominator both approach zero, but again, the general expressions for A and S cannot be solutions to the system in such cases.Does this mean that the system of ODEs has no solution if the parameters are in the region for which β-αθ(p+θ)= 0? No, because the original system can be rewritten imposing the condition β-αθ(p+θ) = 0 from the outset. This gives rise to a “special case” solution to the system. It can be solved using the same procedure of uncoupling, applying DSolve, calculating the constants of integration from the initial conditions, and recovering the original variables. The special case solutions are given by:
2
β-αθ(p+θ)
2
m
2
m
2
m
2
m
2
m
In[]:=
Aspecial=θ(p+θ)q0(p+2θ)+θ(1+pt+2tθ)-θ-S0--;
t(p+θ)-t(p+2θ)
θφ
3
(p+2θ)
-tθ
t(p+θ)2q0+tθ
θφ
2(p+2θ)
θθ+q0(p+2θ)
θφ
3
(p+2θ)
The numerical values we reported in the paper were all calculated using the general solutions. Section SI 3.4 of the paper describes the tests we carried out to show that these numerical values correspond to parameter values that are not in the subspace of the special solution.