operator_discretization_finite_differences.tex

\documentclass[11pt]{etk-article}
\usepackage{pstool} 
\usepackage{etk-bib}
\pdfmetadata{}{}{}{}

\begin{document}
\title{Discretizing Continuous-Time Operators with Finite Differences}
\author{Jesse Perla\\UBC}
\date{\today}
\maketitle
 Here, we expand on details for how to discretize linear and nonlinear operators for use in other methods.\footnote{These notes expand on Ben Moll's superb notes in \url{http://www.princeton.edu/~moll/HACTproject/}. Also thanks to the excellent RA work by Sev Hou and Steven Zhang.}  Most of the notes will focus on time-homogeneity examples.\footnote{For more background on differential operators, see \url{https://en.wikipedia.org/wiki/Differential_operator}.}
\paragraph{General Notation} The following expressions are used in the document:
\begin{itemize}
	\item Given an operator (or infinitesimal generator) associated with a particular stochastic process, $\mathcal{A}$.\footnote{See \url{https://en.wikipedia.org/wiki/Infinitesimal_generator_(stochastic_processes)} for some formulas and interpretation for diffusions, and \url{https://en.wikipedia.org/wiki/Transition_rate_matrix} for the generator of continuous-time Markov chains.}  The purpose of these notes is to discretize $\mathcal{A}$ on this grid using finite differences.  Crucially, it is necessary to put boundary-values into the discretized operator as well. 
	\item Where appropriate, we will define the adjoint of the operator $\mathcal{A}$ as $\mathcal{A}^*$.  This is useful when trying to solve for the evolution of distributions (rather than solving for the value functions).
	\item For a given variable $q$, define the notation $q^{-} \equiv \min\set{q,0}$ and $q^{+} \equiv \max\set{q,0}$, which will be useful for defining finite-differences with an upwind scheme.\footnote{For more details on the notation and the upwind scheme, see \url{http://www.princeton.edu/~moll/HACTproject/HACT_Numerical_Appendix.pdf}.}   This can apply to vectors as well. For example, $q_i^{-} = q_i$ if $q_i < 0$ and $0$ if $q_i > 0$, and $q_i^{-} \equiv \set{q^{-}_i}_{i=1}^{I}$.
	\item Finally,  derivatives are denoted by the operator $\D$ and univariate derivatives such as $\D[x]v(x) \equiv v'(x)$.	
\end{itemize}

\paragraph{Grids} Start with the simplification of a univariate function of $x$,
\begin{itemize}
	\item In the univariate case with $I$ points, $\set{x_i}_{i=1}^I$ with $x_1 = \underline{x}$ and $x_I = \bar{x}$ when $x \in [\underline{x}, \bar{x}]$.  After discretizing, we will denote the grid with the variable name, i.e. $x \equiv \set{x_i}_{i=1}^I$
\item For the grid $x$, denote the forward and backward distance between grid points at node $i$ as,
\begin{align}
	\Delta_{i,+} &\equiv x_{i+1} - x_i\label{eq:Delta-i-plus}\\
	\Delta_{i,-} &\equiv x_i - x_{i-1}\label{eq:Delta-i-minus}	
\end{align}
In the simple case of a uniform grid, $\Delta \equiv \Delta_{i,+} = \Delta_{i,-} = x_{i+1}$ for all $i$
\item When we discreteize a function, use the function name without arguments to denote the vector.  i.e. $v(x)$ discretized on a grid $\set{x_i}_{i=1}^{I}$ is $v \equiv \set{v(x_i)}_{i=1}^I \in \R^I$
\end{itemize}

\section{Time-Independent Univariate Diffusions}\label{sec:time-independent-univariate-diffusion}
Assume that the time variable does not directly enter into anything important in the problem (e.g. drifts, payoffs, stopping values, boundary values, or discount rates).  This is also assuming that the drift, variance, and boundary values are not a (direct) function of any $v_i$.\footnote{However, in the case of controlled-diffusions a numerical method is to use the last iterations $v^{n-1}_i$ as part of the $\mu^n_i$ definition for the equations in $v^n_i$.  This process has a linear $\mathcal{A}$ at any given iteration.}  Otherwise, the operator would be nonlinear.  In the case of a linear $\mathcal{A}$, the resulting discretized operator will be a (heavily sparse) matrix $A \in \R^{I \times I}$, which is especially convenient to solve.

\subsection{Stochastic Process and Boundary Values}
Take a diffusion process for $x_t$ according to the following stochastic difference equation,
\begin{align}
\diff x_t = \mu(x_t)\diff t + \sigma(x_t)\diff\mathbb{W}_t\label{eq:x-SDE}
\end{align}
where $\mathbb{W}_t$ is Brownian motion and $x \in (\underline{x}, \bar{x})$ where $-\infinity \leq \underline{x} < \bar{x} \leq \infinity$.\footnote{We are being a little sloppy with notation $x_t$ being exactly at the bounds because it is a measure $0$ event in these setups.}  The functions $\mu(\cdot)$ and $\sigma(\cdot)$ are left general, but it is essential that they are time-homogeneous.

The partial differential operator (infinitesimal generator) associated with the stochastic process in \cref{eq:x-SDE} is\footnote{When time enters payoffs, stopping values, etc. the generator is instead on a function of both $t$ an $x$.  See \cref{sec:time-independent-univariate-diffusion}}
\begin{align}
	\mathcal{A} &\equiv \mu(x)\D[x] + \frac{\sigma(x)^2}{2}\D[xx]\label{eq:A-generator-univariate-diffusion}
\end{align}

We will consider several types of boundary values for an operator on a generic function $v(x)$.  Because the process is univariate, we require boundary conditions at both sides of $x$, and assume that the boundary and $\mu(\cdot), \sigma(x)^2$ are independent of time.  Given that the operator is on a function $v(x)$, consider:
\begin{itemize}
	\item \textbf{Absorbing Barrier:}   $v(\underline{x}) = \underline{v}$ and/or $v(\bar{x}) = \bar{v}$ for constants $\underline{v}$ and/or $\bar{v}$.  This is called a ``Dirichlet'' boundary condition in PDEs. 
	\item \textbf{Reflecting Barrier:} $\D[x]v(\underline{x}) = v'(\underline{x})=0$ and/or $v'(\bar{x}) = 0$.  This is called a ``Neumann'' boundary condition in PDEs. 
	\item \textbf{Constrained Process:}  In many cases, bounds on the variable are not essential to the process, but really just necessary to discretize.  In those cases, adding in an artificial ``reflecting barrier'' at some large $\bar{x}$ (or small $\bar{x}$) seems to be the best approach to constraining the process for numerical methods.\footnote{See \cite{KushnerDupuis2001} section 1.4, 5.7, and 11.1 for a discussion.}  The hope is that the introduction of this boundary condition would have little impact on the solution except very close to the artificial boundary.
	\item \textbf{Transversality:} Problem specific?  But the hope is that using a ``Constrained Process'' with a reflecting barrier converges to the correct solution in most of the domain.
\end{itemize}

To distinguish between the boundary values we will consider, define
\begin{align}
	R_{\underline{x}} &\equiv \begin{cases}
		1 & \text{ if $x$ is reflected at $\underline{x}$}\\
		0 & \text{ if there is a absorbing barrier with value $\underline{v}$}
	\end{cases}\label{eq:R-x-min}\\
	R_{\bar{x}} &\equiv \begin{cases}
	1 & \text{ if $x$ is reflected at $\bar{x}$}\\
	0 & \text{ if there is a absorbing barrier with value $\bar{v}$}
\end{cases}\label{eq:R-x-max}
\end{align}	
\paragraph{What Boundary Value to Assume?}
In general, a reflection for a large $\bar{x}$ and small $\underline{x}$ is the most appropriate.  The exceptions would tend to have direct economics associated with a death or uncontrolled stopping.\footnote{For example, even if you have firms optimally exiting, you could use a small reflection at $\bar{x}$ for the variational methods in \url{optimal_stopping.pdf} as long as $\bar{x}$ is below any binding level.}

\subsection{Upwind Finite Difference Discretization for $A$ with a Uniform Grid}
Denote the vector of drifts and variances as $\mu \equiv \set{\mu(x_i)}_{i=1}^I$ and $\sigma^2 \equiv \set{\sigma(x_i)^2}_{i=1}^I$.  For a uniform grid with $\Delta \equiv x_{i+1} - x_i$.  Define the vectors, $X,Y,Z \in \R^I$ such that,
\begin{align}
	X &= - \frac{\mu^{-}}{\Delta}+ \frac{\sigma^{2}}{2 \Delta^2}\label{eq:X} \\
	Y &= - \frac{\mu^{+}}{\Delta} + \frac{\mu^{-}}{\Delta^2}- \frac{\sigma^{2}}{\Delta^2}\label{eq:Y} \\
	Z &= \frac{\mu^{+}}{\Delta} + \frac{\sigma^{2}}{2 \Delta^2}\label{eq:Z}
\end{align}

With these, the matrix $A$ is constructed as.
\begin{align}
A &\equiv \begin{bmatrix}
Y_1 + R_{\underline{x}} X_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+R_{\bar{x}}Z_I\\
\end{bmatrix}\in\R^{I\times I}\label{eq:A}\\
b &\equiv \begin{bmatrix}(1-R_{\underline{x}})X_1\underline{v} \\ 0 \\ \vdots \\ 0 \\  (1-R_{\bar{x}})Z_I\bar{v}\end{bmatrix}^T\in\R^I \label{eq:b}
\intertext{where,}
\mathcal{A}v(\set{x_i}) &\approx A v + b
\end{align}

Note that the special cases for the $1$st and $I$th row correspond to the boundary values and are adjusted depending on the type of boundary value. Also, note that the $b$ vector is all zeros if the boundary values are reflective and/or the matched value $\bar{v}$ or $\underline{v}$ is $0$).

\paragraph{Intensity Matrix of a Continuous-Time Markov Chain}
Note that with reflecting barriers on both sides, the sum for each row of $A$ is $X_i +Y_i + Z_i = 0$ for all $i = 1, \ldots I$ from \cref{eq:X,eq:Y,eq:Z}, and it can be shown that the diagonal is always negative.  Hence, we interpret the $A$ matrix fulfills the requirements of an intensity matrix (or infinitesimal generator matrix) of a continuous-time Markov chain.\footnote{See \url{https://en.wikipedia.org/wiki/Transition_rate_matrix} and \url{https://en.wikipedia.org/wiki/Markov_chain}.}  Furthermore, the Markov-chain has the convenient property that jumps are only between adjacent states (i.e. if at $x_i$ then can only jump to $x_{i-1}$ and $x_{i+1}$), which is extremely sparse.

When the boundary value is absorbing at $\underline{x}$ or $\bar{x}$, the interpretation is a little trickier since the sum of all elements is not $0$.  Essentially, we need to take into account ``death'' since the agents are no longer reflected within the boundaries at all times.  This is the reason for a non-zero $b$ matrix.


\paragraph{Interior of $A$:}
To better understand the construction of $A$ and $b$, look at individual rows of $A$ with the ODE.

In the interior ($1 < i < I$), the discretization of \cref{eq:A-generator-univariate-diffusion}
\begin{align}
\mathcal{A} v(x_i) &\approx \underbrace{\left(v_i-v_{i-1}\right)\frac{\mu_i^{-}}{\Delta}+ \left(v_{i+1}-v_i\right)\frac{\mu_i^{+}}{\Delta}}_{\text{Upwind Scheme}}  + \frac{\sigma_i^2}{2} \frac{v_{i+1} - 2 v_i + v_{i-1}}{\Delta^2}\label{eq:A-generator-univariate-diffusion-interior}\\
\intertext{The upwind scheme chooses either forward or backward differences, depending on the sign of the drift.  Collecting terms, we see the derivation for the definitions in \cref{eq:X,eq:Y,eq:Z}}
&= \underbrace{\left(-\frac{\mu_i^{-}}{\Delta} +\frac{\sigma_i^2}{2 \Delta^2}\right)}_{\equiv X_i}v_{i-1} + \underbrace{\left(\frac{\mu_i^{-}}{\Delta} - \frac{\mu_i^{+}}{\Delta}-\frac{\sigma_i^2}{\Delta^2}\right)}_{\equiv Y_i}v_i + \underbrace{\left(\frac{\mu_i^{+}}{\Delta} + \frac{\sigma_i^2}{2 \Delta^2}\right)}_{\equiv Z_i}v_{i+1}\label{eq:A-collected-interior}
\end{align}


\paragraph{Boundary Value at $\underline{x}$:}
As $i =1$, given the ``ghost node'' $x_0$, the discretized operator from \cref{eq:A-collected-interior} is
\begin{align}
\mathcal{A} v(x_1) &\approx \left(-\frac{\mu_1^{-}}{\Delta} +\frac{\sigma_1^2}{2 \Delta^2}\right)v_0 + \left(\frac{\mu_1^{-}}{\Delta} - \frac{\mu_1^{+}}{\Delta}-\frac{\sigma_1^2}{\Delta^2}\right)v_1 + \left(\frac{\mu_1^{+}}{\Delta} + \frac{\sigma_1^2}{2 \Delta^2}\right)v_2\label{eq:A-collected-left}\\
\intertext{In the case of the boundary value $v(\underline{x}) = \underline{v}$, subsitute for $v_0 =  \underline{v}$ to find,}
&\approx \underbrace{\overbrace{\left(-\frac{\mu_1^{-}}{\Delta} +\frac{\sigma_1^2}{2 \Delta^2}\right)}^{\equiv X_1}\underline{v}}_{\equiv b_1} +  \underbrace{\left(\frac{\mu_1^{-}}{\Delta} - \frac{\mu_1^{+}}{\Delta}-\frac{\sigma_1^2}{\Delta^2}\right)}_{\equiv Y_1}v_1 + \underbrace{\left(\frac{\mu_1^{+}}{\Delta} + \frac{\sigma_1^2}{2 \Delta^2}\right)}_{\equiv Z_1}v_2\\
\intertext{Alternatively, if the boundary value is $v'(\underline{x}) = 0$, take \cref{eq:A-collected-left} and use $v'(\underline{x}) \approx \frac{v_1-v_0}{\Delta} = 0, \implies v_1 = v_0$, so}
\mathcal{A} v(x_1) &\approx \left(\underbrace{\left(-\frac{\mu_1^{-}}{\Delta} +\frac{\sigma_1^2}{2 \Delta^2}\right)}_{\equiv X_1} + \underbrace{\left(\frac{\mu_1^{-}}{\Delta} - \frac{\mu_1^{+}}{\Delta}-\frac{\sigma_1^2}{\Delta^2}\right)}_{\equiv Y_1}\right)v_1 + \underbrace{\left(\frac{\mu_1^{+}}{\Delta} + \frac{\sigma_1^2}{2 \Delta^2}\right)}_{\equiv Z_1}v_2
\end{align}

\paragraph{Boundary Value at $\bar{x}$:}
As $i=I$, given the ``ghost node'' $x_{I+1}$, from \cref{eq:A-collected-interior}
\begin{align}
\mathcal{A} v(x_I)&\approx \left(-\frac{\mu_I^{-}}{\Delta} +\frac{\sigma_I^2}{2 \Delta^2}\right)v_{I-1} + \left(\frac{\mu_I^{-}}{\Delta} - \frac{\mu_I^{+}}{\Delta}-\frac{\sigma_I^2}{\Delta^2}\right)v_I + \left(\frac{\mu_I^{+}}{\Delta} + \frac{\sigma_I^2}{2 \Delta^2}\right)v_{I+1}\\
\intertext{For the absorbing barrier, substitute for $v(x_{I+1}) = \bar{v}$,}
\mathcal{A} v(x_I)&\approx \underbrace{\left(-\frac{\mu_I^{-}}{\Delta} +\frac{\sigma_I^2}{2 \Delta^2}\right)}_{\equiv X_I} v_{I-1} + \underbrace{\left(\frac{\mu_I^{-}}{\Delta} - \frac{\mu_I^{+}}{\Delta}-\frac{\sigma_I^2}{\Delta^2}\right)}_{\equiv Y_I} v_I + \underbrace{\overbrace{\left(\frac{\mu_I^{+}}{\Delta} + \frac{\sigma_I^2}{2 \Delta^2}\right)}^{\equiv Z_I}\bar{v}}_{\equiv b_I}\\
\intertext{For a reflecting barrier, the boundary value $v'(\bar{x}) \approx \frac{v_{I+1}-v_I}{\Delta} = 0, \implies v_{I+1} = v_I$,}
\mathcal{A} v(x_I)&\approx \underbrace{\left(-\frac{\mu_I^{-}}{\Delta} +\frac{\sigma_I^2}{2 \Delta^2}\right)}_{\equiv X_I} v_{I-1} + \left(\underbrace{\left(\frac{\mu_I^{-}}{\Delta} - \frac{\mu_I^{+}}{\Delta}-\frac{\sigma_I^2}{\Delta^2}\right)}_{\equiv Y_I} + \underbrace{\left(\frac{\mu_I^{+}}{\Delta} + \frac{\sigma_I^2}{2 \Delta^2}\right)}_{\equiv Z_I} \right) v_I 
\end{align}
\paragraph{Simple Case with $\mu(x) < 0$ and Reflecting Barrier}
In the special case of $\mu(\bar{x}) < 0$, the upwind direction is trivial: $\mu^{-}= \mu$ and $\mu^{+} = 0$.  Specializing \cref{eq:X,eq:Y,eq:Z},
\begin{align}
	X &= - \frac{\mu}{\Delta}+ \frac{\sigma^{2}}{2 \Delta^2}\label{eq:X-backwards} \\
	Y &= \frac{\mu}{\Delta} - \frac{\sigma^{2}}{\Delta^2}\label{eq:Y-backwards} \\
	Z &= \frac{\sigma^{2}}{2 \Delta^2}\label{eq:Z-backwards}
\end{align}
and the formula for the discretized operator is,
\begin{align}
	\mathcal{A} v(x_i)  &= \left(v_i-v_{i-1}\right)\frac{\mu}{\Delta} + \frac{\sigma^2}{2} \frac{\left(v_{i+1} - 2 v_i + v_{i-1}\right)}{\Delta^2},\quad \text{ for } i = 2,\ldots I-1\\
	\intertext{With the reflecting barrier at $\underline{x}$ and $\bar{x}$,}
	\mathcal{A} v(x_1) &= -\frac{\sigma^{2}}{2 \Delta^2}v_1 + \frac{\sigma^{2}}{2 \Delta^2} v_2\\
\mathcal{A} v(x_I)&\approx \left(-\frac{\mu}{\Delta} +\frac{\sigma^2}{2 \Delta^2}\right) v_{I-1} + \left(\frac{\mu}{\Delta} - \frac{\sigma^2}{2 \Delta^2}\right) v_I 
\end{align}




\subsection{Upwind Finite Difference Discretization for $A$ with a Non-Uniform Grid}
For non-uniform grid, $\Delta_{+} \equiv x_{i+1} - x_{i}$ and $\Delta_{-} \equiv x_{i} - x_{i-1}$\\
Based on the recommendation of second derivative choice from \url{http://www.princeton.edu/~moll/HACTproject/HACT_Numerical_Appendix.pdf}, a good candidate approximation is
\begin{equation}
v_i^{''} = \frac{\Delta_{-}v_{i+1} - (\Delta_{-}+\Delta_{+})v_i+\Delta_{+}v_{i-1}}{\frac{1}{2}(\Delta_{+}+\Delta_{-})\Delta_{-}\Delta_{+}} \label{eq:v''}
\end{equation}
Then we can define $X$, $Y$ and $Z$ (after simplification) as 
\begin{align}
X% &= -\frac{\mu^{-}}{\Delta_{-}} + \frac{\sigma^2}{2}\frac{\frac{\Delta_{+}}{\frac{1}{2}(\Delta_{+}+\Delta_{-})}}{\Delta_{+}\Delta_{-}}\\
&= -\frac{\mu^{-}}{\Delta_{-}} +\frac{\sigma^{2}}{\Delta_{-}(\Delta_{+}+\Delta_{-})}\label{eq:X-non-uniform}\\
Y &= -\frac{\mu^{+}}{\Delta_{+}}+\frac{\mu^{-}}{\Delta_{-}} -\frac{\sigma^2}{\Delta_{+}\Delta_{-}}\label{eq:Y-non-uniform}\\
Z %&= \frac{\mu^{+}}{\Delta_{+}} + \frac{\sigma^2}{2}\frac{\frac{\Delta_{-}}{\frac{1}{2}(\Delta_{+}+\Delta_{-})}}{\Delta_{+}\Delta_{-}}\\
&= \frac{\mu^{+}}{\Delta_{+}} + \frac{\sigma^{2}}{\Delta_{+}(\Delta_{+}+\Delta_{-})}\label{eq:Z-non-uniform}
\end{align}
%Turned off.  More trouble than worth
%For the stability of numerical calculation, we can multiply every term in $A$ by $\Delta_{-}$ to obtain $\Delta_{-}A$. So we can redefine $X, Y $and $Z$ as 
%\begin{align}
%X &= -\mu^{-} + \frac{\sigma^2}{(\Delta_{+}+\Delta_{-})} \label{eq:X-non-uniform}\\
%Y &= -\mu^{+}\frac{\Delta_{-}}{\Delta_{+}}+\mu^{-} -\frac{\sigma^2}{\Delta_{+}} \label{eq:Y-non-uniform}\\
%Z &=\mu^{+}\frac{\Delta_{-}}{\Delta_{+}} + \frac{\Delta_{-}}{\Delta_{+}}\frac{\sigma^2}{(\Delta_{+}+\Delta_{-})} \label{eq:Z-non-uniform}
%\end{align}
We treat $\Delta_{1, -} = \Delta_{1, +}$ and $\Delta_{I, +} = \Delta_{I, -}$, due to ghost points, $x_0$ and $x_{I+1}$ on both boundaries.  \\
%This isn't really well define in terms of matrix multiplication/etc.
%In terms of the structure of \ref{eq:A} and \ref{eq:b}, we can similarly write 
%\begin{equation}
%\mathcal{A}v(\set{x_i}) \approx \frac{A}{\Delta_{-}} v + \frac{b}{\Delta_{-}}
%\end{equation}


\paragraph{Calculating in Matrix Form for Reflecting Barriers}
First, calculate the first difference vectors $\Delta_{-}, \Delta_{+} \in \R^I$.  From the $x \in \R^I$ vector, this can be calculated in using the first-difference function in matlab/etc. as the identical structure to \cref{eq:A}
\begin{align}
\Delta_{-} &\equiv \begin{bmatrix} x_2 - x_1 \\
\text{diff}(x)
\end{bmatrix}\\
\Delta_{+} &\equiv \begin{bmatrix} \text{diff}(x)\\
x_I - x_{I-1}
\end{bmatrix}
\end{align}
Next, with this definition \cref{eq:X-non-uniform,eq:Y-non-uniform,eq:Z-non-uniform} hold directly as vectors.  With these, the matrix $A$ is constructed (assuming reflecting boundaries) as.
\begin{align}
A &\equiv \begin{bmatrix}
Y_1 + X_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+Z_I\\
\end{bmatrix}\in\R^{I\times I}\label{eq:A-non-uniform}
\end{align}
In terms of the structure of \cref{eq:A,eq:b}, we can similarly write
%\footnote{Note that $\mathbf{I}\Delta_{-}$ creates diagonal of $\Delta_{i,-}$}
\begin{align}
%\mathbf{I}\Delta_{-} \mathcal{A}v(\set{x_i}) &\approx A v\\
%\intertext{However, this can be more succintly written with the definition,}
%\mathbf{\Delta}_{-} &\equiv \mathbf{I}\Delta_{-} = \text{diag}(\Delta_{-})\\
%\mathbf{\Delta}_{+} &\equiv \mathbf{I}\Delta_{+} = \text{diag}(\Delta_{+})\\
%\intertext{With this and the above,}
%\mathbf{\Delta}_{-} \mathcal{A}v(\set{x_i}) &\approx A v
%\intertext{Or,}
\mathcal{A}v(\set{x_i}) &\approx  A v
\end{align}




\paragraph{Interior of $A$:}
To better understand the construction of $A$ and $b$, look at individual rows of $A$ with the ODE.

In the interior ($1 < i < I$), the discretization of \cref{eq:A-generator-univariate-diffusion}
\begin{align}
\mathcal{A} v(x_i) &\approx \underbrace{\left(v_i-v_{i-1}\right)\frac{\mu_i^{-}}{\Delta_{i,-}}+ \left(v_{i+1}-v_i\right)\frac{\mu_i^{+}}{\Delta_{i,+}}}_{\text{Upwind Scheme}}  + \sigma_i^2 \frac{\Delta_{i,-}v_{i+1} - (\Delta_{i,+}+\Delta_{i,-}) v_i + \Delta_{i,+}v_{i-1}}{\Delta_{i,-}\Delta_{i,+}(\Delta_{i,+}+\Delta_{i,-})}\label{eq:A-generator-univariate-diffusion-with-nonuniform-grid-interior}\\
\intertext{The upwind scheme chooses either forward or backward differences, depending on the sign of the drift.  Collecting terms, we see the derivation for the definitions in \cref{eq:X-non-uniform,eq:Y-non-uniform,eq:Z-non-uniform}}
&= \underbrace{\left(-\frac{\mu_i^{-}}{\Delta_{i,-}} +\frac{\sigma_i^2}{\Delta_{i,-}(\Delta_{i,+}+\Delta_{i,-})}\right)}_{\equiv X_i}v_{i-1} + \underbrace{\left(\frac{\mu_i^{-}}{\Delta_{i,-}} - \frac{\mu_i^{+}}{\Delta_{i,+}}-\frac{\sigma_i^2}{\Delta_{i,-}\Delta_{i,+}}\right)}_{\equiv Y_i}v_i + \underbrace{\left(\frac{\mu_i^{+}}{\Delta_{i,+}} + \frac{\sigma_i^2}{\Delta_{i,+}(\Delta_{i,+}+\Delta_{i,-})}\right)}_{\equiv Z_i}v_{i+1}\label{eq:A-collected-interior-with-nonuniform-grid}
\end{align}

\paragraph{Boundary Value at $\underline{x}$:}
As $i =1$, given the ``ghost node'' $x_0$, we recall the assumption that $\Delta_{1,-} = \Delta_{1,+}$.  Hence, the discretized operator from \cref{eq:A-collected-interior-with-nonuniform-grid} is
\begin{align}
\mathcal{A} v(x_1) = &\approx \left(-\frac{\mu_1^{-}}{\Delta_{1,+}} +\frac{\sigma_1^2}{2 \Delta_{1,+}^2}\right)v_0 + \left(\frac{\mu_1^{-} - \mu_1^{+}}{\Delta_{1,+}}-\frac{\sigma_1^2}{\Delta_{1,+}^2}\right) v_1 + \left(\frac{\mu_1^{+}}{\Delta_{1,+}} + \frac{\sigma_1^2}{2 \Delta_{1,+}^2}\right)v_2\label{eq:A-collected-with-nonuniform-grid-left}\\
\intertext{If the boundary value is $v'(\underline{x}) = 0$, take \cref{eq:A-collected-with-nonuniform-grid-left} and use $v'(\underline{x}) \approx \frac{v_1-v_0}{\Delta} = 0, \implies v_1 = v_0$, so}
 &\approx \left(\underbrace{\left(\frac{-\mu_1^{-}}{\Delta_{1,-}} +\frac{\sigma_1^2}{\Delta_{1,-}(\Delta_{1,+}+\Delta_{1,-})}\right)}_{\equiv X_1} + \underbrace{\left(\mu_1^{-} - \frac{\mu_1^{+}}{\Delta_{1,+}}-\frac{\sigma_1^2}{\Delta_{1,+}\Delta_{1,-}}\right)}_{\equiv Y_1}\right)v_1 + \underbrace{\left(\frac{\mu_1^{+}}{\Delta_{1,+}} + \frac{\sigma_1^2}{2 \Delta_{1,+}\Delta_{1,-}}\right)}_{\equiv Z_1}v_2\label{eq:reflecting-left-nonuniform}
\end{align}

\paragraph{Boundary Value at $\bar{x}$:}
As $i=I$, given the ``ghost node'' $x_{I+1}$, and the assumption  $\Delta_{I,+} = \Delta_{I,-}$, from \cref{eq:A-collected-interior},
\begin{align}
\mathcal{A}  v(x_I)&\approx \left(-\frac{\mu_I^{-}}{\Delta_{I,-}} +\frac{\sigma_I^2}{2 \Delta_{I,-}^2}\right)v_{I-1} + \left(\frac{\mu_I^{-} - \mu_I^{+}}{\Delta_{I,-}}-\frac{\sigma_I^2}{\Delta_{I,-}^2}\right)v_I + \left(\frac{\mu_I^{+}}{\Delta_{I,-}} + \frac{\sigma_I^2}{2 \Delta_{I,-}^2}\right)v_{I+1}\\
\intertext{For a reflecting barrier, the boundary value $v'(\bar{x}) \approx \frac{v_{I+1}-v_I}{\Delta} = 0, \implies v_{I+1} = v_I$,}
&\approx \underbrace{\left(-\frac{\mu_I^{-}}{\Delta_{I,-}} +\frac{\sigma_I^2}{\Delta_{I,-}(\Delta_{I,+}+\Delta_{I,-})}\right)}_{\equiv X_I} v_{I-1}\nonumber\\
&+ \left(\underbrace{\left(\frac{\mu_I^{-}}{\Delta_{I,-}} - \frac{\mu_I^{+}}{\Delta_{I,+}}-\frac{\sigma_I^2}{\Delta_{I,-}\Delta_{I,+}}\right)}_{\equiv Y_I} + \underbrace{\left(\frac{\mu_I^{+}}{\Delta_{I,+}} + \frac{\sigma_I^2}{\Delta_{I,+}(\Delta_{I,+}+\Delta_{I,-})}\right)}_{\equiv Z_I} \right) v_I \label{eq:reflecting-right-nonuniform}
\end{align}
\paragraph{Simple Case with $\mu(x) < 0$ and Reflecting Barrier}
In the special case of $\mu(\bar{x}) < 0$, the upwind direction is trivial: $\mu^{-} = \mu$ and $\mu^{+} = 0$.  Specializing \cref{eq:X-non-uniform,eq:Y-non-uniform,eq:Z-non-uniform,eq:reflecting-left-nonuniform,eq:reflecting-right-nonuniform}
\begin{align}
	X_i &= - \frac{\mu_i}{\Delta_{i,-}}+ \frac{\sigma_i^{2}}{\Delta_{i,-}(\Delta_{i,+}+\Delta_{i,-})}\label{eq:X-non-uniform-backwards} \\
	Y_i &= \frac{\mu_i}{\Delta_{i,-}} - \frac{\sigma_i^{2}}{\Delta_{i,+}\Delta_{i,-}}\label{eq:Y-non-uniform-backwards} \\
	Z_i &= \frac{\sigma_i^{2}}{\Delta_{i,+}(\Delta_{i,+}+\Delta_{i,-})}\label{eq:Z-non-uniform-backwards}
\end{align}
%and the formula for the discretized operator is,
%\begin{align}
%	\Delta_{-} \mathcal{A} v(x_i)  &= \left(v_i-v_{i-1}\right)\mu + \sigma^2 \frac{\Delta_{-}v_{i+1} - (\Delta_{+}+\Delta_{-}) v_i + \Delta_{+}v_{i-1}}{\Delta_{+}(\Delta_{+}+\Delta_{-})},\quad \text{ for } i = 2,\ldots I-1\\
%	\intertext{With the reflecting barrier at $\underline{x}$ and $\bar{x}$,}
%	\Delta_{-} \mathcal{A} v(x_1) &= -\frac{\Delta_{-}}{\Delta_{+}}\frac{\sigma^{2}}{\Delta_{+}+\Delta_{-}}v_1 + \frac{\Delta_{-}}{\Delta_{+}}\frac{\sigma^{2}}{\Delta_{+}+\Delta_{-}} v_2\\
%\Delta_{-} \mathcal{A} v(x_I)&\approx -\left(\mu -\frac{\sigma^2}{\Delta_{-}+\Delta_{+}}\right) v_{I-1} + \left(\mu - \frac{\sigma^2}{\Delta_{+}+\Delta_{-}}\right) v_I 
%\end{align}
Where $\Delta_{1,-} = \Delta_{1,+}$ and $\Delta_{I,+} = \Delta_{I,-}$ as required.

\subsection{Stationary Distribution of Reflected Univariate Diffusions}
Given a model without any birth/death/stopping, let the pdf of the distribution of $x$ be,
\begin{align}
	\D[t]f(t,x) &= \mathcal{A}^* f(t,x),\, \text{with appropriate BVs and }\label{eq:inhomogenous-kfe-univariate-diffusion}\\
	1 &= \int_0^1 f(t,x)\diff x,\label{eq:inhomogenous-kfe-univariate-diffusion-integral} \, \text{ for all $t$}
	\intertext{where $\mathcal{A}^*$ is the adjoint operator of the $\mathcal{A}$ operator of the stochastic process.  In the stationary set, $\D[t]f(t,x) = 0$ for all $x$.  That is,}
	0 &= \mathcal{A}^*f(x)\\
	1 &= \int_0^1 f(x) \diff x
	\intertext{If the linear operator $\mathcal{A}$ is discretized as a matrix $A$, then it can be shown that the discretized $\mathcal{A}^*$ is the transpose $A^T$.  To find the discretized PDF of the distribution, $f \in \R^I$, solve}
	0 &= A^T f\label{eq:linear-operator-stationary}\\
	1 &= \omega \cdot f\label{eq:linear-operator-sum}
\end{align}
where $\omega \in \R^I$ is a quadrature weight to ensure it integrates to $1$.  For crude simplicity, could use a vector of $1$s.  Alternatively, $\omega$ could be trapezoidal or simpsons rule weights (depending on if the $x$ grid is uniform).
\paragraph{As an Eigenvalue Problem} Note that \cref{eq:linear-operator-stationary} is an eigenvalue problem.  Therefore, we can find the stationary distribution can be found as the eigenvector of $A^T$ associated with the zero eigenvalue, and rescales to ensure \cref{eq:linear-operator-sum} holds.  This is useful since $A^T$ is always singular, and hence \cref{eq:linear-operator-stationary} cannot be solved as a simple linear system.

\paragraph{As a Linear System} Alternatively, define the following, after multiplying \cref{eq:linear-operator-stationary,eq:linear-operator-sum} by $\Delta$ and stacking,
\begin{align}
y &\equiv \begin{bmatrix}0 \\ \vdots \\ 0 \\ 1\end{bmatrix}\in\R^{I + 1}\\
X &\equiv \begin{bmatrix}
	A^T\\
	\omega
\end{bmatrix}\in \R^{(I+1)\times I}
\intertext{Then, the stationary distribution, $f\in\R^I$, solves the following linear system}
X f &= y\label{eq:linear-system}
\end{align}
Note that this is an overdetermined system, but the rank of $X$ is $I$ rather than $I+1$.  Hence, direct solvers for non-square systems should solve it.  However, in practice when the system is large and sparse, it may be best solved as a linear least squares problem---with the understanding that the residual should be asymptotically zero.\footnote{Consider adding a $f \geq 0$ inequality constraint if necessary}

\begin{align}
	\min_f& \norm{y - X f}^2\label{eq:KFE-least-squares}
\end{align}

\subsection{Solving Simple Time-Independent Linear HJBE}\label{sec:simple-HJBE}
The following works for either a uniform or nonuniform grid.

In the simple case of an uncontrolled linear stochastic process with no optimal stopping, if an agent has payoff $u(x)$ with discounting $\rho > 0$, then the following HJBE determines the value
\begin{align}
	\rho v(x) &= u(x) + \mathcal{A}v(x)\\
	\intertext{subject to the appropriate boundary conditions.  Discretizing the grid, let $u \equiv \set{u(x_i)}_{i=1}^I\in\R^I$, then the discretized version is,}
	\rho v &= u + A v\\
	\intertext{Or,}
	(\rho  \mathbf{I} - A) v &= u\label{eq:discretized-linear-HJBE}
	\intertext{This simple linear system can be solved for $v$.  In principle, the HJBE and the KFE could be solved at the same time (here they are uncoupled, but in principle there could be direct connections in a nonlinear setup) by stacking $q \equiv \begin{bmatrix}v & f\end{bmatrix}^T\in\R^{2 I}$ \cref{eq:discretized-linear-HJBE,eq:KFE-least-squares}.  Define,}
y &\equiv \begin{bmatrix}u \\0 \\ \vdots \\ 0 \\ 1\end{bmatrix}\in\R^{2 I + 1}\\
X &\equiv \begin{bmatrix}
	\rho \mathbf{I} - A & \mathbf{0}_{I\times I}\\
	\mathbf{0}_{I\times I} & A^T\\
	\mathbf{0}_{1\times I} & \omega
\end{bmatrix}\in \R^{(2 I+1)\times 2 I}	
\intertext{Then, the joint solution to the problems is a $q$,}
X q &= y
	 \intertext{As above, we would expect $X$ to have only rank $2 I$, and hence this is an overdetermined problem with a unique solution.  Alternatively we could solve it as a linear least squares problem,}
	\min_{q}& \norm{y - X q}^2\label{eq:KFE-HJBE-least-squares}	
\end{align}
%
%\subsection{Solving Simple Time-Independent Linear HJBE with a Nonuniform Grid}\label{sec:simple-HJBE-nonuniform}
%To derive the equivalent to \cref{eq:discretized-linear-HJBE}
%\begin{align}
%\rho v &= u + A v\\
%\text{or,}
%( \rho \mathbf{\Delta}_{-} - A) v &= \mathbf{\Delta}_{-} u\label{eq:discretized-linear-HJBE-non-uniform}
%\end{align}
%The equation in \cref{eq:discretized-linear-HJBE-non-uniform} can be solved for $v$ on the non-uniform grid.

\section{Time-Dependent Univariate Diffusions}\label{sec:time-dependent-univariate-diffusion}
This section implements a version of \cref{sec:time-dependent-univariate-diffusion}, but where important parameters/settings/etc. vary with time.  Here, assume $t \in [0, T]$, where $T$ may be a very large number if the system converges.  The SDE for the state is now \begin{align}
	\diff x_t = \mu(t, x_t)\diff t + \sigma(t, x_t)\diff\mathbb{W}_t\label{eq:x-t-SDE}
\end{align}
and the infinitesimal generator is 
\begin{align}
	\mathcal{A} &\equiv \mu(t,x)\D[x] + \frac{\sigma(t,x)^2}{2}\D[xx] + \D[t]\label{eq:A-generator-univariate-diffusion-t}
\end{align}
Even if the $\mu(\cdot)$ and $\sigma(\cdot)$ are independent of time, the discretized generator used for solving the problem may be dependent on time due to boundary values or changes outside of the matrix itself.


\paragraph{Boundary Values}
As before, we will start by assuming reflecting boundaries at some $\underline{x}$ and $\bar{x}$, so that the boundary values are,
\begin{align}
\D[x]v(t,\underline{x}) &= 0\\
\D[x]v(t,\bar{x}) &= 0
\end{align}

Assume that for all $t \geq T$ the process (and payoffs) all converge.  i.e. $\mu(t,x) = \mu(T,x)$ for all $t \geq T$, with similar convergence for $\sigma(t,x)$ and payoffs $u(t,x)$.  Hence, for the boundary value in the time dimension, we could use the solution to the stationary problem.  For example, lets say that $\bar{v}(x)$ is the solution to stationary setup, then the boundary value is
\begin{align}
v(T,x) &= \bar{v}(x)
\intertext{However, the more natural approach is to simply say that the equation is in a steady state (i.e., the time derivative is $0$)}
\D[t]v(t,x) &= 0,\, \text{ for all } t \geq T\label{eq:D-steady-state-T}
\end{align} 
The second approach is ideal because we don't have to solve for the stationary solution separately.   Moreover, the strong suspicion is that using the stationary boundary value on the derivative will essentially solve for the stationary equilibrium as part of the system of equations.  Basically, if you look at the bottom right corner of the stacked system of equations after imposing the boundary value, it should be exactly the system we solved to find the stationary setup.
\subsection{Upwind Discretization with Time Variation}
The grid on $x \in [\underline{x},\bar{x}]$ is as discussed before, with step-sizes $\Delta_{+}$ and $\Delta_{-}$ (and $\Delta$ if uniform grid).  The grid on time $t \in [0,T]$ has step sizes $h_{+}$ and $h_{-}$ (and just $h$ if uniform). The gridpoints are index $n=1,\ldots N$ for time, where $t_1 = 0$ and $t_N = T$.  That is,
\begin{align}
h_{n,+} &\equiv t_{n+1} - t_n\label{eq:h-i-plus}\\
h_{n,-} &\equiv t_n - t_{n-1}\label{eq:h-i-minus}	
\end{align}  
\paragraph{Discretizing Drift, Volatility, and Payoff}
As $\mu(t,x), \sigma(t,x),$ and any payoff $u(t,x)$ may now be time varying, we can discretize them as,
\begin{align}
	\mu^n_i &\equiv \mu(t_n, x_i)\\
	\sigma^n_i &\equiv \sigma(t_n, x_i)\\
	u^n_i &\equiv u(t_n, x_i)\\
	\intertext{Stack these up to form a vector for a given time period.  For example,}
	\mu^n &\equiv \begin{bmatrix} \mu^n_1 &  \mu^n_2 & \ldots & \mu^n_I\end{bmatrix}^{\intercal}\in \R^I \label{discretize}
	\intertext{Then stack these up over time,}
	\mu &\equiv \begin{bmatrix} \mu^{1\intercal} &  \mu^{2\intercal} & \ldots&  \mu^{N\intercal}\end{bmatrix}^{\intercal}\in \R^{N I}
\end{align}
Define $\sigma^n$, $\sigma$, etc. similarly.

Following \cref{eq:X-non-uniform,eq:Y-non-uniform,eq:Z-non-uniform} we can define for each time step $X^n$, $Y^n$ and $Z^n$ (after simplification) as 
\begin{align}
X^n &\equiv -\frac{\mu^{n,-}}{\Delta_{-}} +\frac{\left(\sigma^n\right)^2}{\Delta_{-}(\Delta_{+}+\Delta_{-})}\label{eq:X-non-uniform-t}\\
Y^n &\equiv -\frac{\mu^{n,+}}{\Delta_{+}}+\frac{\mu^{n,-}}{\Delta_{-}} -\frac{\left(\sigma^n\right)^2}{\Delta_{+}\Delta_{-}}\label{eq:Y-non-uniform-t}\\
Z^n &\equiv \frac{\mu^{n,+}}{\Delta_{+}} + \frac{\left(\sigma^n\right)^2}{\Delta_{+}(\Delta_{+}+\Delta_{-})}\label{eq:Z-non-uniform-t}
\end{align}

Next, make a version of \cref{eq:A} to make a block for a given time,
\begin{align}
A^n &\equiv \begin{bmatrix}
Y^n_1 + X^n_1 & Z^n_1 & 0 & \cdots & \cdots & \cdots & 0 \\
X^n_2 & Y^n_2 & Z^n_2 & 0 & \ddots& & \vdots \\
0 & \ddots & \ddots & \ddots & \ddots &  & \vdots \\
\vdots & &\ddots & \ddots & \ddots & \ddots  & \vdots \\
\vdots & & & \ddots & X^n_{I-1} & Y^n_{I-1}  & Z^n_{I-1} \\
0 & \cdots & \cdots & \cdots & 0 & X^n_I & Y^n_I+Z^n_I\\
\end{bmatrix}\in\R^{I\times I}\label{eq:A-t}
\end{align}

Define the following diagonal matrices from the time-stepping for forward vs. backward differences,
\begin{align}
D_{-}^n &\equiv \frac{1}{h_{n,-}} \mathbf{I}\\
D_{+}^n &\equiv \frac{1}{h_{n,+}} \mathbf{I}
\end{align}

\paragraph{Implicit vs. Explicit Time Stepping}
The key choice when discretizing the time dimension with a PDE is whether to use a derivative that is implicit vs. explicit in time.\footnote{The tradeoff, typically, is that explicit methods are faster (because you can step forward in time given a solution to $v_n^i$) but they can be unstable (i.e., choosing the wrong $\Delta/h$ ratio and it diverges, even if both are small!).  Implicit methods are much more stable and often have unconditional convergence, but require solving systems of equations.
	
	For economic applications, the speed advantage of explicit methods may be less useful because this setup is rarely an initial value problem (i.e. $v^1_i$ is not given).  On the other hand, the boundary/initial condition at $t=0$ is unclear with when using an implicit method.}
\begin{itemize}
	\item Explicit: $\D[t]v(t_n,x_i) \approx \frac{v^{n+1}_i - v^n_i}{h_{n,+}}$
	\item Implicit: $\D[t]v(t_n,x_i) \approx \frac{v^n_i - v^{n-1}_i}{h_{n,-}}$
\end{itemize}


\subsection{Algebra for Explicit Time Stepping}
If we use explicit time stepping, then the discretization for uniform and non-uniform grids is a modification of  \cref{eq:A-generator-univariate-diffusion-with-nonuniform-grid-interior,eq:A-generator-univariate-diffusion-interior} in the interior $0 < n < N$ and $1 < i < I$:

\begin{align}
\mathcal{A} v(t_n, x_i) &\approx \underbrace{\left(v^n_i-v^n_{i-1}\right)\frac{\mu_i^{n,-}}{\Delta}+ \left(v^n_{i+1}-v^n_i\right)\frac{\mu_i^{n,+}}{\Delta}}_{\text{Upwind Scheme}}  + \frac{(\sigma^n_i)^2}{2} \frac{v^n_{i+1} - 2 v^n_i + v^n_{i-1}}{\Delta^2} + \underbrace{\frac{v^{n+1}_i - v^n_i}{h}}_{\text{Explicit}} \label{eq:A-generator-univariate-diffusion-interior-t}\\
\intertext{And with a non-uniform grid}
\mathcal{A} v(t_n, x_i) &\approx \underbrace{\left(v^n_i-v^n_{i-1}\right)\frac{\mu_i^{n, -}}{\Delta_{i,-}}+ \left(v^n_{i+1}-v^n_i\right)\frac{\mu_i^{n, +}}{\Delta_{i,+}}}_{\text{Upwind Scheme}}  + (\sigma^n_i)^2 \frac{\Delta_{i,-}v^n_{i+1} - (\Delta_{i,+}+\Delta_{i,-}) v^n_i + \Delta_{i,+}v^n_{i-1}}{\Delta_{i,-}\Delta_{i,+}(\Delta_{i,+}+\Delta_{i,-})}+ \underbrace{\frac{v^{n+1}_i - v^n_i}{h_{n,+}}}_{\text{Explicit}}\label{eq:A-generator-univariate-diffusion-with-nonuniform-grid-interior-t}
\end{align}

\paragraph{Interior of $A$:}
In the interior ($1 < i < I$), the discretization of \cref{eq:A-generator-univariate-diffusion-interior-t} and \cref{eq:A-generator-univariate-diffusion-with-nonuniform-grid-interior-t}
\begin{align}
\mathcal{A} v(t_n,x_i) & = \underbrace{\left(-\frac{\mu_i^{n,-}}{\Delta} +\frac{(\sigma_i^n)^2}{2 \Delta^2}\right)}_{\equiv X_i}v_{i-1}^n + \underbrace{\left(\frac{\mu_i^{n,-}-\mu_i^{n,+}}{\Delta}-\frac{(\sigma_i^n)^2}{\Delta^2}\right)}_{\equiv Y_i}v_i^n + \underbrace{\left(\frac{\mu_i^{n,+}}{\Delta} + \frac{(\sigma_i^n)^2}{2 \Delta^2}\right)}_{\equiv Z_i}v_{i+1}^{n}-\frac{1}{h}v_{i}^{n}+\frac{1}{h}v_{i}^{n+1}\label{eq:A-collected-interior-with-uniform-grid-t}
\intertext{And with a non-uniform grid}
\mathcal{A} v(t_n,x_i) & = \underbrace{\left(-\frac{\mu_i^{n,-}}{\Delta_{i,-}} +\frac{(\sigma_i^n)^2}{\Delta_{i,-}(\Delta_{i,+}+\Delta_{i,-})}\right)}_{\equiv X_i}v_{i-1}^n + \underbrace{\left(\frac{\mu_i^{n,-}}{\Delta_{i,-}}-\frac{\mu_i^{n,+}}{\Delta_{i,+}}-\frac{(\sigma_i^n)^2}{\Delta_{i,-}\Delta_{i,+}}\right)}_{\equiv Y_i}v_i^n + \\
&\underbrace{\left(\frac{\mu_i^{n,+}}{\Delta_{i,+}}+ \frac{(\sigma_i^n)^2}{\Delta_{i,+}(\Delta_{i,+}+\Delta_{i,-})}\right)}_{\equiv Z_i}v_{i+1}^n-\frac{1}{h_{n,+}}v_{i}^{n}+\frac{1}{h_{n,+}}v_{i}^{n+1}\label{eq:A-collected-interior-with-nonuniform-grid-t}
\end{align}


\paragraph{Boundary Value at $\underline{x}$:}
As $i =1$, given the ``ghost node'' $x_0$, the discretized operator from \cref{eq:A-collected-interior-with-uniform-grid-t} is
\begin{align}
\mathcal{A} v(t_n,x_1) &\approx \left(-\frac{\mu_1^{n,-}}{\Delta} +\frac{(\sigma_1^n)^2}{2 \Delta^2}\right)v_0^n + \left(\frac{\mu_1^{n,-}-\mu_1^{n,+}}{\Delta}-\frac{(\sigma_1^n)^2}{\Delta^2}\right)v_1^n + \left(\frac{\mu_1^{n,+}}{\Delta} + \frac{(\sigma_1^n)^2}{2 \Delta^2}\right)v_2^n-\frac{1}{h}v_{1}^{n}+\frac{1}{h}v_{1}^{n+1}\label{eq:A-collected-with-uniform-grid-left-t}\\
\intertext{Use $\partial_x v(t_n, \underline{x}) \approx \frac{v_1^n-v_0^n}{\Delta} = 0, \implies v_1^n = v_0^n$, so}
 &\approx \left(\underbrace{\left(-\frac{\mu_1^{n,-}}{\Delta} +\frac{(\sigma_1^n)^2}{2 \Delta^2}\right)}_{\equiv X_1} + \underbrace{\left(\frac{\mu_1^{n,-}-\mu_1^{n,+}}{\Delta}-\frac{(\sigma_1^n)^2}{\Delta^2}\right)}_{\equiv Y_1}\right)v_1^n + \underbrace{\left(\frac{\mu_1^{n,+}}{\Delta} + \frac{(\sigma_1^n)^2}{2 \Delta^2}\right)}_{\equiv Z_1}v_2^n-\frac{1}{h}v_{1}^{n}+\frac{1}{h}v_1^{n+1}
\intertext{And with a non-uniform grid, we recall the assumption that $\Delta_{1,-} = \Delta_{1,+}$.}
\mathcal{A} v(t_n,x_1) &\approx \left(-\frac{\mu_1^{n,-}}{\Delta_{1,+}} +\frac{(\sigma_1^n)^2}{2 \Delta_{1,+}^2}\right)v_0^n + \left(\frac{\mu_1^{n,-} - \mu_1^{n,+}}{\Delta_{1,+}}-\frac{(\sigma_1^n)^2}{\Delta_{1,+}^2}\right) v_1^n + \left(\frac{\mu_1^{n,+}}{\Delta_{1,+}} + \frac{(\sigma_1^n)^2}{2 \Delta_{1,+}^2}\right)v_2^n-\frac{1}{h_{n,+}}v_{1}^n+\frac{1}{h_{n,+}}v_1^{n+1}\label{eq:A-collected-with-nonuniform-grid-left-t}\\
 &\approx \left(\underbrace{\left(\frac{-\mu_1^{n,-}}{\Delta_{1,-}} +\frac{(\sigma_1^n)^2}{\Delta_{1,-}(\Delta_{1,+}+\Delta_{1,-})}\right)}_{\equiv X_1} + \underbrace{\left(\frac{\mu_1^{n,-} - \mu_1^{n,+}}{\Delta_{1,+}}-\frac{(\sigma_1^n)^2}{\Delta_{1,+}\Delta_{1,-}}\right)}_{\equiv Y_1}\right)v_1^n + \\ \nonumber
& \underbrace{\left(\frac{\mu_1^{n,+}}{\Delta_{1,+}} + \frac{(\sigma_1^n)^2}{2 \Delta_{1,+}\Delta_{1,-}}\right)}_{\equiv Z_1}v_2^n-\frac{1}{h_{n,+}}v_1^n+\frac{1}{h_{n,+}}v_1^{n+1}\label{eq:reflecting-left-nonuniform-explicit-t}
\end{align}

\paragraph{Boundary Value at $\bar{x}$:}
As $i=I$, 
\begin{align}
\mathcal{A} v(t_n,x_I)&\approx \left(-\frac{\mu_I^{n,-}}{\Delta} +\frac{(\sigma_I^n)^2}{2 \Delta^2}\right)v_{I-1}^n + \left(\frac{\mu_I^{n,-}-\mu_I^{n,+}}{\Delta}-\frac{(\sigma_I^n)^2}{\Delta^2}\right)v_I^n + \left(\frac{\mu_I^{n,+}}{\Delta} + \frac{(\sigma_I^n)^2}{2 \Delta^2}\right)v_{I+1}^n-\frac{1}{h}v_{I}^{n}+\frac{1}{h}v_{I}^{n+1}\\
\intertext{For a reflecting barrier, the boundary value $\partial_x v(t_n, \bar{x}) \approx \frac{v_{I+1}^n-v_I^n}{\Delta} = 0, \implies v_{I+1}^n = v_I^n$,}
&=\underbrace{\left(-\frac{\mu_I^{n,-}}{\Delta} +\frac{(\sigma_I^n)^2}{2 \Delta^2}\right)}_{\equiv X_I} v_{I-1}^n + \left(\underbrace{\left(\frac{\mu_I^{n,-}-\mu_I^{n,+}}{\Delta}-\frac{(\sigma_I^n)^2}{\Delta^2}\right)}_{\equiv Y_I} + \underbrace{\left(\frac{\mu_I^{n,+}}{\Delta} + \frac{(\sigma_I^n)^2}{2 \Delta^2}\right)}_{\equiv Z_I} \right) v_I^n -\frac{1}{h}v_{I}^n+\frac{1}{h}v_{I}^{n+1}
\intertext{And with a non-uniform grid, we recall the assumption that $\Delta_{I,+} = \Delta_{I,-}$.}
\mathcal{A}  v(t_n,x_I)&\approx \left(-\frac{\mu_I^{n,-}}{\Delta_{I,-}} +\frac{(\sigma_I^n)^2}{2 \Delta_{I,-}^2}\right)v_{I-1}^n + \left(\frac{\mu_I^{n,-} - \mu_I^{n,+}}{\Delta_{I,-}}-\frac{(\sigma_I^n)^2}{\Delta_{I,-}^2}\right)v_I^n + \left(\frac{\mu_I^{n,+}}{\Delta_{I,-}} + \frac{(\sigma_I^n)^2}{2 \Delta_{I,-}^2}\right)v_{I+1}^n-\frac{1}{h_{n,+}}v_I^n+\frac{1}{h_{n,+}}v_I^{n+1}\\
&\approx \underbrace{\left(-\frac{\mu_I^{n,-}}{\Delta_{I,-}} +\frac{(\sigma_I^n)^2}{\Delta_{I,-}(\Delta_{I,+}+\Delta_{I,-})}\right)}_{\equiv X_I} v_{I-1}^n\nonumber\\
&+ \left(\underbrace{\left(\frac{\mu_I^{n,-}-\mu_I^{n,+}}{\Delta_{I,+}}-\frac{(\sigma_I^n)^2}{\Delta_{I,-}\Delta_{I,+}}\right)}_{\equiv Y_I} + \underbrace{\left(\frac{\mu_I^{n,+}}{\Delta_{I,+}} + \frac{(\sigma_I^n)^2}{\Delta_{I,+}(\Delta_{I,+}+\Delta_{I,-})}\right)}_{\equiv Z_I} \right) v_I^n-\frac{1}{h_{n,+}}v_I^n+\frac{1}{h_{n,+}}v_I^{n+1} \label{eq:reflecting-right-nonuniform-explicit-t}
\end{align}

\subsection{Explicit Time Steps as a Matrix}

Define the matrix of zeros as $\mathbf{0} \equiv 0_{I\times I}$. Make the guess that the discretized operator is,
\begin{align}
A &\equiv \begin{bmatrix}
	A^1 - D_{+}^1	& D_{+}^1		& \mathbf{0}			 & \ldots 			& \mathbf{0}\\
	\mathbf{0} 		&A^2 - D_{+}^2 & D_{+}^2 				 & \ldots			 & \mathbf{0}\\	
	\vdots 			& \vdots 		& \vdots				& \vdots			&\vdots\\
	\mathbf{0} 		&\ldots 	 	&\mathbf{0}				& A^{N-1} - D_{+}^{N-1} & D_{+}^{N-1}\\
	\mathbf{0} 		&\ldots 		&						& \mathbf{0}	&A^N
	\end{bmatrix}\in\R^{(N I)\times (N I)}
\end{align}
Why is there no subtraction of the $D_{+}$ at time $T$?  It comes from the boundary condition that $\D[t]v(T,x) = 0$.  If things are setup properly, the bottom right of the setup should be solvable on its own, and exactly replicate the stationary setup.  Also notice that $A$ nests the stationary setup if there is only a single time-period.

\paragraph{Algebra on matrix $A$:}

Discretizing $v$ the same way as \cref{discretize}:

\begin{align}
v^n &\equiv \begin{bmatrix} v^n_1 &  v^n_2 & \ldots & v^n_I\end{bmatrix}^{\intercal}\in \R^I 
\intertext{Then stack these up over time,}
v &\equiv \begin{bmatrix} v^{1\intercal} &  v^{2\intercal} & \ldots&  v^{N\intercal}\end{bmatrix}^{\intercal}\in \R^{N I}
\end{align}

\paragraph{Interior with $1\leq n<N$:}\footnote{Notice here the case $n=1$ is the same as $1\leq n<N$ because we use explicit time, which means time only moves forwards, where $h_{1,+}=t_{2}-t_{1}$} Define new matrices $M^n$ for simplification:

\begin{align}
M^n &\equiv \begin{bmatrix} \mathbf{0} &  \ldots & A^n-D_{+}^n & D_{+}^n &\ldots &\mathbf{0}\end{bmatrix}\in \R^{I\times N I} 
\end{align}

 From \cref{eq:A-t} and \cref{eq:D-steady-state-T} we have:

\begin{align}
M^n v &\equiv \left(A^n-D_{+}^n\right)v^n+D_{+}^n v^{n+1}\\
=& \begin{bmatrix}
\left(Y^n_1 + X^n_1\right)v_1^n-\frac{1}{h_{n,+}}v_1^n & Z^n_1v_2^{n} & 0  & \cdots & \cdots & 0 \\
X^n_2v_1^n & \left(Y^n_2-\frac{1}{h_{n,+}}\right)v_2^n & Z^n_2v_3^{n} & \ddots& & \vdots \\
0 & \ddots & \ddots & \ddots &  & \vdots \\
\vdots &  & \ddots & \ddots & \ddots  & \vdots \\
\vdots & & \ddots & X^n_{I-1}v_{I-2}^n & \left(Y^n_{I-1}-\frac{1}{h_{n,+}}\right)v_{I-1}^n  & Z^n_{I-1}v_I^n \\
0 & \cdots & \cdots & 0 & X^n_Iv_{I-1}^{n} & \left(Y^n_I+Z^n_I-\frac{1}{h_{n,+}}\right)v_{I}^{n}\\
\end{bmatrix}\\
&+\begin{bmatrix} \frac{1}{h_{n,+}}v_1^{n+1} &\mathbf{0} &\ldots &\ldots &\mathbf{0}\\
\mathbf{0} & \frac{1}{h_{n,+}}v_2^{n+1} & \ldots &\ldots &\mathbf{0}\\
\vdots & \vdots & \vdots &\vdots & \vdots\\
\mathbf{0} & \ldots & \ldots & \frac{1}{h_{n,+}} v_{I-1}^{n+1} &\mathbf{0}\\
\mathbf{0} & \ldots &\ldots & \mathbf{0} & \frac{1}{h_{n,+}} v_{I}^{n+1} \end{bmatrix}
\end{align}

This is consistent with the results implied by (82)-(92).

\paragraph{Boundary with $n=N$:} In this case as we use the boundary condition $\D[t]v(T,x) = 0$, which implies $\frac{v_i^{N+1}-v_i^N}{h_{N,+}}=0$. From \cref{eq:A-collected-interior-with-uniform-grid-t} and \cref{eq:A-collected-interior-with-nonuniform-grid-t} we have:

\begin{align}
\mathcal{A} v(t_N,x_i) & = \underbrace{\left(-\frac{\mu_i^{N,-}}{\Delta} +\frac{(\sigma_i^N)^2}{2 \Delta^2}\right)}_{\equiv X_i}v_{i-1}^N + \underbrace{\left(\frac{\mu_i^{N,-}-\mu_i^{N,+}}{\Delta}-\frac{(\sigma_i^N)^2}{\Delta^2}\right)}_{\equiv Y_i}v_i^N + \underbrace{\left(\frac{\mu_i^{N,+}}{\Delta} + \frac{(\sigma_i^N)^2}{2 \Delta^2}\right)}_{\equiv Z_i}v_{i+1}^{N}
\intertext{And with a non-uniform grid}
\mathcal{A} v(t_N,x_i) & = \underbrace{\left(-\frac{\mu_i^{N,-}}{\Delta_{i,-}} +\frac{(\sigma_i^N)^2}{\Delta_{i,-}(\Delta_{i,+}+\Delta_{i,-})}\right)}_{\equiv X_i}v_{i-1}^N + \underbrace{\left(\frac{\mu_i^{N,-}}{\Delta_{i,-}}-\frac{\mu_i^{N,+}}{\Delta_{i,+}}-\frac{(\sigma_i^N)^2}{\Delta_{i,-}\Delta_{i,+}}\right)}_{\equiv Y_i}v_i^N + \\
&\underbrace{\left(\frac{\mu_i^{N,+}}{\Delta_{i,+}}+ \frac{(\sigma_i^N)^2}{\Delta_{i,+}(\Delta_{i,+}+\Delta_{i,-})}\right)}_{\equiv Z_i}v_{i+1}^N
\end{align}

\noindent This explains why in matrix $A$ we don't need the subtraction of $D_{+}$ at time $T$.\\
\noindent Using the matrix $A$ given above, we can solve the following,

\begin{align}
	\rho v(t,x) &= u(t,x) + \mathcal{A}v(t,x)\\
	\intertext{With the discretized $u\in\R^{N I}$, use this $A$ matrix.}
	\rho v &= u + A v
\intertext{subject to the appropriate boundary conditions.  Discretizing the grid, let $u \equiv \set{u(x_i)}_{i=1}^I\in\R^I$, then the discretized version is,}
(\rho  \mathbf{I} - A) v &= u\label{eq:discretized-linear-HJBE-dynamic}
\end{align}
%
%%Turned off until we can figure out if this is correct
%\subsection{Implicit in Time}
%\textbf{TODO:}\footnote{If I had to guess, implicit would end up something like:
%\begin{align}
%A &\equiv \begin{bmatrix}
%A^1 ???	& 0		& \mathbf{0}			 & \ldots 			& \mathbf{0}\\
%D_{-}^2 		&A^1 - D_{-}^2 & \mathbf{0}				 & \ldots			 & \mathbf{0}\\	
%\vdots 			& \vdots 		& \vdots				& \vdots			&\vdots\\
%\mathbf{0} 		&\ldots 	 	& D_{-}^{N-1}			& A^{N-1} - D_{-}^{N-1} &\mathbf{0}	\\
%\mathbf{0} 		&\ldots 		&						& \mathbf{0}	&A^N
%\end{bmatrix}\in\R^{(N I)\times (N I)}
%\end{align}
%The tough part to reason out is the $n=1$ setup at $t=0$, because you can't go backwards...  It is possible that this is simply not a good approach to solving the problem, and that the other approach is necessary?}
%The issue with implicit is what to do with the $t=0$ case?  Need to be very careful with the steady-state there.  It also may be trick to do the $t=T$ case as we may need ghost nodes to apply a boundary value?(e.g. maybe we have to add in the $N=1$ as the old stationary corner we have in the explicit method?\\
%\subsubsection{Algebra for Implicit Time Stepping}
%If we use implicit time stepping, then the discretization for uniform and non-uniform grids is a modification of  \cref{eq:A-generator-univariate-diffusion-with-nonuniform-grid-interior,eq:A-generator-univariate-diffusion-interior} in the interior $0 < n < N$ and $1 < i < I$:
%
%\begin{align}
%\mathcal{A} v(t_n, x_i) &\approx \underbrace{\left(v^n_i-v^n_{i-1}\right)\frac{\mu_i^{n,-}}{\Delta}+ \left(v^n_{i+1}-v^n_i\right)\frac{\mu_i^{n,+}}{\Delta}}_{\text{Upwind Scheme}}  + \frac{(\sigma^n_i)^2}{2} \frac{v^n_{i+1} - 2 v^n_i + v^n_{i-1}}{\Delta^2} + \underbrace{\frac{v^{n}_i - v^{n-1}_i}{h}}_{\text{Implicit}} \label{eq:A-generator-univariate-diffusion-interior-implicit-t}\\
%\intertext{And with a non-uniform grid}
%\mathcal{A} v(t_n, x_i) &\approx \underbrace{\left(v^n_i-v^n_{i-1}\right)\frac{\mu_i^{n, -}}{\Delta_{i,-}}+ \left(v^n_{i+1}-v^n_i\right)\frac{\mu_i^{n, +}}{\Delta_{i,+}}}_{\text{Upwind Scheme}}  + (\sigma^n_i)^2 \frac{\Delta_{i,-}v^n_{i+1} - (\Delta_{i,+}+\Delta_{i,-}) v^n_i + \Delta_{i,+}v^n_{i-1}}{\Delta_{i,-}\Delta_{i,+}(\Delta_{i,+}+\Delta_{i,-})}+ \underbrace{\frac{v^{n}_i - v^{n-1}_i}{h_{n,-}}}_{\text{Implicit}}\label{eq:A-generator-univariate-diffusion-with-nonuniform-grid-interior-implicit-t}
%\end{align}
%\paragraph{Interior of $A$ with implicit time stepping:}
%In the interior ($1 < i < I$), with a uniform grid, the discretization of \cref{eq:A-generator-univariate-diffusion-interior-implicit-t} and \cref{eq:A-generator-univariate-diffusion-with-nonuniform-grid-interior-implicit-t}
%\begin{align}
%\mathcal{A} v(t_n,x_i) & = \underbrace{\left(-\frac{\mu_i^{n,-}}{\Delta} +\frac{(\sigma_i^n)^2}{2 \Delta^2}\right)}_{\equiv X_i}v_{i-1}^n + \underbrace{\left(\frac{\mu_i^{n,-}-\mu_i^{n,+}}{\Delta}-\frac{(\sigma_i^n)^2}{\Delta^2}\right)}_{\equiv Y_i}v_i^n + \underbrace{\left(\frac{\mu_i^{n,+}}{\Delta} + \frac{(\sigma_i^n)^2}{2 \Delta^2}\right)}_{\equiv Z_i}v_{i+1}^{n}-\frac{1}{h}v_{i}^{n-1}+\frac{1}{h}v_{i}^{n}\label{eq:A-collected-interior-with-uniform-grid-implicit-t}
%\intertext{And with a non-uniform grid}
%\mathcal{A} v(t_n,x_i) & = \underbrace{\left(-\frac{\mu_i^{n,-}}{\Delta_{i,-}} +\frac{(\sigma_i^n)^2}{\Delta_{i,-}(\Delta_{i,+}+\Delta_{i,-})}\right)}_{\equiv X_i}v_{i-1}^n + \underbrace{\left(\frac{\mu_i^{n,-}}{\Delta_{i,-}}-\frac{\mu_i^{n,+}}{\Delta_{i,+}}-\frac{(\sigma_i^n)^2}{\Delta_{i,-}\Delta_{i,+}}\right)}_{\equiv Y_i}v_i^n + \nonumber\\
%&\underbrace{\left(\frac{\mu_i^{n,+}}{\Delta_{i,+}}+ \frac{(\sigma_i^n)^2}{\Delta_{i,+}(\Delta_{i,+}+\Delta_{i,-})}\right)}_{\equiv Z_i}v_{i+1}^n-\frac{1}{h_{n,-}}v_{i}^{n-1}+\frac{1}{h_{n,-}}v_{i}^{n}\label{eq:A-collected-interior-with-nonuniform-grid-implicit-t}
%\end{align}
%
%\paragraph{Boundary Value at $\underline{x}$:}
%As $i =1$, given the ``ghost node'' $x_0$, the discretized operator from \cref{eq:A-collected-interior-with-uniform-grid-implicit-t} is
%\begin{align}
%\mathcal{A} v(t_n,x_1) &\approx \left(-\frac{\mu_1^{n,-}}{\Delta} +\frac{(\sigma_1^n)^2}{2 \Delta^2}\right)v_0^n + \left(\frac{\mu_1^{n,-}-\mu_1^{n,+}}{\Delta}-\frac{(\sigma_1^n)^2}{\Delta^2}\right)v_1^n + \left(\frac{\mu_1^{n,+}}{\Delta} + \frac{(\sigma_1^n)^2}{2 \Delta^2}\right)v_2^n-\frac{1}{h}v_{1}^{n-1}+\frac{1}{h}v_{1}^{n}\label{eq:A-collected-with-uniform-grid-left-implicit-t}\\
%\intertext{Use $\partial_x v(t_n, \underline{x}) \approx \frac{v_1^n-v_0^n}{\Delta} = 0, \implies v_1^n = v_0^n$, so}
% &\approx \left(\underbrace{\left(-\frac{\mu_1^{n,-}}{\Delta} +\frac{(\sigma_1^n)^2}{2 \Delta^2}\right)}_{\equiv X_1} + \underbrace{\left(\frac{\mu_1^{n,-}-\mu_1^{n,+}}{\Delta}-\frac{(\sigma_1^n)^2}{\Delta^2}\right)}_{\equiv Y_1}\right)v_1^n + \underbrace{\left(\frac{\mu_1^{n,+}}{\Delta} + \frac{(\sigma_1^n)^2}{2 \Delta^2}\right)}_{\equiv Z_1}v_2^n-\frac{1}{h}v_{1}^{n-1}+\frac{1}{h}v_1^{n}
%\intertext{And with a non-uniform grid, we recall the assumption that $\Delta_{1,-} = \Delta_{1,+}$.}
%\mathcal{A} v(t_n,x_1) &\approx \left(-\frac{\mu_1^{n,-}}{\Delta_{1,+}} +\frac{(\sigma_1^n)^2}{2 \Delta_{1,+}^2}\right)v_0^n + \left(\frac{\mu_1^{n,-} - \mu_1^{n,+}}{\Delta_{1,+}}-\frac{(\sigma_1^n)^2}{\Delta_{1,+}^2}\right) v_1^n + \left(\frac{\mu_1^{n,+}}{\Delta_{1,+}} + \frac{(\sigma_1^n)^2}{2 \Delta_{1,+}^2}\right)v_2^n-\frac{1}{h_{n,-}}v_{1}^{n-1}+\frac{1}{h_{n,-}}v_1^{n}\label{eq:A-collected-with-nonuniform-grid-left-implicit-t}\\
% &\approx \left(\underbrace{\left(\frac{-\mu_1^{n,-}}{\Delta_{1,-}} +\frac{(\sigma_1^n)^2}{\Delta_{1,-}(\Delta_{1,+}+\Delta_{1,-})}\right)}_{\equiv X_1} + \underbrace{\left(\frac{\mu_1^{n,-} - \mu_1^{n,+}}{\Delta_{1,+}}-\frac{(\sigma_1^n)^2}{\Delta_{1,+}\Delta_{1,-}}\right)}_{\equiv Y_1}\right)v_1^n + \\ \nonumber
%& \underbrace{\left(\frac{\mu_1^{n,+}}{\Delta_{1,+}} + \frac{(\sigma_1^n)^2}{2 \Delta_{1,+}\Delta_{1,-}}\right)}_{\equiv Z_1}v_2^n-\frac{1}{h_{n,-}}v_1^{n-1}+\frac{1}{h_{n,-}}v_1^{n}\label{eq:reflecting-left-nonuniform-implicit-t}
%\end{align}
%
%\paragraph{Boundary Value at $\bar{x}$:}
%As $i=I$, 
%\begin{align}
%\mathcal{A} v(t_n,x_I)&\approx \left(-\frac{\mu_I^{n,-}}{\Delta} +\frac{(\sigma_I^n)^2}{2 \Delta^2}\right)v_{I-1}^n + \left(\frac{\mu_I^{n,-}-\mu_I^{n,+}}{\Delta}-\frac{(\sigma_I^n)^2}{\Delta^2}\right)v_I^n + \left(\frac{\mu_I^{n,+}}{\Delta} + \frac{(\sigma_I^n)^2}{2 \Delta^2}\right)v_{I+1}^n-\frac{1}{h}v_{I}^{n-1}+\frac{1}{h}v_{I}^{n}\\
%\intertext{For a reflecting barrier, the boundary value $\partial_x v(t_n, \bar{x}) \approx \frac{v_{I+1}^n-v_I^n}{\Delta} = 0, \implies v_{I+1}^n = v_I^n$,}
%&=\underbrace{\left(-\frac{\mu_I^{n,-}}{\Delta} +\frac{(\sigma_I^n)^2}{2 \Delta^2}\right)}_{\equiv X_I} v_{I-1}^n + \left(\underbrace{\left(\frac{\mu_I^{n,-}-\mu_I^{n,+}}{\Delta}-\frac{(\sigma_I^n)^2}{\Delta^2}\right)}_{\equiv Y_I} + \underbrace{\left(\frac{\mu_I^{n,+}}{\Delta} + \frac{(\sigma_I^n)^2}{2 \Delta^2}\right)}_{\equiv Z_I} \right) v_I^n -\frac{1}{h}v_{I}^{n-1}+\frac{1}{h}v_{I}^{n}
%\intertext{And with a non-uniform grid, we recall the assumption that $\Delta_{I,+} = \Delta_{I,-}$.}
%\mathcal{A}  v(t_n,x_I)&\approx \left(-\frac{\mu_I^{n,-}}{\Delta_{I,-}} +\frac{(\sigma_I^n)^2}{2 \Delta_{I,-}^2}\right)v_{I-1}^n + \left(\frac{\mu_I^{n,-} - \mu_I^{n,+}}{\Delta_{I,-}}-\frac{(\sigma_I^n)^2}{\Delta_{I,-}^2}\right)v_I^n + \left(\frac{\mu_I^{n,+}}{\Delta_{I,-}} + \frac{(\sigma_I^n)^2}{2 \Delta_{I,-}^2}\right)v_{I+1}^n-\frac{1}{h_{n,-}}v_I^{n-1}+\frac{1}{h_{n,-}}v_I^{n}\\
%&\approx \underbrace{\left(-\frac{\mu_I^{n,-}}{\Delta_{I,-}} +\frac{(\sigma_I^n)^2}{\Delta_{I,-}(\Delta_{I,+}+\Delta_{I,-})}\right)}_{\equiv X_I} v_{I-1}^n\nonumber\\
%&+ \left(\underbrace{\left(\frac{\mu_I^{n,-}-\mu_I^{n,+}}{\Delta_{I,+}}-\frac{(\sigma_I^n)^2}{\Delta_{I,-}\Delta_{I,+}}\right)}_{\equiv Y_I} + \underbrace{\left(\frac{\mu_I^{n,+}}{\Delta_{I,+}} + \frac{(\sigma_I^n)^2}{\Delta_{I,+}(\Delta_{I,+}+\Delta_{I,-})}\right)}_{\equiv Z_I} \right) v_I^n-\frac{1}{h_{n,-}}v_I^{n-1}+\frac{1}{h_{n,-}}v_I^{n} \label{eq:reflecting-right-nonuniform-implicit-t}
%\end{align}
%\subsubsection{Implicit Time Steps as a Matrix}
%On the left hand side, since we do not have the ghost time point $t_0$, we can assume $h_{1, -} = h_{1, +}$, which is similar with the logic we used to set $\Delta_{1, -} = \Delta_{1, +}$. Similarly, we are not able to have value points $\{v_i^0\}_{i=1}^I$ as well. We may assmue $v_i^1 - v_i^0 = v_i^2 - v_i^1$. The right end can be formulated by the same way as what we do for interiors. \\
%Define the matrix of zeros as $\mathbf{0} \equiv 0_{I\times I}$. Make the guess that the discretized operator is,
%\begin{align}
%A &\equiv \begin{bmatrix}
%	A^1 - D_{-}^1	&D_{-}^1			& \mathbf{0}			 & \ldots 			& \mathbf{0}\\
%	-D_{-}^2		 &A^2 + D_{-}^2 &\mathbf{0}  				 & \ldots			 & \mathbf{0}\\	
%	\vdots 			& \vdots 		& \vdots				& \vdots			&\vdots\\
%	\mathbf{0} 		&\ldots 	 	& -D_{-}^{N-1}				& A^{N-1} + D_{-}^{N-1} &\mathbf{0}\\
%	\mathbf{0} 		&\ldots 		&\ldots						& -D_{-}^{N}	&A^N+D_{-}^{N}
%	\end{bmatrix}\in\R^{(N I)\times (N I)}  \label{A-with-implicit-t}
%\end{align}
%When we are applying the boundary condition $\partial_t v(T, x) = 0$, A becomes
%\begin{align}
%A &\equiv \begin{bmatrix}
%	A^1 - D_{-}^1	&D_{-}^1			& \mathbf{0}			 & \ldots 			& \mathbf{0}\\
%	-D_{-}^2		 &A^2 + D_{-}^2 &\mathbf{0}  				 & \ldots			 & \mathbf{0}\\	
%	\vdots 			& \vdots 		& \vdots				& \vdots			&\vdots\\
%	\mathbf{0} 		&\ldots 	 	& -D_{-}^{N-1}				& A^{N-1} + D_{-}^{N-1} &\mathbf{0}\\
%	\mathbf{0} 		&\ldots 		&\ldots						& \mathbf{0}	&A^N
%	\end{bmatrix}\in\R^{(N I)\times (N I)}
%\end{align}
%
%
%Why is there no subtraction of the $D_{+}$ at time $T$?  It comes from the boundary condition that $\D[t]v(T,x) = 0$.  If things are setup properly, the bottom right of the setup should be solvable on its own, and exactly replicate the stationary setup.  Also notice that $A$ nests the stationary setup if there is only a single time-period. 
%
%\paragraph{Algebra on matrix $A$:}
%Instead of losing generality, following equations are formulated in terms of \eqref{A-with-implicit-t}. Discretizing $v$ the same way as \cref{discretize}:
%
%\begin{align}
%v^n &\equiv \begin{bmatrix} v^n_1 &  v^n_2 & \ldots & v^n_I\end{bmatrix}^{\intercal}\in \R^I 
%\intertext{Then stack these up over time,}
%v &\equiv \begin{bmatrix} v^{1\intercal} &  v^{2\intercal} & \ldots&  v^{N\intercal}\end{bmatrix}^{\intercal}\in \R^{N I}
%\end{align}
%
%\paragraph{Interior with $1< n < N$ and Right end with $n = N$:}\footnote{Notice here the case $n=N$ is the same as $1\leq n<N$ because we are using implicit time, which means time only moves backwards, where $h_{N,-}=t_{N}-t_{N-1}$} Define new matrices $M^n$ for simplification:
%
%\begin{align}
%M^n &\equiv \begin{bmatrix} \mathbf{0} & \ldots &-D_{-}^n  & A^n+D_{-}^n &\ldots &\mathbf{0}\end{bmatrix}\in \R^{I\times N I} 
%\end{align}
%
% From \cref{eq:A-t} and \cref{eq:D-steady-state-T} we have:
%
%\begin{align}
%M^n v &\equiv \left(A^n+D_{-}^{n}\right)v^n-D_{-}^n v^{n-1}\\
%=& \begin{bmatrix}
%\left(Y^n_1 + X^n_1+\frac{1}{h_{n,-}}\right)v_1^{n} & Z^n_1v_2^{n} & 0  & \cdots & \cdots & 0 \\
%X^n_2v_1^n & \left(Y^n_2+\frac{1}{h_{n,-}}\right)v_2^n & Z^n_2v_3^{n} & \ddots& & \vdots \\
%0 & \ddots & \ddots & \ddots &  & \vdots \\
%\vdots &  & \ddots & \ddots & \ddots  & \vdots \\
%\vdots & & \ddots & X^n_{I-1}v_{I-2}^n & \left(Y^n_{I-1}+\frac{1}{h_{n,-}}\right)v_{I-1}^n  & Z^n_{I-1}v_I^n \\
%0 & \cdots & \cdots & 0 & X^n_Iv_{I-1}^{n} & \left(Y^n_I+Z^n_I+\frac{1}{h_{n,-}}\right)v_{I}^{n}\\
%\end{bmatrix} \nonumber \\
%&-\begin{bmatrix} \frac{1}{h_{n,-}}v_1^{n-1} &\mathbf{0} &\ldots &\ldots &\mathbf{0}\\
%\mathbf{0} & \frac{1}{h_{n,-}}v_2^{n-1} & \ldots &\ldots &\mathbf{0}\\
%\vdots & \vdots & \vdots &\vdots & \vdots\\
%\mathbf{0} & \ldots & \ldots & \frac{1}{h_{n,-}} v_{I-1}^{n-1} &\mathbf{0}\\
%\mathbf{0} & \ldots &\ldots & \mathbf{0} & \frac{1}{h_{n,-}} v_{I}^{n-1} \end{bmatrix}
%\end{align}
%
%\paragraph{Boundary with $n=1$:} Based on assumptions we proposed above, we have $\frac{v_i^1-v_i^0}{h_{1, -}} = \frac{v_i^2-v_i^1}{h_{1, +}}$. From \cref{eq:A-collected-interior-with-uniform-grid-implicit-t} and \cref{eq:A-collected-interior-with-nonuniform-grid-implicit-t} we have:
% 
% \begin{align}
%\mathcal{A} v(t_1,x_i) & = \underbrace{\left(-\frac{\mu_i^{1,-}}{\Delta} +\frac{(\sigma_i^1)^2}{2 \Delta^2}\right)}_{\equiv X_i}v_{i-1}^1 + \underbrace{\left(\frac{\mu_i^{1,-}-\mu_i^{1,+}}{\Delta}-\frac{(\sigma_i^1)^2}{\Delta^2}\right)}_{\equiv Y_i}v_i^1 + \underbrace{\left(\frac{\mu_i^{1,+}}{\Delta} + \frac{(\sigma_i^1)^2}{2 \Delta^2}\right)}_{\equiv Z_i}v_{i+1}^{1}+\frac{v_{i}^{2}-v_{i}^{1}}{h}
%\intertext{And with a non-uniform grid}
%\mathcal{A} v(t_1,x_i) & = \underbrace{\left(-\frac{\mu_i^{1,-}}{\Delta_{i,-}} +\frac{(\sigma_i^1)^2}{\Delta_{i,-}(\Delta_{i,+}+\Delta_{i,-})}\right)}_{\equiv X_i}v_{i-1}^1 + \underbrace{\left(\frac{\mu_i^{1,-}}{\Delta_{i,-}}-\frac{\mu_i^{1,+}}{\Delta_{i,+}}-\frac{(\sigma_i^1)^2}{\Delta_{i,-}\Delta_{i,+}}\right)}_{\equiv Y_i}v_i^1 + \nonumber\\
%&\underbrace{\left(\frac{\mu_i^{1,+}}{\Delta_{i,+}}+ \frac{(\sigma_i^1)^2}{\Delta_{i,+}(\Delta_{i,+}+\Delta_{i,-})}\right)}_{\equiv Z_i}v_{i+1}^1+\frac{v_{i}^{2}-v_{i}^{1}}{h_{1,+}}
%\end{align}




\bibliography{etk-references}

\end{document}