Numerical Schemes

Finite Difference Policy Iteration

Discretize the HJB on a grid. The theory of viscosity solutions [@crandall1983viscosity, @user2013users] guarantees convergence of monotone schemes [@barles1991convergence]. Practical implementations follow [@forsyth2007numerical, @kushner2001numerical].

At each time step \(t_n\):

  1. Policy evaluation: solve \(A(u^{(k)}) V^{(k)} = V^{n+1} - c(u^{(k)})\) (linear system).
  2. Policy improvement [@howard1960dynamic]: \(u_i^{(k+1)} \in \arg\max_u \lbrace [A(u)V^{(k)}]_i + c(u)_i \rbrace\).
  3. Iterate until \(\|u^{(k+1)} - u^{(k)}\| < \text{tol}\).

Upwind Discretization

First derivatives:

\[D_x V \approx \begin{cases} \dfrac{V(x + \Delta x) - V(x)}{\Delta x}, & \text{drift}(x, u) > 0 \\[8pt] \dfrac{V(x) - V(x - \Delta x)}{\Delta x}, & \text{drift}(x, u) < 0 \end{cases}\]

Second derivatives (central difference):

\[D_x^2 V \approx \frac{V(x + \Delta x) - 2V(x) + V(x - \Delta x)}{\Delta x^2}\]

The CFL condition [@courant1928partiellen] constrains the time step for explicit schemes. Rannacher smoothing [@rannacher1984finite] is used for Crank-Nicolson startup.

Linear solvers: SOR, Thomas (tridiagonal), LAPACK dgtsv. Schemes: Implicit, Explicit, Crank-Nicolson, Strang ADI.

Implementation: solver::numeric::finite_difference

BSDE Least-Squares Monte Carlo

Backward SDEs provide a probabilistic representation of the HJB solution [@pardoux1990adapted, @el1997backward]. The regression-based scheme follows [@gobet2005empirical].

Forward SDE (Euler-Maruyama)

The forward process is simulated under the optimal control \(u^*\):

\[dX_t = b(t, X_t, u^*_t)\,dt + \sigma(t, X_t)\,dW_t\]

Discretized as

\[X_{n+1} = X_n + b(t_n, X_n, u^*_n)\Delta t + \sigma(t_n, X_n)\Delta W_n.\]

Backward SDE

For the controlled problem with running reward \(f\) and terminal cost \(g\), the value \(Y_t = V(t, X_t)\) satisfies the nonlinear Feynman-Kac relation [@pardoux1990adapted, @el1997backward]:

\[-dY_t = f(t, X_t, u^*_t)\,dt - Z_t\,dW_t, \qquad Y_T = g(X_T),\]

with \(Z_t = \sigma(t, X_t)^\top D_x V(t, X_t)\).

Equivalently

\[V(t, x) = \mathbb{E}\!\left[g(X_T) + \int_t^T f(s, X_s, u^*_s)\,ds \;\middle\|\; X_t = x\right].\]

The backward driver is the running reward \(f\) alone, not the full HJB generator \(f + \mathcal{L}^u V\). The infinitesimal generator \(\mathcal{L}^u V = b\cdot\nabla V + \tfrac12\mathrm{tr}(\sigma\sigma^T D^2 V)\) is already accounted for by simulating the forward SDE under the control; adding it again in the backward step double-counts the drift and diffusion.

LSMC Algorithm

  1. Simulate \(M\) forward paths under the optimal control.
  2. Set \(Y_N^{(m)} = g(X_N^{(m)})\) (terminal condition).
  3. For \(n = N-1, \dots, 0\):
    • Compute the regression target \(\text{Target}^{(m)} = Y_{n+1}^{(m)} + f(t_n, X_n^{(m)}, u^*_n)\Delta t\).
    • Regress \(\text{Target}^{(m)}\) onto basis functions \(\psi_k(X_n^{(m)})\): \[c^n = \arg\min_c \sum_{m=1}^M \left(\text{Target}^{(m)} - \sum_{k=1}^K c_k \psi_k(X_n^{(m)})\right)^2\]
    • Set \(Y_n^{(m)} = \sum_k c_k^n \psi_k(X_n^{(m)})\).

Basis types: Power, Hermite, Chebyshev, Laguerre. Feature count for degree-2: \(K = \frac{(N+1)(N+2)}{2}\) (6 for N=2, 10 for N=3).

See [@han2018solving, @han2017deep] for extensions with neural network approximation (deep BSDE).

Implementation: solver::numeric::bsde