Nonlinear Black–Scholes: Implementing the Glover–Duck–Newton Full-Feedback PDE

1. Background: why the linear model is not enough

The classical Black–Scholes model prices a European option by assuming the underlying asset price follows geometric Brownian motion and that the hedging trades of any single agent have no effect on the market. Under these assumptions the option price V(S,t) satisfies the linear PDE

where r is the risk-free rate, q the dividend yield, and \sigma the constant volatility. This equation is linear in V: the coefficient \tfrac{1}{2}\sigma^2 S^2 on the second derivative term does not depend on the solution itself.

In reality, a sufficiently large trader — a fund delta-hedging a sizable options position — can move the market with their own trades. When the trader buys stock to delta-hedge a long call, that buying pressure pushes the stock price up, which changes the delta, which requires buying more stock, and so on. This feedback loop breaks the no-impact assumption and must be modelled explicitly.

Glover, Duck & Newton (2009) studied a class of such models in detail. The key quantity driving the nonlinearity is the option's gamma, V_{SS} = \partial^2 V / \partial S^2, which measures the convexity of the option value with respect to the stock price. Large gamma means the hedging strategy changes rapidly with S, amplifying the price-impact feedback loop.

2. The full-feedback nonlinear PDE

Following Schönbucher & Wilmott (2000) and the canonical form studied in Glover, Duck & Newton (2009, eq. 2.10), when the hedger's trading strategy is based on the actual delta of the modified option price (full feedback), the option price satisfies the fully nonlinear PDE

The only structural difference from classical Black–Scholes is the denominator (1 - \lambda\,V_{SS})^2 inside the diffusion term. In this model, \lambda \geq 0 is a constant dimensionless liquidity parameter (Schönbucher & Wilmott 2000; Glover et al. 2009, p. 5): it does not carry an explicit factor of S. When \lambda = 0 the denominator equals 1 and the equation reduces to classical Black–Scholes exactly.

The parameter \lambda controls the strength of price impact. As \lambda increases, the effective volatility

is amplified when V_{SS} > 0 (positive gamma, normal for long options) and compressed when V_{SS} < 0. The denominator vanishes — a singularity corresponding to an infinitely illiquid market — when V_{SS} = 1/\lambda. Glover et al. (2009, Section 5) show that this singularity can be induced even by smooth payoffs, which is one of the principal cautionary findings of their paper.

The nonlinearity lives entirely in the diffusion coefficient. Defining

the PDE takes the same structural form as linear Black–Scholes,

with q(S) = rS and r(S) = -r, except that \tilde{p} now depends on V_{SS}. That circular dependency is what makes the equation nonlinear and requires an iterative solver at each timestep.

3. The existing solver: structure of parab.m

The baseline code solves linear Black–Scholes using a Crank–Nicolson finite-difference scheme on a nonuniform spatial grid. The key features relevant to our extension are:

• Spatial discretisation. The functions cfd2mat.m and nuni.m build second-order centred-difference matrices A_2, A_1, A_0 and an extended matrix A_{2,\text{ext}} that maps the full-grid solution vector to second derivatives at interior nodes. Support for three boundary condition modes (Direxact, Dirdisc, Linear) is already in place.
• Crank–Nicolson timestepping. At each step the elliptic operator is evaluated as L = pA_2 + qA_1 + rA_0 with coefficient vectors pvct, qvct, rvct. The implicit-side matrix is A_0 - \theta \Delta t \, L and the explicit-side multiplier is A_0 + (1-\theta)\Delta t \, L.
• American option support. penalty_solver.m enforces the early-exercise constraint by adding a large diagonal penalty to the linear system wherever the solution falls below the payoff.
• Rannacher / adaptive timestepping. The first few steps use implicit (backward-Euler, \theta = 1) half-steps to damp oscillations caused by non-smooth initial data, then switch to Crank–Nicolson (\theta = 1/2). An adaptive step-size controller adjusts \Delta t based on a per-component change ratio.

All of this infrastructure works equally for the nonlinear model — we only need to change how pvct is computed at each timestep, and add a short iteration loop around the linear solve.

4. Solving the nonlinearity: Picard iteration

4.1 Why the equation cannot be solved directly

After Crank–Nicolson discretisation in time, each timestep requires solving a system of the form

The matrix A_0 - \theta \Delta t \, L(\mathbf{V}^{n+1}) on the left depends on \mathbf{V}^{n+1} through the coefficient vector \tilde{p}, which involves V_{SS}^{n+1} = A_{2,\text{ext}}\,\mathbf{V}^{n+1}. So the system is nonlinear in \mathbf{V}^{n+1} and cannot be solved with a single sparse factorisation.

