HJBE_discretization.tex

\documentclass[11pt]{etk-article}
\usepackage{pstool} 
\usepackage{etk-bib}
\pdfmetadata{}{}{}{}
\externaldocument[FD:]{operator_discretization_finite_differences}

\begin{document}
\title{Solving HJBE with Finite Differences}
\author{Jesse Perla\\UBC}
\date{\today}
\maketitle
 Here, we expand on details for how to discretize HJBE with a control of the drift.\footnote{See \url{operator_discretization_finite_differences.pdf} for more details on the discretization of a linear diffusion operator with finite-differences, and general notation (with equation numbers in that document prefaced by \textit{FD}). Thanks to Sari Wang for superb Research Assistance.}  In particular, this set of notes will focus on solving the neoclassical growth model first with a deterministic and then a stochastic TFP.
 
 \section{Neoclassical Growth}
 This solves the simple deterministic neoclassical growth model.
 \subsection{Value Function with Capital Reversibility}
 
 Take a standard neoclassical growth model with capital $k$, consumption $c$, production $f(k)$, utility $u(c)$, depreciation rate $\delta$, and discount rate $\rho$.\footnote{
 This builds on \url{http://www.princeton.edu/~moll/HACTproject/HACT_Additional_Codes.pdf} Section 2.2, with code in \url{http://www.princeton.edu/~moll/HACTproject/HJB_NGM_implicit.m}}
The law of motion for capital in this setup is,
\begin{align}
	\D[t] k(t) &= f(k) - \delta k - c\label{eq:lom-capital}\\
	\intertext{With this, the standard HJBE for the value of capital $k$ is,}
	\rho v(k) &= \max_{c}\set{u(c) + \left(f(k) - \delta k - c\right)v'(k)}\label{eq:HJBE-neoclassical-growth}
	\intertext{Assume an interior $c$ and envelope conditions, then taking the FOC of the choice is,}
	u'(c) &= v'(k)
	\intertext{Assume a functional form $u(c) = \frac{c^{1-\gamma}}{1-\gamma}$ for the utility and $f(k)=Ak^\alpha$ for production. Then the FOC can be inverted such that}
	c &= \left(v'(k) \right)^{-\frac{1}{\gamma}}\label{eq:c-neoclassical-growth}
	\intertext{And,}
	u(c) &= \frac{\left(v'(k)\right)^{-\frac{1-\gamma}{\gamma}}}{1-\gamma}\label{eq:u-c-neoclassical-growth}
\intertext{If \cref{eq:c-neoclassical-growth,eq:u-c-neoclassical-growth} were substituted back into \cref{eq:HJBE-neoclassical-growth}, we would have a nonlinear ODE in just $k$.  Define the following, which will be important in the finite difference scheme}
	\mu(k) &\equiv f(k) - \delta k - c\label{eq:mu-neoclassical-growth}
	\intertext{Then with \cref{eq:c-neoclassical-growth} this can be defined as a function of the $v(\cdot)$ function,}
	\mu(k;v) &\equiv f(k) - \delta k - \left(v'(k) \right)^{-\frac{1}{\gamma}}\label{eq:mu-v-neoclassical-growth}
	\intertext{Then with \cref{eq:HJBE-neoclassical-growth,eq:mu-v-neoclassical-growth,eq:u-c-neoclassical-growth}}
	\rho v(k) &= \frac{\left(v'(k)\right)^{-\frac{1-\gamma}{\gamma}}}{1-\gamma} + \mu(k;v) v'(k)\label{eq:HJBE-neoclassical-growth-mu}
\end{align}
The HJBE \cref{eq:HJBE-neoclassical-growth-mu} is a nonlinear ODE in $v(k)$.

\paragraph{Steady State}  The goal is for the finite-difference scheme to find the steady-state on its own.  However, in this example we have an equation to verify the solution:
\begin{align}
k^* &= \left(\frac{\alpha A}{\rho + \delta}\right)^{\frac{1}{1-\alpha}}\label{eq:steady-state-k}\\
c^* &= A \left(k^*\right)^{\alpha} - \delta  \left(k^*\right)^{\alpha}\label{eq:steady-state-c}
\end{align}
\subsection{Finite Differences and Discretization}
In order to solve this problem, we will need to use the appropriate ``upwind'' direction for the finite differences given the sign of the drift in \cref{eq:mu-v-neoclassical-growth}.  
\paragraph {Two basic approaches:} 
\begin{itemize}
	\item Fix the nonlinear $v'(k)$ function so that the operator becomes linear, solve it as a sparse linear system, and then iterate on the $v(k)$ solution.
	\item Solve it directly as a nonlinear system of equations.
\end{itemize}

\paragraph {Setup}
\begin{itemize}
\item Define a uniform grid with $I$ discrete points  $\set{k_i}_{i=1}^I$,  for some small $k_1 \in (0, k^*)$ and large $k_I > k^*$, with distance between grid points $\Delta \equiv k_i - k_{i-1}$ for all $i$. After discretizing, we will denote the grid with the variable name, i.e. $k \equiv \set{k_i}_{i=1}^I$.
Further define the notations $\underline{k} \equiv k_1$, $\overline{k} \equiv k_I$ and $v_i \equiv v(k_i)$ for simplicity. 
\item When we discretize a function, use the function name without arguments to denote the vector.  i.e. $v(k)$ discretized on a grid $\set{k_i}_{i=1}^{I}$ is $v \equiv \set{v(k_i)}_{i=1}^I \in \R^I$.
\item When referring to a variable $\mu$, define the notation $\mu^{-} \equiv \min\set{\mu,0}$ and $\mu^{+} \equiv \max\set{\mu,0}$. This can apply to vectors as well. For example, $\mu_i^{-} = \mu_i$ if $\mu_i < 0$ and $0$ if $\mu_i > 0$, and $\mu^{-} \equiv \set{\mu^{-}_i}_{i=1}^{I}$.

\item To discretize the derivative at $k_i$, consider both forwards and backwards differences,
\begin{align}
	v'_F(k_i) &\approx \frac{v_{i+1} - v_i}{\Delta}\label{eq:forward-diff}\\
	v'_B(k_i) &\approx \frac{v_i - v_{i-1}}{\Delta}\label{eq:backward-diff}
	\intertext{subject to the state constraints} 
		v'_{F}(\overline{k}) = (f(\overline{k}) - \delta \overline{k})^{-\gamma}\label{eq:forward-constraint}\\
		v'_{B}(\underline{k}) = (f(\underline{k}) - \delta \underline{k})^{-\gamma}\label{eq:backward-constraint}
	\intertext {Our finite difference approximation to the HJB equation \cref{eq:HJBE-neoclassical-growth-mu} is then}
	\rho v(k_i) &= \frac{\left(v'(k_i)\right)^{-\frac{1-\gamma}{\gamma}}}{1-\gamma} + \mu(k_i;v) v'(k_i)
\end{align}
\end{itemize}
 Following the ``upwind scheme'', we will use the forwards approximation whenever there is a positive drift above the steady state level of capital $k^*$, and the backwards approximation when below $k^*$. 

\subsection{Iterating on a Linear System of Equations}
This section solves the problem through a series of iterations on a linear system of equations.
\paragraph{Iterative Method }Pick an $h$, which is used in the iterative process.\footnote{\textbf{TODO:} I believe this is the Howard algorithm.  In effect, I think that this is like using an explicit time step in a PDE, where we are forward iterating fixing the policy?  Because of this, there is a high likelihood that it is not unconditionally stable.}  We first make an initial guess $v^0 = (v_1^0, \dots v_I^0)$. Then for each consecutive iteration $n = 0, 1, 2, ...$,  update $v^n$ using the following equation
\begin{align}
\frac{v_i^{n+1}-v_i^{n}}{h} + \rho v_i^{n+1} = u(c_i^n) + (v_i^{n+1})^{'}  \underbrace{\left[f(k_i) - \delta k_i - c_i^n\right]}_{\equiv \mu_i}\label{eq:fd-approx-capital}
\end{align}
Using both the forwards and backwards difference approximations \cref{eq:forward-diff} and \cref{eq:backward-diff}, find consumption levels $c_F^n$ and $c_B^n$ as given by \cref{eq:c-neoclassical-growth} and compute savings,
\begin{align}
\mu^n_{i,F} = f(k_i) - \delta k_i - \left({v_F^n}'(k_i) \right)^{-\frac{1}{\gamma}}\\
\mu^n_{i,B} = f(k_i) - \delta k_i - \left({v_B^n}'(k_i) \right)^{-\frac{1}{\gamma}}
\end{align}

Depending upon the sign of the drift, a choice is made between using forward or backward differences. Let $\bold{1}_{\{ \cdot \}}$ be an indicator function\footnote{For now, assume that if the case $\mu_{F}>0$  and $\mu_{B}<0$ should arise, take $\bold{1}_{\{ \mu_{i,F} > 0\}} $ to be 1 and $\bold{1}_{ \{\mu_{i,B} < 0\}}\ $ to be 0.}and  $\bar{v_i}' $ be the derivative at the steady state, given in this example by $\bar{v_i}' = (f(k_i)-\delta k_i)^{-\gamma}$. Then the approximation of the derivative is 
\begin{align}
{v_i^n}'={v_{i,F}^n}'\bold{1}_{\{ \mu_{i,F} > 0\}}+{v_{i,B}^n}\bold{1}_{\{ \mu_{i,B}<0 \}}+{\bar{v}_i^n}'\bold{1}_{\{ \mu_{i,F}<0<\mu_{i,B} \}}\label{eq:dv-upwind}
\end{align}

Define the vectors $X, Y, Z \in \R^{I} $ such that 
\begin{align}
	X &= -\frac {\mu^{-} _B}{\Delta}\label{eq:X-delta} \\
	Y &= -\frac {\mu^{+} _F}{\Delta} + \frac {\mu^{-} _B}{\Delta}\label{eq:Y-delta} \\
	Z &= \frac {\mu^{+} _F}{\Delta}\label{eq:Z-delta}
\end{align}
%
%%For algebraic simplicity later on, multiply every term by  $\Delta$ and obtain 
%\begin{align}
%	X &= -({\mu^{n} _B})^{-}\label{eq:X} \\
%	Y &= -({\mu^{n} _F})^{+} + ({\mu^{n} _B})^{-}\label{eq:Y} \\
%	Z &= ({\mu^{n} _F})^{+}\label{eq:Z}
%\end{align}

With these, construct the sparse matrix $A^n$
\begin{align}
A^n &\equiv \begin{bmatrix}
Y_1 & Z_1 & 0 & \cdots & \cdots & \cdots & 0 \\
X_2 & Y_2 & Z_2 & 0 & \ddots& & \vdots \\
0 & \ddots & \ddots & \ddots & \ddots &  & \vdots \\
\vdots & &\ddots & \ddots & \ddots & \ddots  & \vdots \\
\vdots & & & \ddots & X_{I-1} & Y_{I-1}  & Z_{I-1} \\
0 & \cdots & \cdots & \cdots & 0 & X_I & Y_I\\
\end{bmatrix}\in\R^{I\times I}\label{eq:A}
\end{align}

By substituting the approximation ${v^{n}}'$ found in \cref{eq:dv-upwind} into equation \cref{eq:c-neoclassical-growth}, define the utility vector
\begin{align}
u(c^n) &\equiv u\left(\left({v^n}'\right)^{-\frac{1}{\gamma}}\right)\\
\intertext{We can then write as a system of equations in matrix form and solve for $v^{n+1}$}
B^n &\equiv \left(\rho + \frac{1}{h}\right)I - A^n\\
b^n &\equiv u(c^n) + \frac{v^n}{h}\\
B^{n}v^{n+1} &= {b^n}
\end{align}



When $v^{n+1}$ is sufficiently close in value to $v^n$, the algorithm is complete. Otherwise, update the value of $v^n$ and repeat the previous steps for the next iteration. Keep in mind that the solution may not be unconditionally stable for an arbitrary $h$.


\subsection{Nonlinear System of Equations}
This section will focus on solving the problem directly as a nonlinear system of equations. Returning to equation \cref{eq:HJBE-neoclassical-growth-mu}, we have 
\begin{align}
\rho v_i &= \frac{\left(v'_i\right)^{-\frac{1-\gamma}{\gamma}}}{1-\gamma} + \mu(k_i;v) v'_i
\end{align}
for every $k_i$, $i=1,2,...I$.

\bibliography{etk-references}

\end{document}