\[\begin{split}\begin{align*}
d \log N_t =
\begin{cases}
\left( \lambda_1 + \lambda_2 Y_t \right) \mathcal{E}_t \left[ \bar{\theta}\, dt + \varsigma\, dW_t \right] + \dfrac{ \lambda_2 |\varsigma|^2 \left( \mathcal{E}_t \right)^2 }{2}\, dt, & \text{if } t \leq \tau \\[2ex]
\left[ \lambda_1 + \lambda_2 \widehat{Y}_t + \lambda_3(\ell)\left( \widehat{Y}_t - \bar{y} \right) \right] \mathcal{E}_t \left[ \bar{\theta}\, dt + \varsigma\, dW_t \right] \\
\quad + \dfrac{ \left[ \lambda_2 + \lambda_3(\ell) \right] |\varsigma|^2 \left( \mathcal{E}_t \right)^2 }{2}\, dt, & \text{if } t > \tau
\end{cases}
\end{align*}\end{split}\]
1.1 HJB equations
The value function can be written as
\[{V}(X_t) = {\hat V}(X_t^1) - {\hat N}_t\]
where it is straightforward to verify the additive separability in the
logarithm of damages. This separability simplifies our numerical
solutions.
1.1.1 Post-tech-post-damage HJB
The optimal choice of emission is greater than \(\beta \alpha K_t\),
after technology jump occurs. As temperature and damage remains
constant, emission only affect the production function. The optimal
production function is:
\[\max_{{\mathcal E}_t} \alpha K_t \left (1 - \phi_0\left( \left(1 - \frac {{\mathcal E}_t}{\beta \alpha K_t} \right){\mathbf 1}_{\{{\mathcal E}_t < \beta \alpha K_t\} } \right)^{\phi_1} \right) = \alpha K_t\]
In this case, the only state variable, control space and distortion
space become
\[\begin{split}\begin{align*}
x &= \{ \hat{k} \}\\
\Phi &= \{ i^k \}\\
\Gamma &= \{h^k\}
\end{align*}\end{split}\]
Compute \(\hat{v}^{\ell,L}\) for \(\ell = 1, ..., L-1\)
conditioned on both a technology jump and a damage jump occurring,by
solving HJB equation
\[\begin{split}\begin{align*}
0= & \max_{i^k}\min_{{h^k}} \left(\frac{\delta}{1-\rho}\right)\left[\left(\frac{\alpha k -i^k}{\exp (\hat{v}^{\ell,L})} \right)^{1-\rho}-1\right] \\
& +\frac{\partial \hat{v}^{\ell,L}}{\partial \hat{k}}\left[\mu_k+\frac{i^k}{k}-\frac{\kappa}{2} \left(\frac{i^k}{k}\right)^2-\frac{\left|\sigma_k\right|^2}{2}+\sigma_k {h^k}\right]+\frac{\partial^2 \hat{v}^{\ell,L}}{\partial \hat{k} \, \partial \hat{k}'}\frac{\left|\sigma_k\right|^2}{2} \\
& +\xi^k \frac{\left|{h^k}\right|^2}{2}
\end{align*}\end{split}\]
If \(\rho =1\), the first term of the HJB becomes
\[\delta \log ( \alpha k -i^k ) - \delta \hat{v}^{\ell,L}\]
1.1.2 Post Technology and Pre Damage HJB
Compute the value function \(\hat{v}^L\) assuming that only a
technology jump has been realized. This value function incorporates the
possibility of jumping to one of \(L-1\) possible damage states,
however, because the temperature anomaly remains constant at the point
where any incremental curvature has no impact on the damages this damage
curve realization is inconsequential. Therefore, we can ignore the
damage curve intensities and the associated continuation values for this
computation.
In this case, the state space, control space and distortion space become
\[\begin{split}\begin{align*}
x &= \{ \hat{k} \}\\
\Phi &= \{ i^k \}\\
\Gamma &= \{h^k\}
\end{align*}\end{split}\]
We have HJB for post technology and post damage jump as follows
\[\begin{split}\begin{align*}
0= & \max_{i^k}\min_{{h^k}} \left(\frac{\delta}{1-\rho}\right)\left[\left(\frac{\alpha k -i^k}{\exp (\hat{v}^L)} \right)^{1-\rho}-1\right] \\
& +\frac{\partial \hat{v}^{L}}{\partial \hat{k}}\left[\mu_k+\frac{i^k}{k}-\frac{\kappa}{2} \left(\frac{i^k}{k}\right)^2-\frac{\left|\sigma_k\right|^2}{2}+\sigma_k {h^k}\right]+\frac{\partial^2 \hat{v}^{ L}}{\partial \hat{k} \, \partial \hat{k}'} \frac{\left|\sigma_k\right|^2}{2} \\
& +\xi^k \frac{\left|{h^k}\right|^2}{2}
\end{align*}\end{split}\]
1.1.3 Pre Technology and Post Damage HJB
Compute the values functions \(\hat{v}^{\ell}\) assuming that only a
damage jump has been realized for \(\ell = 1,..., L-1.\) These
values functions depend on the entire state vector \(X\) and have
one possible jump state which is the technology discovery with intensity
\({\mathcal J}^L\). The continuation value for the jump is
\(\hat{v}^{\ell,L}\) viewed as a function of \(x\) for
\(\ell=1,...,L-1.\)
\[\begin{split}\begin{align*}
x &= \{ \hat{k}, y, \hat{r}, \hat{n} \}\\
\Phi &= \{ i^k, i^r, e \}\\
\Gamma &= \{{h^k}, {h^y}, {h^r}, g\}
\end{align*}\end{split}\]
After plugging this simplification into our HJB equation and removing
common terms, we are left with the following simplified HJB to solve:
\[\begin{split}\begin{align*}
& 0=\max_{i^k, i^r, e} \min_{{h^k}, {h^y}, {h^r}, g} \left(\frac{\delta}{1-\rho}\right)\left[\left(\frac{\alpha k -i^k-i^r-\alpha k \phi_0(z)\left[1-\frac{e}{\beta_t \alpha k }\right]^{\phi_1}}{\exp (\hat{v}^\ell)} \right)^{1-\rho}-1\right] \\
& +\frac{\partial \hat{v}^\ell}{\partial \hat{k}}\left[\mu_k+\frac{i^k}{k}-\frac{\kappa}{2} \left(\frac{i^k}{k}\right)^2-\frac{\left|\sigma_k\right|^2}{2}+\sigma_k {h^k}\right]+\frac{\partial^2 \hat{v}^\ell }{\partial \hat{k} \partial \hat{k}'} \frac{\left|\sigma_k\right|^2}{2} \\
& +\frac{\partial \hat{v}^\ell}{\partial \hat{y}}\left( \bar{\theta}+\varsigma {h^y}\right) e+\frac{\partial^2 \hat{v}^\ell}{\partial y \partial y'} \frac{|\varsigma|^2}{2} e^2 \\
& -\left(\left[\lambda_1+\lambda_2 y+\lambda_3(y-\bar{y})\right]\left( \bar{\theta}+\varsigma {h^y}\right) e+\left(\lambda_2+\lambda_3\right) \frac{|\varsigma|^2}{2} e^2\right) \\
& +\frac{\partial \hat{v}^\ell}{\partial \hat{r} }\left(-\zeta+\psi_0\left(i^r\right)^{\psi_1} \exp \left(-\psi_1 \log r\right)-\frac{\left|\sigma_r\right|^2}{2}+\sigma_r {h^r}\right)+\frac{\partial^2 \hat{v}^\ell}{\partial \hat{r} \partial \hat{r}'}\frac{\left|\sigma_r\right|^2}{2} \\
& +\xi^g \mathcal{J}_g (1-g +g \log g )+\mathcal{J}_g g \left(\hat{v}^{\ell,L}-\hat{v}^\ell \right) \\
& +\xi^k \frac{\left|{h^k}\right|^2}{2}+\xi^c \frac{\left|{h^y}\right|^2}{2}+\xi^r \frac{\left|{h^r}\right|^2}{2} \\
&
\end{align*}\end{split}\]
If \(\rho =1\), the first term of the HJB becomes
\[\delta \log ( \alpha k -i^k-i^r-\alpha k \phi_0(z)\left[1-\frac{e}{\beta_t \alpha k }\right]^{\phi_1} ) - \delta \hat{v}^\ell\]
1.1.4 Pre-tech-pre-damage HJB
Compute \(\hat{v}\) prior to any jumps occurring. This value
function has two possible types of jumps, either a technology jump or a
damage curvature jump. The continuation value for the technology jump is
\(\hat{v}^L\), and the potential continuation values for the damage
curvature jump are the set of \(\hat{v}^{\ell}\) for
\(\ell = 1,..., L-1.\)
\[\begin{split}\begin{align*}
x &= \{ k, y,r,n \}\\
\Phi &= \{ i^k, i^r, e \}\\
\Gamma &= \{{h^k}, {h^y}, {h^r}, g, f\}
\end{align*}\end{split}\]
After plugging this simplification into our HJB equation and removing
common terms,
\[\begin{split}\begin{align*}
0 = & \max_{i^k, i^r, e} \, \min_{{h^k}, {h^y}, {h^r}, g, f} \, \frac{\delta}{1-\rho} \left(\left(\frac{\alpha k-i^{k}-i^{r}-\alpha k \phi_0 \left(1-\frac{e}{\beta \alpha k}\right)^{\phi_1}}{\exp(\hat{v})} \right)^{1-\rho}-1 \right) \\
& + \frac{\partial \hat{v}}{\partial \hat{k}} \left( -\mu_{k}+ \frac{i^{k}}{k}-\frac{\kappa}{2}\left(\frac{I^{k}}{k}\right)^{2}-\frac{|\sigma_{k}|^{2}}{2} + \sigma_k h^k \right) + \frac{\partial^2 \hat{v}}{\partial \hat{k} \, \partial \hat{k}'}\frac{|\sigma_{k}|^{2}}{2} \\
& + \frac{\partial \hat{v}}{\partial y} e \left( \bar{\theta}+\varsigma h^y \right) + \frac{\partial^2 \hat{v}}{\partial y \, \partial y'}\frac{|\varsigma|^{2}}{2}e^{2} - \left( (\lambda_{1}+\lambda_{2}y) e \left( \bar{\theta}+\varsigma h^y \right) +\lambda_{2}\frac{|\varsigma|^{2}}{2}e^{2} \right) \\
& + \frac{\partial \hat{v}}{\partial \hat{r}} \left( -\zeta + \psi_{0}(i^{r})^{\psi_{1}}\exp( -\psi_{1} \hat{r})-\frac{|\sigma_{r}|^{2}}{2}+\sigma_{r} h^r \right) +\frac{\partial^2 \hat{v}}{\partial \hat{r} \, \partial \hat{r}'}\frac{|\sigma_{r}|^{2}}{2} \\
& +\xi^g \mathcal{J}_g (1-g +g \log g )+\mathcal{J}_g \cdot g \cdot \left(\hat{v}^L -\hat{v}\right) \\
&+\xi^d \mathcal{J}_n \sum_{\ell} \pi^\ell (1-f^\ell +f^\ell \log f^\ell ) \\
&+\mathcal{J}_n \sum_{\ell } \pi^\ell f^\ell \cdot \left(\hat{v}^\ell-\hat{v}\right) \\
&+\xi^k \frac{\left|{h^k}\right|^2}{2}+\xi^c \frac{\left|{h^y}\right|^2}{2}+\xi^r \frac{\left|{h^r}\right|^2}{2}
\end{align*}\end{split}\]
To solve HJB equations, we first run below code in
two-capital-climate-change/master
/master_zero_shock.sh.
Make sure you give right command-line arguments.
We solve four types of HJB equations sequentially.
First, solve one post-tech-post-damage HJB. As after technology jump
occurs, the curvature of damage function does not appear in HJB
equations.
Second we solve one post-tech-pre-damage and twenty
pre-tech-post-damage HJB conditional on post-tech-post-damage value
function.
Finally, we solve pre-tech-pre-damage HJB given post-tech-pre-damage
and pre-tech-post-damage value functions.
In
Postdamage.sh,
we solve post_damage_post_tech and post-damage-pre-tech value functions
and controls.
Post_damage_post_tech
section solves post-damage-post-tech HJB.
Post-damage-pre-tech
section solves Post-damage-pre-tech HJB. In order to make sure our
results are stable, we first randomly pick initial values and then use
the first result to resolve the HJB.
Postdamage_sub.sh
is aimed at further improving computational efficiency. The solutions
obtained from post_damage.py serve as baseline solutions for
Postdamage_sub.py to resolve the HJB equations.
In
Predamage.sh,
we solve pre_damage_post_tech and pre-damage-pre-tech value functions
and controls.
Pre_damage_post_tech
section solves pre-damage-post-tech HJB.
Pre-damage-pre-tech
section solves Pre-damage-pre-tech HJB.
1.2 Computation method
In this section, we explain how did we solve HJB equation.
1.2.1 Policy Iteration
For simplicity, I denote the control set and distortion set:
\[\begin{split}\begin{align*}
\Phi^n &= \{ i_k^{n}, i_j^{n}, \mathcal{E}^{n} \} \\
\Gamma^n &=\{ h_k^{n}, h_y^{n}, h_j^{n}, g^{n}, f_\ell^{n} \}
\end{align*}\end{split}\]
Algorithm: Solving the HJB Equation via Policy Iteration
\[\begin{split}\begin{align*}
\textbf{Input:} &\ \text{Initial guess for value function } \hat{v}^0, \epsilon = 10^{-7} \\
&\text{Initialize } n = 0, \hat{v}^n = \hat{v}^0 \\
\textbf{while} &\ |\hat{v}^{n+1} - \hat{v}^n| \geq \epsilon \text{ do:} \\
&\ \quad \text{Step 1: Solve for optimal actions} \Phi^{n+1} \text{ by maximization} \\
&\ \quad \quad \text{Cobweb algorithm is applied here:} \\
&\ \quad \quad \Phi^{n+1} = \Phi(\hat{v}^n, \Phi^{n}, \Gamma^{n}) \\
&\ \quad \text{Step 2: Solve for optimal probability distortions } \Gamma^{n+1} \text{ by minization}\\
&\ \quad \quad \Gamma^{n+1} = \Gamma(\hat{v}^n, \Phi^{n+1}, \Gamma^{n}) \\
&\ \quad \text{Step 3: Update value function } \hat{v}^{n+1} \text{ by minimization}\\
&\ \quad \quad \hat{v}^{n+1} = V(\hat{v}^n, \Phi^{n+1}, \Gamma^{n+1}) \\
&\ \quad \text{Step 4: Check for convergence} \\
&\ \quad \quad\text{If } |\hat{v}^{n+1} - \hat{v}^n| < \epsilon \text{ then stop, otherwise continue.} \\
\textbf{Return:} &\ \hat{v}^* \\
\end{align*}\end{split}\]
1.2.2 Updating Rules \(\Phi^{n+1} = \Phi(\hat{v}^n,\Phi^{n},\Gamma^{n})\)
In solving HJB equations, we often encounter complex, highly non-linear
equations that do not admit analytical solutions. To address this
challenge, iterative numerical methods like the
Cobweb
algorithm are employed to approximate the optimal control variables.
The Cobweb algorithm works by:
Starting with an initial guess for the control variable.
Computing the corresponding values in the equations.
Updating the control variable based on the discrepancies observed.
Repeating the process until the control variable converges to a
stable value.
For example, we update for \(i_k\) for pre damage pre technology
HJB, using the first-order condition:
\[\delta \left( \frac{\alpha k - i_k - i_j - \alpha k \phi_0(z) \left[1 - \frac{\mathcal{E}}{\beta_t \alpha k}\right]^{\phi_1}}{\exp(\hat{v})} \right)^{-\rho} \frac{1}{\exp(\hat{v})} = \frac{\partial \hat{v}}{\partial \log k} \left(1 - \kappa i_k\right)\]
Since this equation is highly non-linear and does not admit an
analytical solution, we use the
Cobweb
algorithm to iteratively update the actions. For each iteration
\(n\), the update is:
\[\begin{align}
\hat{i}_k^{n+1} = \frac{1}{\kappa}-\frac{1}{\kappa}\delta \left( \frac{\alpha k - i_k^n - i_j - \alpha k \phi_0(z) \left[1 - \frac{\mathcal{E}}{\beta_t \alpha k}\right]^{\phi_1}}{\exp(\hat{v})} \right)^{-\rho} \frac{1}{\exp(\hat{v})} \frac{1}{\frac{\partial \hat{v}}{\partial \log k}} \end{align}\]
The updated
action
\(i_k^{n+1}\) is computed using a relaxation parameter \(\chi\):
\[i_k^{n+1} = \chi i_k^n + (1 - \chi) \hat{i}_k^{n+1}\]
1.2.3 Updating Rules \(\Gamma^{n+1} = \Gamma(\hat{v}^n,\Phi^{n+1},\Gamma^{n} )\)
Every distortion has analytical solution. For example, we solve for
\(h_k\), and the same logic applies to \(h_y, h_j, g, f_l\). The
first-order condition for \(h_k\) is:
\[\frac{\partial \hat{v}}{\partial \log k} \sigma_k = - \xi_k h_k\]
Given the value function \(v^n\), we
update
the distortion \(h_k^{n+1}\) as follows:
\[h_k^{n+1} = - \frac{1}{\xi_k} \frac{\partial \hat{v}^n}{\partial \log k} \sigma_k\]
1.2.4 Solve Linear PDE Equation
Updating value functions, given the state variables and controls, is
solving a linear PDE system. To mitigate the potential instability of
the non-linear HJB, we add a false transcient (time) dimension and solve
it until convergence. Here we use Petsc
to solve the PDE system, so we show how to rewrite the PDE and call
Petsc package.
For example, in pre-tech-pre-damage case with \(\rho\neq 1\), we can
write the HJB into the form:
\[A \hat{v}
+B_{\hat{k}} \frac{\partial \hat{v}}{\partial \hat{k}}
+B_{y}\frac{\partial \hat{v}}{\partial y}
+B_{\hat{r}} \frac{\partial \hat{v}}{\partial \hat{r}}
+C_{\hat{k}} \frac{\partial^2 \hat{v}}{\partial \hat{k} \, \partial \hat{k}'}
+C_{y} \frac{\partial^2 \hat{v}}{\partial y \, \partial y'}
+C_{\hat{r}} \frac{\partial^2 \hat{v}}{\partial \hat{r} \, \partial \hat{r}'}
+D =0\]
First and Second order partial derivatives
\[\left\{\frac{\partial \hat{v}}{\partial \hat{k}},\frac{\partial \hat{v}}{\partial y}, \frac{\partial \hat{v}}{\partial \hat{r}}\right\}, \quad, \left\{ \frac{\partial^2 \hat{v}}{\partial \hat{k} \, \partial \hat{k}'}, \frac{\partial^2 \hat{v}}{\partial y \, \partial y'}, \frac{\partial^2 \hat{v}}{\partial \hat{r} \, \partial \hat{r}'} \right\}\]
The coefficient before Value function:
\[A = - \mathcal{J}_g \cdot g-\mathcal{J}_n \sum_{\ell } \pi^\ell f^\ell\]
Coefficient of first order partial derivatives:
\[B_{\hat{k}} = -\mu_{k}+ \frac{i^{k}}{k}-\frac{\kappa}{2}\left( i^{k} \right)^{2}-\frac{|\sigma_{k}|^{2}}{2} + \sigma_k h^k\]
\[B_{y} =e \left( \bar{\theta}+\varsigma h^y \right)\]
\[B_{\hat{r}} = -\zeta + \psi_{0}(i^{r})^{\psi_{1}}\exp( -\psi_{1} \hat{r})-\frac{|\sigma_{r}|^{2}}{2}+\sigma_{r} h^r\]
Coefficient of second order partial derivatives:
\[C_{\hat{k}} = \frac{|\sigma_{k}|^{2}}{2},\quad C_{y} = \frac{|\varsigma|^{2}}{2}e^{2},\quad C_{\hat{r}} = \frac{|\sigma_{r}|^{2}}{2}\]
\[\begin{split}\begin{align*}
D = & \frac{\delta}{1-\rho} \left(\left(\frac{\alpha k-i^{k}-i^{r}-\alpha k \phi_0 \left(1-\frac{e}{\beta \alpha k}\right)^{\phi_1}}{\exp(\hat{v})} \right)^{1-\rho}-1 \right) \\
& - \left( (\lambda_{1}+\lambda_{2}y) e \left( \bar{\theta}+\varsigma h^y \right) +\lambda_{2}\frac{|\varsigma|^{2}}{2}e^{2} \right) \\
& +\xi^g \mathcal{J}_g (1-g +g \log g )+\mathcal{J}_g \cdot g \cdot \hat{v}^L \\
&+\xi^d \mathcal{J}_n \sum_{\ell} \pi^\ell (1-f^\ell +f^\ell \log f^\ell ) \\
&+\mathcal{J}_n \sum_{\ell } \pi^\ell f^\ell \cdot \hat{v}^\ell \\
&+\xi^k \frac{\left|{h^k}\right|^2}{2}+\xi^c \frac{\left|{h^y}\right|^2}{2}+\xi^r \frac{\left|{h^r}\right|^2}{2}
\end{align*}\end{split}\]
1.2.5 Finite Difference Schemes
\[\begin{split}\begin{align*}
\left(\frac{\partial f}{\partial x}\right)_i = \frac{f_{i+1} - f_{i-1}}{2 \Delta x} \\
\left(\frac{\partial^2 f}{\partial x^2}\right)_i =\frac{f_{i+1} + f_{i-1} - 2f_i}{\Delta x^2}
\end{align*}\end{split}\]
\[\begin{split}\begin{align*}
\left(\frac{\partial f}{\partial x}\right)_0 =\frac{f_{1} - f_{0}}{\Delta x} \\
\left(\frac{\partial^2 f}{\partial x^2}\right)_0 =\frac{f_{2} + f_{0} - 2f_{1}}{\Delta x^2}
\end{align*}\end{split}\]
\[\begin{split}\begin{align*}
\left(\frac{\partial f}{\partial x}\right)_{N-1} =\frac{f_{N-1} - f_{N-2}}{\Delta x} \\
\left(\frac{\partial^2 f}{\partial x^2}\right)_{N-1}=\frac{f_{N-1} + f_{N-3} - 2f_{N-2}}{\Delta x^2}
\end{align*}\end{split}\]
Below two functions are two finite difference functions we used in
solving HJB equations.