4.2 Why we chose Picard, not Newton

Two standard iterative approaches exist for this type of problem:

We therefore use Picard iteration, also called the lagged-coefficient method in the PDE literature.

4.3 The combined Picard + penalty algorithm

Rather than nesting a penalty loop inside a Picard loop, both are collapsed into a single fixed-point iteration. At each iteration k, the nonlinear diffusion coefficient \tilde{p}^{(k)} and the penalty matrix P^{(k)} are both lagged from the previous iterate \mathbf{V}^{(k)}, and a single linear system is solved:

Convergence is declared when the solution change falls below tolerance and the penalty active set (the set of nodes where V < \psi) stops changing between iterates.

Let \mathbf{V}^n denote the solution at time t^n. At each timestep we iterate:

1. Initial guess. Set V_{SS}^{(0)} = A_{2,\text{ext}}\,\mathbf{V}^n and P^{(0)} = \mathrm{Large}\cdot\mathbf{1}[V^n < \psi] — gamma and penalty active set both from the previous timestep. This is an O(\Delta t)-accurate starting point.

2. Update diffusion coefficient. Compute the nonlinear coefficient at every grid node using the lagged gamma (constant \lambda, no S factor):

3. Rebuild the operator. Form L^{(k)} = \tilde{P}^{(k)} A_2 + Q A_1 + R A_0 and assemble AA_{nl}^{(k)} = A_0 - \theta \Delta t \, L^{(k)}.

4. Solve the combined linear system. A single sparse solve simultaneously advances the iterate and enforces the early-exercise constraint.

5. Update gamma and penalty. Set V_{SS}^{(k+1)} = A_{2,\text{ext}}\,\mathbf{V}^{(k+1)} and P^{(k+1)} = \mathrm{Large}\cdot\mathbf{1}[V^{(k+1)} < \psi].

6. Check convergence. Stop if \|\mathbf{V}^{(k+1)} - \mathbf{V}^{(k)}\|_\infty / \max(1, \|\mathbf{V}^{(k+1)}\|_\infty) < \varepsilon and the active set has not changed, otherwise return to step 2.

In practice, convergence in 2–3 iterations is the norm for \lambda \in [0.01, 0.10], confirming the analysis of Koleva and Vulkov (2013).

4.4 Why the previous-timestep gamma is a good initial guess

The smoothness of the solution in time is the key. Assuming a sufficiently small timestep,

so the initial Picard error is already O(\Delta t). The lagged-coefficient operator also contracts with a factor proportional to \lambda, which is small in practice. Koleva and Vulkov (2013, Section 3) prove convergence of this iteration for \lambda small enough that the denominator 1 - \lambda V_{SS} stays bounded away from zero — the same singularity condition identified by Glover et al. (2009) in their analytical study.

5. Code changes: what was added and why

Three new files were written. The original six files (parab.m, penalty_solver.m, cfd2mat.m, rhscfd2pb.m, nuni.m, fb_quad.m) are unmodified.

5.1 nlbs_diffcoeff.m — the nonlinear coefficient

This function computes the nonlinear diffusion coefficient \tilde{p} at every node on the full grid. Following Glover et al. (2009, eq. 2.10), \lambda is constant and there is no explicit S factor in the denominator:

function p_nl = nlbs_diffcoeff(x, sigma, lambda, Vss)
  % Base coefficient (linear BS): ½σ²S²
  p_nl = 0.5 * sigma^2 .* x.^2;
  if lambda == 0, return; end

  % Denominator: 1 - λ·Vss  (constant λ, Glover et al. 2009 eq. 2.10)
  % Note: no S factor here — contrast with Frey (1998): 1 - λ̂·S·Vss
  denom = 1 - lambda .* Vss;

  % Safety clamp: prevent division by zero near the singularity Vss = 1/λ
  denom_min = 0.01;
  bad = abs(denom) < denom_min;
  if any(bad)
    denom(bad) = sign(denom(bad) + eps) * denom_min;
  end

  % Nonlinear diffusion coefficient (eq. 2.10)
  p_nl = p_nl ./ (denom.^2);
end

When lambda = 0 the function returns the standard linear coefficient immediately. The denominator clamp prevents division by zero when \lambda V_{SS} \to 1.

5.2 nlbs_picard_step.m — the iteration driver

This function encapsulates the combined Picard + penalty loop for a single Crank–Nicolson timestep. Its signature is:

