Modules:

The module src includes functions used for solving ODEs and simulating trajectories of the project.

model.py

Functions that solve the HJBs in the draft paper.

src.model.minimize_g(y_grid, numy_bar, ems_star, φ_list, args, ε=3, tol=1e-06, max_iter=3000)

compute jump model with ambiguity over climate models

src.model.minimize_π(y_grid, numy_bar, ems_star, φ_list, args, with_damage=False, ε=2, tol=1e-07, max_iter=3000)

compute jump model with ambiguity over climate models

src.model.solve_baseline(y_grid, num_stop, ems_star, φ_list, args, ε=2, tol=1e-08, max_iter=3000)

compute jump model with ambiguity over climate models

src.model.solve_hjb_y(y, model_args=(), v0=None, ε=1.0, tol=1e-08, max_iter=10000, print_iteration=True)

Solve the HJB that is only related to y. y is the accumulative change of temperature.

Parameters

• y ((N,) ndarray) – An evenly spaced grid of y.

• model_args (tuple of model inputs:) –

  The first 9 elements are floats for \(η, δ, σ_y, \bar{y}, γ_1, γ_2, γ_3, ξ_w, ξ_a\).

  The last two are (L,) ndarrays for \(\{θ_\ell\}_{\ell=1}^L, \{\pi^a_\ell\}_{\ell=1}^L\).

• v0 ((N,) ndarray, default=None) – Initial guess of the value function

• ϵ (float) – Step size of the false transient method.

• tol (float) – Tolerance level of the hjb solver.

• max_iter (int) – Maximum iterations of the false transient method.

• print_iteration (bool) – Print the information of each iteration if True.

Returns

res –

dict: {

v(N,) ndarray

Value function, \(\phi(y)\).

dvdy(N,) ndarray

First order derivative of the value function, \(\frac{d\phi(y)}{dy}\).

dvddy(N,) ndarray

Second order derivative of the value function, \(\frac{d^2\phi(y)}{dy^2}\).

e_tilde(N,) ndarray

\(\tilde{e}\) on the grid of y.

πc(M, N) ndarray

Distorted probabilities of \(\{θ_\ell\}_{\ell=1}^L\).

h(N,) ndarray

Implied drift distortion.

y(N,) ndarray

Grid of y.

model_argstuple

Model parameters, see model_args above in Parameters.

}

Return type

dict of ndarrays

src.model.solve_hjb_y_jump(y, model_args=(), v0=None, ε=0.5, tol=1e-08, max_iter=10000, print_iteration=True)

Solve the HJB with y prior to jump. We impose certainty equivalent as a boundary condition to the equation instead of adding the jump intensity term.

Parameters

• y ((N,) ndarray) – An evenly spaced grid of y.

• model_args (tuple of model inputs:) –

  The ten float elements are for \(η, δ, σ_y, \bar{y}, γ_1, γ_2, γ_3, ξ_w, ξ_p, ξ_a\);

  Two (L,) ndarrays are for \(\{θ_\ell\}_{\ell=1}^L\), \(\{\pi^a_\ell\}_{\ell=1}^L\);

  Two (M,N) ndarrays are for \(\{\phi_m(y)\}_{m=1}^M\), \(\{\pi_p^m\}_{m=1}^M\).

• v0 ((N,) ndarray) – Initial guess of the value function

• ϵ (float) – Step size of the false transient method.

• tol (float) – Tolerance level of the hjb solver.

• max_iter (int) – Maximum iterations of the false transient method.

• print_iteration (bool) – Print the information of each iteration if True.

Returns

res –

dict: {

v: (N,) ndarray

Value function \(\phi(y)\)

dvdy: (N,) ndarray

First order derivative of the value function, \(\frac{d\phi(y)}{dy}\)

dvddy: (N,) ndarray

Second order derivative of the value function.

e_tilde(N,) ndarray

\(\tilde{e}\) on the grid of y.

h(N,) ndarray

Implied drift distortion.

πc(M, N) ndarray

Distorted probabilities of θ.

g(K, N) ndarray

Change in damage probability and intensity.

πd(K, N) ndarray

Distorted probabilities of damage functions.

bcfloat

The boundary condition that we impose on the HJB.

y(N,) ndarray

Grid of y.

model_argstuple

Model parameters.

}

Return type

dict of ndarrays

src.model.solve_hjb_y_jump_old(y, model_args=(), v0=None, ε=0.5, tol=1e-08, max_iter=10000, print_iteration=True)

