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

\begin{document}
\title{Discretizing Continuous-Time Operators with Finite Differences}
\author{<NAME>\\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/}. }  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 &= - \mu^{-}+ \frac{\sigma^{2}}{2 \Delta}\label{eq:X} \\
	Y &= - \mu^{+} + \mu^{-}- \frac{\sigma^{2}}{\Delta}\label{eq:Y} \\
	Z &= \mu^{+} + \frac{\sigma^{2}}{2 \Delta}\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 \frac{A}{\Delta} v + \frac{b}{\Delta}
\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}
\Delta \mathcal{A} v(x_i) &\approx \underbrace{\left(v_i-v_{i-1}\right)\mu_i^{-}+ \left(v_{i+1}-v_i\right)\mu_i^{+}}_{\text{Upwind Scheme}}  + \frac{\sigma_i^2}{2} \frac{v_{i+1} - 2 v_i + v_{i-1}}{\Delta}\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(-\mu_i^{-} +\frac{\sigma_i^2}{2\Delta}\right)}_{\equiv X_i}v_{i-1} + \underbrace{\left(\mu_i^{-} - \mu_i^{+}-\frac{\sigma_i^2}{\Delta}\right)}_{\equiv Y_i}v_i + \underbrace{\left(\mu_i^{+} + \frac{\sigma_i^2}{2\Delta}\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}
\Delta \mathcal{A} v(x_1) &\approx \left(-\mu_1^{-} +\frac{\sigma_1^2}{2\Delta}\right)v_0 + \left(\mu_1^{-} - \mu_1^{+}-\frac{\sigma_1^2}{\Delta}\right)v_1 + \left(\mu_1^{+} + \frac{\sigma_1^2}{2\Delta}\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(-\mu_1^{-} +\frac{\sigma_1^2}{2\Delta}\right)}^{\equiv X_1}\underline{v}}_{\equiv b_1} +  \underbrace{\left(\mu_1^{-} - \mu_1^{+}-\frac{\sigma_1^2}{\Delta}\right)}_{\equiv Y_1}v_1 + \underbrace{\left(\mu_1^{+} + \frac{\sigma_1^2}{2\Delta}\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}
\Delta \mathcal{A} v(x_1) &\approx \left(\underbrace{\left(-\mu_1^{-} +\frac{\sigma_1^2}{2\Delta}\right)}_{\equiv X_1} + \underbrace{\left(\mu_1^{-} - \mu_1^{+}-\frac{\sigma_1^2}{\Delta}\right)}_{\equiv Y_1}\right)v_1 + \underbrace{\left(\mu_1^{+} + \frac{\sigma_1^2}{2\Delta}\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}
\Delta \mathcal{A} v(x_I)&\approx \left(-\mu_i^{-} +\frac{\sigma_I^2}{2\Delta}\right)v_{I-1} + \left(\mu_i^{-} - \mu_i^{+}-\frac{\sigma_I^2}{\Delta}\right)v_I + \left(\mu_i^{+} + \frac{\sigma_I^2}{2\Delta}\right)v_{I+1}\\
\intertext{For the absorbing barrier, substitute for $v(x_{I+1}) = \bar{v}$,}
\Delta \mathcal{A} v(x_I)&\approx \underbrace{\left(-\mu_i^{-} +\frac{\sigma_I^2}{2\Delta}\right)}_{\equiv X_I} v_{I-1} + \underbrace{\left(\mu_i^{-} - \mu_i^{+}-\frac{\sigma_I^2}{\Delta}\right)}_{\equiv Y_I} v_I + \underbrace{\overbrace{\left(\mu_i^{+} + \frac{\sigma_I^2}{2\Delta}\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$,}
\Delta \mathcal{A} v(x_I)&\approx \underbrace{\left(-\mu_i^{-} +\frac{\sigma_I^2}{2\Delta}\right)}_{\equiv X_I} v_{I-1} + \left(\underbrace{\left(\mu_i^{-} - \mu_i^{+}-\frac{\sigma_I^2}{\Delta}\right)}_{\equiv Y_I} + \underbrace{\left(\mu_i^{+} + \frac{\sigma_I^2}{2\Delta}\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 &= - \mu+ \frac{\sigma^{2}}{2 \Delta}\label{eq:X-backwards} \\
	Y &= \mu - \frac{\sigma^{2}}{\Delta}\label{eq:Y-backwards} \\
	Z &= \frac{\sigma^{2}}{2 \Delta}\label{eq:Z-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_i + \frac{\sigma_i^2}{2} \left(v_{i+1} - 2 v_i + v_{i-1}\right),\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{\sigma^{2}}{2\Delta}v_1 + \frac{\sigma^{2}}{2\Delta} v_2\\
\Delta \mathcal{A} v(x_I)&\approx \left(-\mu_I +\frac{\sigma_I^2}{2\Delta}\right) v_{I-1} + \left(\mu_I - \frac{\sigma_I^2}{2\Delta}\right) v_I 
\end{align}

\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 &= \frac{A^T}{\Delta} 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 multiply by $\Delta$ (which does not change problem for a $0$ eigenvalue) and then 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 \\ \Delta\end{bmatrix}\in\R^{I + 1}\\
X &\equiv \begin{bmatrix}
	A^T\\
	\Delta \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}
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,}
	\Delta \rho v &= \Delta u + A v\\
	\intertext{Or,}
	(\Delta \rho \mathbf{I} - A) v &= \Delta 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}\Delta u \\0 \\ \vdots \\ 0 \\ \Delta 1\end{bmatrix}\in\R^{2 I + 1}\\
X &\equiv \begin{bmatrix}
	\rho \Delta  \mathbf{I} - A & \mathbf{0}_{I\times I}\\
	\mathbf{0}_{I\times I} & A^T\\
	\mathbf{0}_{1\times I} & \Delta \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}

\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.  In that case, the SDE is
\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.


\bibliography{etk-references}

\end{document}<commit_msg>Put in placeholder for non-uniform grid<commit_after>\documentclass[11pt]{etk-article}
\usepackage{pstool} 
\usepackage{etk-bib}
\pdfmetadata{}{}{}{}

\begin{document}
\title{Discretizing Continuous-Time Operators with Finite Differences}
\author{<NAME>\\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/}. }  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}\label{sec:upwind-uniform}
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 &= - \mu^{-}+ \frac{\sigma^{2}}{2 \Delta}\label{eq:X} \\
	Y &= - \mu^{+} + \mu^{-}- \frac{\sigma^{2}}{\Delta}\label{eq:Y} \\
	Z &= \mu^{+} + \frac{\sigma^{2}}{2 \Delta}\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 \frac{A}{\Delta} v + \frac{b}{\Delta}\label{eq:A-operator-approx}
\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}
\Delta \mathcal{A} v(x_i) &\approx \underbrace{\left(v_i-v_{i-1}\right)\mu_i^{-}+ \left(v_{i+1}-v_i\right)\mu_i^{+}}_{\text{Upwind Scheme}}  + \frac{\sigma_i^2}{2} \frac{v_{i+1} - 2 v_i + v_{i-1}}{\Delta}\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(-\mu_i^{-} +\frac{\sigma_i^2}{2\Delta}\right)}_{\equiv X_i}v_{i-1} + \underbrace{\left(\mu_i^{-} - \mu_i^{+}-\frac{\sigma_i^2}{\Delta}\right)}_{\equiv Y_i}v_i + \underbrace{\left(\mu_i^{+} + \frac{\sigma_i^2}{2\Delta}\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}
\Delta \mathcal{A} v(x_1) &\approx \left(-\mu_1^{-} +\frac{\sigma_1^2}{2\Delta}\right)v_0 + \left(\mu_1^{-} - \mu_1^{+}-\frac{\sigma_1^2}{\Delta}\right)v_1 + \left(\mu_1^{+} + \frac{\sigma_1^2}{2\Delta}\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(-\mu_1^{-} +\frac{\sigma_1^2}{2\Delta}\right)}^{\equiv X_1}\underline{v}}_{\equiv b_1} +  \underbrace{\left(\mu_1^{-} - \mu_1^{+}-\frac{\sigma_1^2}{\Delta}\right)}_{\equiv Y_1}v_1 + \underbrace{\left(\mu_1^{+} + \frac{\sigma_1^2}{2\Delta}\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}
\Delta \mathcal{A} v(x_1) &\approx \left(\underbrace{\left(-\mu_1^{-} +\frac{\sigma_1^2}{2\Delta}\right)}_{\equiv X_1} + \underbrace{\left(\mu_1^{-} - \mu_1^{+}-\frac{\sigma_1^2}{\Delta}\right)}_{\equiv Y_1}\right)v_1 + \underbrace{\left(\mu_1^{+} + \frac{\sigma_1^2}{2\Delta}\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}
\Delta \mathcal{A} v(x_I)&\approx \left(-\mu_i^{-} +\frac{\sigma_I^2}{2\Delta}\right)v_{I-1} + \left(\mu_i^{-} - \mu_i^{+}-\frac{\sigma_I^2}{\Delta}\right)v_I + \left(\mu_i^{+} + \frac{\sigma_I^2}{2\Delta}\right)v_{I+1}\\
\intertext{For the absorbing barrier, substitute for $v(x_{I+1}) = \bar{v}$,}
\Delta \mathcal{A} v(x_I)&\approx \underbrace{\left(-\mu_i^{-} +\frac{\sigma_I^2}{2\Delta}\right)}_{\equiv X_I} v_{I-1} + \underbrace{\left(\mu_i^{-} - \mu_i^{+}-\frac{\sigma_I^2}{\Delta}\right)}_{\equiv Y_I} v_I + \underbrace{\overbrace{\left(\mu_i^{+} + \frac{\sigma_I^2}{2\Delta}\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$,}
\Delta \mathcal{A} v(x_I)&\approx \underbrace{\left(-\mu_i^{-} +\frac{\sigma_I^2}{2\Delta}\right)}_{\equiv X_I} v_{I-1} + \left(\underbrace{\left(\mu_i^{-} - \mu_i^{+}-\frac{\sigma_I^2}{\Delta}\right)}_{\equiv Y_I} + \underbrace{\left(\mu_i^{+} + \frac{\sigma_I^2}{2\Delta}\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 &= - \mu+ \frac{\sigma^{2}}{2 \Delta}\label{eq:X-backwards} \\
	Y &= \mu - \frac{\sigma^{2}}{\Delta}\label{eq:Y-backwards} \\
	Z &= \frac{\sigma^{2}}{2 \Delta}\label{eq:Z-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_i + \frac{\sigma_i^2}{2} \left(v_{i+1} - 2 v_i + v_{i-1}\right),\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{\sigma^{2}}{2\Delta}v_1 + \frac{\sigma^{2}}{2\Delta} v_2\\
\Delta \mathcal{A} v(x_I)&\approx \left(-\mu_I +\frac{\sigma_I^2}{2\Delta}\right) v_{I-1} + \left(\mu_I - \frac{\sigma_I^2}{2\Delta}\right) v_I 
\end{align}

\subsection{Upwind Finite Difference Discretization for $A$ with a Non-Uniform Grid}\label{sec:upwind-non-uniform}
Redoing section \cref{sec:upwind-uniform} with a nonuniform grid, from \cref{eq:Delta-i-plus,eq:Delta-i-minus} define the vectors,
 \begin{align}
\Delta_{+}&\equiv \set{\Delta_{i,+}}_{i=1}^I = \set{x_{i+1} - x_i}_{i=1}^{I}\label{eq:Delta-plus}\\
\Delta_{-}&\equiv \set{\Delta_{i,-}}_{i=1}^I = \set{x_i - x_{i-1}}_{i=1}^I\label{eq:Delta-minus}	
\end{align}
Note here that for simplicity of notation the $\Delta_{+}$ node goes beyond $i=I$ and the $\Delta_{-}$ goes before $i=1$.  Hence, it is convenient to specify ``ghost-nodes'' $x_0$ and $x_{I+1}$.  While we could set them to arbitrary values, it is probably easiest to just let $x_0 = x_1 - \Delta_{1,+}$ and $x_{I+1} = x_I + \Delta_{I,-}$.

Only solving the versions with reflecting barriers (i.e. $b=0$), the modified versions of \cref{eq:X,eq:Y,eq:Z,eq:A} are, \textbf{TODO: THESE WERE JUST A FEW GUESSES...}
\begin{align}
X &= - \mu^{-}+ \frac{\sigma^{2}}{\Delta_{+} + \Delta_{-}}????\label{eq:X-nonuniform} \\
Y &= - \mu^{+} + \mu^{-}- \frac{\sigma^{2}}{(\Delta_{+} + \Delta_{-})/2}???\label{eq:-nonuniform} \\
Z &= \mu^{+} + \frac{\sigma^{2}}{\Delta_{+} + \Delta_{-}}???\label{eq:Z-nonuniform}\\
A &= \ldots \label{eq:A-nonuniform}
\end{align}
And the approximation is,\textbf{TODO:} Incorrect, need to take into account dividing by vector $\Delta$.
\begin{align}
\mathcal{A}v(\set{x_i}) &\approx \frac{A}{\Delta} v ???\label{eq:A-operator-approx-nonuniform}
\end{align}

\paragraph{Interior of $A$:} Repeating the analysis of \cref{eq:A-generator-univariate-diffusion-interior}

\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 &= \frac{A^T}{\Delta} 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 multiply by $\Delta$ (which does not change problem for a $0$ eigenvalue) and then 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 \\ \Delta\end{bmatrix}\in\R^{I + 1}\\
X &\equiv \begin{bmatrix}
	A^T\\
	\Delta \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}
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,}
	\Delta \rho v &= \Delta u + A v\\
	\intertext{Or,}
	(\Delta \rho \mathbf{I} - A) v &= \Delta 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}\Delta u \\0 \\ \vdots \\ 0 \\ \Delta 1\end{bmatrix}\in\R^{2 I + 1}\\
X &\equiv \begin{bmatrix}
	\rho \Delta  \mathbf{I} - A & \mathbf{0}_{I\times I}\\
	\mathbf{0}_{I\times I} & A^T\\
	\mathbf{0}_{1\times I} & \Delta \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}

\section{Time-Independent Univariate Diffusions with Non-Uniform Grid}\label{sec:time-independent-univariate-diffusion-nonuniform}
This section does a version of the discretization with a non-uniform grid
\begin{itemize}
	\item Discretization
\end{itemize}
\subsection{}

\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.  In that case, the SDE is
\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.


\bibliography{etk-references}

\end{document}
{'fraction_non_alphanumeric': 0.16987479874077813, 'fraction_numerical': 0.012063537836733713, 'mean_word_length': 4.391112838450576, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 9, 'https://': 10, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4258407', 'n_tokens_mistral': 14679, 'n_tokens_neox': 14121, 'n_words': 4595}