function [uvct1, Vss_new, picard_iters] = nlbs_picard_step(...
    x, Nx, uvct0, rhs_linear, A2, A1, A0, A2ext, ...
    pvct_base, qvct, rvct, ht, theta, ...
    sigma, lambda_fp, bc_mode, ...
    picard_maxiter, picard_tol, ...
    use_american, payoff, american_tol, Large)

The function receives the RHS vector already assembled in parab_nl.m (built with the linear L from the previous timestep), and iteratively updates only the implicit-side operator using the lagged V_{SS}. The penalty matrix P is also lagged and updated each iteration, collapsing what was previously a nested double loop into a single combined loop.

5.3 parab_nl.m — modified main script

This is a lightly modified copy of parab.m with three additions:

Addition 1: the model switch.

% Full-feedback nonlinear BS switch (Glover et al. 2009, eq. 2.10)
use_nlbs  = true;    % false → identical to parab.m
lambda_fp = 0.01;    % constant liquidity parameter; 0 → linear BS
picard_maxiter = 10;
picard_tol     = 1e-6;

Addition 2: the dispatch block. Inside the timestepping loop, the single line uvct1 = AA\rhs (or the penalty_solver call) is replaced by a two-branch conditional:

if use_nlbs && lambda_fp ~= 0
  [uvct1, ~, picard_count] = nlbs_picard_step( ...
      x, Nx, uvct0, rhs, A2, A1, A0, A2ext, ...
      pvct, qvct, rvct, ht_use, theta, ...
      sigma, lambda_fp, bc_mode, ...
      picard_maxiter, picard_tol, ...
      use_american, payoff, american_tol, Large);
  nit(kk, nts(kk)) = picard_count;
else
  if use_american
    [uvct1, iter] = penalty_solver(...);
  else
    uvct1 = AA \ rhs;
  end
end

Addition 3. The pvct comment is updated to clarify it is the base (linear) coefficient; the nonlinear version is computed inside nlbs_picard_step.m.

6. Three bugs found and fixed during integration

Bug 1 — Variable name clash: lambda overwritten by adaptive timestep loop

Bug 2 — Wrong argument passed to penalty_solver inside Picard loop

Bug 3 — Operator precedence error in use_halfstep condition

7. Convergence results

Pointwise errors at S = 100 for the American put under the Glover et al. (2009) full-feedback model with \lambda = 0.01, Crank–Nicolson timestepping, and a nonuniform spatial grid. OrdV, OrdD, OrdG are the computed convergence orders for the option value, delta, and gamma respectively.

8. Summary of changes

9. Status and next steps

• ✓ PDE derivation and alignment with Glover–Duck–Newton (2009, eq. 2.10)
• ✓ nlbs_diffcoeff.m — constant-\lambda nonlinear coefficient
• ✓ nlbs_picard_step.m — combined Picard + penalty iteration
• ✓ parab_nl.m — integrated switch with original solver
• ✓ Three integration bugs identified and resolved
• ✓ Grid convergence study confirming O(h²) in value and delta
• ◯ Quantitative comparison: option prices and early-exercise boundaries for \lambda \in \{0, 0.01, 0.05, 0.10\}
• ◯ Investigation of the V_{SS} = 1/\lambda singularity regime (Glover et al. 2009, Section 5)
• ◯ Full write-up

References

• GDN09 K. J. Glover, P. W. Duck, and D. P. Newton, "On nonlinear models of markets with finite liquidity: Some cautionary notes," SIAM Journal on Applied Mathematics, vol. 70, no. 6, pp. 3252–3271, 2009.
• SW00 P. J. Schönbucher and P. Wilmott, "The feedback effect of hedging in illiquid markets," SIAM Journal of Applied Mathematics, vol. 61, pp. 232–272, 2000.
• FP02 R. Frey and P. Patie, "Risk management for derivatives in illiquid markets: A simulation study," in Advances in Finance and Stochastics, Springer, 2002, pp. 137–159. (Uses asset-dependent \lambda(S)=\hat\lambda S, giving denominator 1-\hat\lambda S\,V_{SS}.)
• Fr98 R. Frey, "Perfect option replication for a large trader," Finance and Stochastics, vol. 2, pp. 115–142, 1998.
• KV13 M. N. Koleva and L. G. Vulkov, "A second-order positivity preserving numerical method for Gamma equation," Applied Mathematics and Computation, vol. 220, pp. 722–734, 2013.
• BS73 F. Black and M. Scholes, "The pricing of options and corporate liabilities," Journal of Political Economy, vol. 81, no. 3, pp. 637–654, 1973.
• CW03 C. C. Christara and K. S. Wu, "Penalty methods for American options," working paper, University of Toronto, 2003.