Solve the HJB with y prior to jump. We impose certainty equivalent as a boundary condition to the equation instead of adding the jump intensity term. :param y: An evenly spaced grid of y. :type y: (N,) ndarray :param model_args: The ten float elements are for \(η, δ, σ_y, \bar{y}, γ_1, γ_2, γ_3, ξ_w, ξ_p, ξ_a\);

Two (L,) ndarrays are for \(\{θ_\ell\}_{\ell=1}^L\), \(\{\pi^a_\ell\}_{\ell=1}^L\);

Two (M,N) ndarrays are for \(\{\phi_m(y)\}_{m=1}^M\), \(\{\pi_p^m\}_{m=1}^M\).

Parameters

• v0 ((N,) ndarray) – Initial guess of the value function

• ϵ (float) – Step size of the false transient method.

• tol (float) – Tolerance level of the hjb solver.

• max_iter (int) – Maximum iterations of the false transient method.

• print_iteration (bool) – Print the information of each iteration if True.

Returns

res –

dict: {

v: (N,) ndarray

Value function \(\phi(y)\)

dvdy: (N,) ndarray

First order derivative of the value function, \(\frac{d\phi(y)}{dy}\)

dvddy: (N,) ndarray

Second order derivative of the value function.

e_tilde(N,) ndarray

\(\tilde{e}\) on the grid of y.

h(N,) ndarray

Implied drift distortion.

πc(M, N) ndarray

Distorted probabilities of θ.

g(K, N) ndarray

Change in damage probability and intensity.

πd(K, N) ndarray

Distorted probabilities of damage functions.

bcfloat

The boundary condition that we impose on the HJB.

y(N,) ndarray

Grid of y.

model_argstuple

Model parameters.

}

Return type

dict of ndarrays

src.model.solve_hjb_z(z, model_args=(), v0=None, ε=0.5, tol=1e-08, max_iter=10000, print_iteration=True)

Solve the HJB that is only related to z.

Parameters

• z ((N,) ndarray) – An evenly spaced grid of z.

• model_args (tuple of model inputs) – Elements are floats for \(η, δ, ρ, μ_2, σ_2, ξ_{1m}\) in the HJB.

• v0 ((N,) ndarray) – Initial guess of the value function

• ϵ (float) – Step size of the false transient method.

• tol (float) – Tolerance level of the hjb solver.

• max_iter (int) – Maximum iterations of the false transient method.

• print_iteration (bool) – Print the information of each iteration if True.

Returns

res –

dict: {

v(N,) ndarray

Value function, \(\zeta(z)\).

dvdz(N,) ndarray

First order derivative of the value function, \(\frac{\partial\zeta(z)}{\partial z}\).

dvddz(N,) ndarray

Second order derivative of the value function, \(\frac{\partial^2\zeta(z)}{\partial z\partial z'}\).

h(N,) ndarray

Implied drift distortion.

z(N,) ndarray

Grid of z.

model_argstuple

Model parameters.

}

Return type

dict of ndarrays

src.model.uncertainty_decomposition(y, model_args=(), e_tilde=None, h=None, πc=None, bc=None, v0=None, ε=0.5, tol=1e-08, max_iter=10000, print_iteration=True)

Solve the PDE with a given \(\tilde{e}\). If h is not given, minimize over h; if πc is not given, minimize over πc; if bc is not given, use certainty equivalent as bc.

Parameters

• y ((N, ) array) – Grid of y, \(y \in [0, \bar{y}]\).

• model_args (tuple of model inputs) – same as model_args in solve_hjb_y_jump().

• e_tilde ((N, ) array, default=None) – \(\tilde{e}\) solution from the HJB with full uncertainty configuration.

• h ((N, ) array, default=None:) – Drift distortion. To consider baseline brownian misspecification, assign np.zeros_like(y).

• πc ((L, N) array, default=None:) – Distorted probabilities of \(\{\omega_\ell\}_{\ell=1}^L\). To consider baseline smooth ambiguity, assign equal weight to coefficients for evey point of y.

• bc (tuple, default=None) – Boundary condition.

• v0 ((N, ) array, default=None) – Initial guess.

• ϵ (float, default=0.5) – Step size.

• tol (float, default=1e-8) – Tolerance level.

• max_iter (int, default=10,000) – Maximum iterations of false transient method.

• print_iteration (bool, default=True) – Print the information of each iteration if True.

Returns

res –

dict: {

v(N, ) ndarrays

Solution of value function for uncertainty decomposition.

dvdy(N, ) ndarrays

First order derivatives

dvddy(N, ) ndarrays

Second order derivatives

y(N, ) ndarray

Grid of y

ME(N, ) ndarray

Marginal value of emission.Computed according to () in

}

Return type

dict of ndarrays

simulation.py

module for simulation

src.simulation.simulate_jump(model_res, θ_list, ME=None, y_start=1, T=100, dt=1)

Simulate temperature anomaly, emission, distorted probabilities of climate models, distorted probabilities of damage functions, and drift distortion. When ME is asigned value, it will also simulate paths for marginal value of emission

Parameters

• model_res (dict) – A dictionary storing solution with misspecified jump process. See solve_hjb_y_jump() for detail.

• θ_list ((N,) ndarray:) – A list of matthew coefficients. Unit: celsius/gigaton of carbon.

• ME ((N,) ndarray) – Marginal value of emission as a function of y.

• y_start (float, default=1) – Initial value of y.

• T (int, default=100) – Time span of simulation.

• dt (float, default=1) – Time interval of simulation.

Returns

simulation_res –

dict: {

yt(T,) ndarray

Temperature anomaly trajectories.

et(T,) ndarray

Emission trajectories.

πct(T, L) ndarray

Trajectories for distorted probabilities of climate models.

πdt(T, M) ndarray

Trajectories for distorted probabilities of damage functions.

ht(T,) ndarray

Trajectories for drift distortion.

if ME is not None, the dictionary will also include

me_t(T,) ndarray

Trajectories for marginal value of emission.

}

Return type

dict of ndarrays

src.simulation.simulate_jump_2(model_res_pre, model_res_post, y_upper, θ_list, ME=None, y_start=1.1, T=100, dt=1)

Simulate temperature anomaly, emission, distorted probabilities of climate models, distorted probabilities of damage functions, and drift distortion. When ME is asigned value, it will also simulate paths for marginal value of emission

Parameters

• model_res (dict) – A dictionary storing solution with misspecified jump process. See solve_hjb_y_jump() for detail.

• θ_list ((N,) ndarray:) – A list of matthew coefficients. Unit: celsius/gigaton of carbon.

• ME ((N,) ndarray) – Marginal value of emission as a function of y.

• y_start (float, default=1) – Initial value of y.

• T (int, default=100) – Time span of simulation.

• dt (float, default=1) – Time interval of simulation.

Returns

simulation_res –

dict: {

yt(T,) ndarray

Temperature anomaly trajectories.

et(T,) ndarray

Emission trajectories.

πct(T, L) ndarray

Trajectories for distorted probabilities of climate models.

πdt(T, M) ndarray

Trajectories for distorted probabilities of damage functions.

ht(T,) ndarray

Trajectories for drift distortion.

if ME is not None, the dictionary will also include

me_t(T,) ndarray

Trajectories for marginal value of emission.

}

Return type

dict of ndarrays

src.simulation.simulate_me(y_grid, e_grid, ratio_grid, θ=0.00186, y_start=1, T=100, dt=1)

simulate trajectories of uncertainty decomposition

\[\log(\frac{ME_{new}}{ME_{baseline}})\times 1000.\]

Parameters

• y_grid ((N, ) ndarray) – Grid of y.

• e_grid ((N, ) ndarray) – Corresponding \(\tilde{e}\) on the grid of y.

• ratio_grid ((N, ) ndarray:) – Corresponding \(\log(\frac{ME_{new}}{ME_{baseline}})\times 1000\) on the grid of y.

• θ (float, default=1.86/1000) – Coefficient used for simulation.

• y_start (floatsimulation) – Initial value of y.

• T (int, default=100) – Time span of simulation.

• dt (float, default=1) – Time interval of simulation. Default=1 indicates yearly simulation.

Returns

• Et ((T, ) ndarray) – Emission trajectory.

• yt ((T, ) ndarray) – Temperature anomaly trajectories.

• ratio_t ((T, ) ndarray) – Uncertainty decomposition ratio trajectories.

solver.py

Functions to numerically solve a nonlinear ODE via false transient scheme.

See also

For more details on false transient: refer to appendix B (TODO:fix to ref).

The ODE we solve is as follows:

The state space descretize into \(y_1, y_2, \dots, y_N\) with equal interval \(\Delta y\). For each \(y_n \in \{y_1, y_2, \dots, y_N\}\) in the grid with , the value function satifies (we use upwinding first order derivative for concern of, and at \(y_1\) we use forward derivative instead):

(1)\[\begin{align} \frac{\phi_{i+1}(y_n) - \phi_{i}(y)}{\epsilon} =& \quad A_n \phi_{i+1}(y) \cr & + B_{n} \frac{\phi_{i+1}(y_n) - \phi_{i+1}(y_{n-1})}{\Delta y} \cr & + C_n \frac{\phi_{i+1}(y_{n+1}) - 2 \phi_{i+1}(y_{n}) + \phi_{i+1}(y_{n-1})}{\Delta y^2} \cr & + D_n \end{align}\]

where \(A_n\), \(B_n\), \(C_n\) and \(D_n\) are coefficients. The exact values are given by the HJB to be solved. Therefore, we construct the following linear system:

(2)\[LHS \cdot \phi_{i+1}(Y) = - D - \frac{1}{\epsilon} \phi_{i}(Y)\]

where \(LHS\) is a \(N\times N\) matrix of coefficients, \(Y = (y_1, y_2, \dots, y_N)'\), and \(D = (D_1, D_2, \dots, D_n)'\).

src.solver.compute_coefficient(LHS, A, B, C, i, dx, ε)

Compute the coefficient of the equation at \(y_i\). Values are computed according to (1).

Parameters

• LHS ((I, I) ndarray) – LHS matrix of the linear system.

• A ((I,) ndarrays) – Coefficient arrays for value function \(\phi(y_i)\).

• B ((I,) ndarrays) – Coefficient arrays for first order derivative.

• C ((I,) ndarrays) – Coefficient arrays for second order derivative.

• i (int) – Index for the grid point.

• dx (float) – Grid step, \(\Delta y\).

• ϵ (float) – False transient step size.

Returns

LHS – Updated LHS matrix of the linear system.

Return type

(I, I) ndarray

src.solver.false_transient(A, B, C, D, v0, ε, dx, bc, impose_bc)

Implement false transient scheme of one iteration, and update from \(\phi_i(y)\) to \(\phi_{i+1}(y)\) according to (2).

See appendix B(TODO: cross-ref link) for further detail.

Parameters

• A – Same as those in compute_coefficient()

• B – Same as those in compute_coefficient()

• C – Same as those in compute_coefficient()

• D – Same as those in compute_coefficient()

• v0 ((I, ) ndarrays) – Value function from last iteration, \(\phi_i(y)\)

• ϵ (float) – Step size of false transient

• dx (float) – Step size of Grid

• bc – See linearize().

• impose_bc – See linearize().

Returns

v – Updated value function, \(\phi_{i+1}(y)\) according to (2).

Return type

(I, ) ndarrays

src.solver.linearize(A, B, C, D, v0, ε, dx, bc, impose_bc)

Construct coefficient matrix of the linear system.

Parameters

• A ((I,) ndarrays) – see compute_coefficient() for further detail.

• B ((I,) ndarrays) – see compute_coefficient() for further detail.

• C ((I,) ndarrays) – see compute_coefficient() for further detail.

• D ((I,) ndarrays) – see compute_coefficient() for further detail.

• v0 ((I,) ndarray) – Value function from last iteration, \(\phi_i(y)\).

• ϵ (float) – False transient step size.

• dx (float) – Grid step size.

• bc (tuple of ndarrays:) –

  Impose v=bc[k] at boundaries.

  Order: lower boundary of x, upper boundary of x,

• impose_bc (tuple of bools) – Order: lower boundary of x, upper boundary of x,

Returns

• LHS ((I, I) ndarray) – LHS of the linear system.

• RHS ((I,) ndarray) – RHS of the linear system.

utilities.py

Functions that facilitate computation.

src.utilities.compute_derivatives(data, order, dlt)

Compute the derivatives of a function.

Parameters

• data ((N, ) ndarrays) – The function whose derivatives to be computed.

• order (int:) –

  1: first order derivative,

  2: second order derivative.

• dlt (float) – Grid step.

Returns

res – Deriative array.

Return type

(N, ) ndarrays

Notes

First order derivatives are computed one-sidedly. Specifically,

\[\begin{split}\frac{d f(x_0)}{dx} \approx & \frac{f(x_1) - f(x_0)}{\Delta x} \\ \frac{d f(x_i)}{dx} \approx & \frac{f(x_i) - f(x_{i-1})}{\Delta x}, i = 1,2,\dots,N\end{split}\]

src.utilities.find_nearest_value(array, value)

Return the index of the element in array that is closest to value.