feat(american-option): log-spot Dynamic Chebyshev variant (Stage F1 front-fixing)

docs/docs/american-option-dynamic-chebyshev.md

@@ -831,8 +831,22 @@ Gamma query onto the **clustered edge of the upper piece**, where a Chebyshev
 second-derivative differentiation matrix is most ill-conditioned (a measured
 ~4000x edge amplification versus the interior). The split even produced a
 non-physical *negative* Gamma one tick above the knot. Raising the node count
-makes both variants worse, because the wide-spot Gauss-Hermite quadrature loses
-finite values at high `n`.
+does not rescue it — the wide `[5, 250]` grid becomes ill-conditioned at high
+`n` and the build returns non-finite values (`n = 321` throws).
+
+The wiggle is concrete (`B0 = 81.868`): the split-piece Gamma rings across the
+knot edge while the global interpolant climbs smoothly.
+
+| spot | global Γ | split Γ |
+| ---: | ---: | ---: |
+| 81.85 | 0.000000 | 0.000000 |
+| 81.90 | 0.025163 | **-0.000404** |
+| 82.00 | 0.025696 | 0.019329 |
+| 82.25 | 0.026978 | 0.011920 |
+| 82.60 | 0.028640 | 0.036449 |
+
+Negative one tick above the knot, a trough at `82.25`, then an overshoot past the
+global by `82.60` — endpoint ringing, not signal.
 
 What *does* help, cheaply, is seeding the terminal backward step with the exact
 one-period European Black-Scholes value (the `ClosedFormTerminalStep` option)
@@ -848,6 +862,25 @@ endpoint sitting on the query (Company, Egorova & Jodar 2014). That is a solver
 re-architecture rather than a knot tweak, and is the direction for a future
 boundary-aware variant.
 
+### A log-spot variant recovers the boundary Gamma
+
+The first step of that front-fixing direction is already informative. Interpolating the continuation in
+`x = log(S)` (a `LogSpot` build option) makes the Gauss–Hermite transition additive, so the images stay
+bounded, and the grid narrow and uniform in `x`, hence far better conditioned than the wide linear
+`[5, 250]` grid. Measured spot-82 Gamma error against the oracle (`0.033689`):
+
+| nodes `n` | linear Γ(82) error | log-spot Γ(82) error |
+| ---: | ---: | ---: |
+| 81 | 0.007993 | 0.009807 |
+| 161 | 0.033689 (linear: spot 82 collapses into the exercise region) | 0.005922 |
+| 321 | non-finite (linear build throws) | 0.002755 |
+
+The log-spot error **falls monotonically with `n`** — the boundary Gamma was resolution-limited, not
+intrinsic — and the log-spot build stays finite at node counts where the linear grid throws. At low `n`
+the log-spot nodes cluster at the domain ends rather than at the boundary, so it trails the linear grid
+there; the gain appears once the grid is refined, and where the linear grid can no longer run. This
+recovers most of the boundary Gamma without the full boundary-tracking (Landau) re-architecture.
+
 ## Accuracy and Speed
 
 The aggregate comparison is useful only after reading the per-candidate

examples/AmericanOptionDynamicChebyshev/DynamicChebyshev.cs

@@ -10,14 +10,11 @@ public sealed record DynamicChebyshevSettings(
     double SpotUpper = 250.0,
     int QuadratureOrder = 8,
     bool ClosedFormTerminalStep = false,
-    // Experimental, and a DOCUMENTED NEGATIVE RESULT: splitting the smooth continuation at the
-    // exercise boundary does not improve Greeks. The kink lives in the price max(payoff, C), which
-    // is applied exactly outside the interpolant; the continuation C is smooth at the boundary, so a
-    // knot there resolves no singularity and only relocates the Gamma query onto the ill-conditioned
-    // Chebyshev piece edge (it yields a non-physical negative Gamma one tick above the knot). Kept,
-    // default-off, for reproducibility; superseded by the planned front-fixing variant. See the
-    // "Why the boundary Gamma is weak" section of the American-option case study.
-    bool BoundarySplit = false);
+    // Stage F1 front-fixing: interpolate the continuation in x = log(S) instead of linear S. The GBM
+    // transition is additive in x, so the Gauss-Hermite images stay bounded; the narrow uniform-in-x
+    // grid is also far better conditioned at high node counts, curing the high-n non-finite build
+    // failure seen on the wide linear [5,250] grid. Default off; the OFF path is bit-identical.
+    bool LogSpot = false);
 
 public sealed record DynamicChebyshevResult(
     double Price,
@@ -156,56 +153,90 @@ public DynamicChebyshevAmericanOptionModel Build(
 
         var stopwatch = Stopwatch.StartNew();
         int buildEvaluations = 0;
-        Func<double, double> nextValue = spot => Payoff(request, spot);
+        // The terminal next-step value is the payoff. In the log-spot frame the value functions are
+        // indexed by x = log(S), so the seed and every step are expressed in x; otherwise in S.
+        Func<double, double> nextValue = settings.LogSpot
+            ? x => Payoff(request, Math.Exp(x))
+            : spot => Payoff(request, spot);
         ChebyshevApproximation? firstContinuation = null;
-        Func<double, double>? firstContinuationFunction = null;
 
         for (int step = settings.ExerciseSteps - 1; step >= 0; step--)
         {
             Func<double, double> valueAtNextStep = nextValue;
             bool closedFormTerminal = settings.ClosedFormTerminalStep
                 && step == settings.ExerciseSteps - 1
                 && request.Right == VanillaOptionRight.Put;
-            Func<double, double> continuationFunction = closedFormTerminal
-                ? spot => EuropeanPutPrice(
-                    spot, request.Strike, dt, request.RiskFreeRate, request.Volatility, request.DividendYield)
-                : spot => ContinuationValue(spot, drift, diffusion, discount, valueAtNextStep);
-            ChebyshevApproximation continuation = BuildContinuationApproximation(
-                request,
-                settings,
-                continuationFunction,
-                maxDerivativeOrder: 2);
-
-            buildEvaluations += continuation.NEvaluations;
-            nextValue = spot => Math.Max(
-                Payoff(request, spot),
-                EvaluateApproximation(continuation, spot, settings, derivativeOrder: 0));
+
+            ChebyshevApproximation continuation;
+            Func<double, double> continuationFunction;
+            if (settings.LogSpot)
+            {
+                // x = log(S): additive (bounded) transition, e^x payoff; the recursion runs in x.
+                continuationFunction = closedFormTerminal
+                    ? x => EuropeanPutPrice(
+                        Math.Exp(x), request.Strike, dt, request.RiskFreeRate, request.Volatility, request.DividendYield)
+                    : x => LogSpotContinuationValue(x, drift, diffusion, discount, valueAtNextStep);
+                continuation = BuildLogSpotContinuationApproximation(
+                    settings, continuationFunction, maxDerivativeOrder: 2);
+
+                buildEvaluations += continuation.NEvaluations;
+                nextValue = x => Math.Max(
+                    Payoff(request, Math.Exp(x)),
+                    EvaluateLogSpotApproximation(continuation, x, settings, derivativeOrder: 0));
+            }
+            else
+            {
+                continuationFunction = closedFormTerminal
+                    ? spot => EuropeanPutPrice(
+                        spot, request.Strike, dt, request.RiskFreeRate, request.Volatility, request.DividendYield)
+                    : spot => ContinuationValue(spot, drift, diffusion, discount, valueAtNextStep);
+                continuation = BuildContinuationApproximation(
+                    request,
+                    settings,
+                    continuationFunction,
+                    maxDerivativeOrder: 2);
+
+                buildEvaluations += continuation.NEvaluations;
+                nextValue = spot => Math.Max(
+                    Payoff(request, spot),
+                    EvaluateApproximation(continuation, spot, settings, derivativeOrder: 0));
+            }
 
             if (step == 0)
             {
                 firstContinuation = continuation;
-                firstContinuationFunction = continuationFunction;
             }
         }
 
         stopwatch.Stop();
 
         Debug.Assert(firstContinuation is not null);
-        Debug.Assert(firstContinuationFunction is not null);
         ChebyshevApproximation firstApprox = firstContinuation!;
-        Func<double, double> firstFunc = firstContinuationFunction!;
 
         Func<double, int, double> continuationCurve;
-        if (settings.BoundarySplit && request.Right == VanillaOptionRight.Put)
+        if (settings.LogSpot)
         {
-            double GlobalContinuation(double spot) =>
-                EvaluateApproximation(firstApprox, spot, settings, derivativeOrder: 0);
-            double PayoffAt(double spot) => Payoff(request, spot);
-            double boundary = FindExerciseBoundary(GlobalContinuation, PayoffAt, settings.SpotLower, request.Strike);
-
-            ChebyshevSpline spline = BuildContinuationSpline(settings, firstFunc, boundary);
-            buildEvaluations += spline.NumPieces * settings.SpotNodeCount;
-            continuationCurve = (spot, order) => EvaluateSpline(spline, spot, settings, order, boundary);
+            // The first continuation is interpolated in x = log(S). Convert S->x once, read the
+            // spectral x-derivatives, and apply the chain rule: Delta = u'(x)/S, Gamma = (u''-u')/S^2.
+            ChebyshevApproximation firstLog = firstApprox;
+            continuationCurve = (spot, order) =>
+            {
+                double x = Math.Log(spot);
+                double u0 = EvaluateLogSpotApproximation(firstLog, x, settings, derivativeOrder: 0);
+                if (order == 0)
+                {
+                    return u0;
+                }
+
+                double u1 = EvaluateLogSpotApproximation(firstLog, x, settings, derivativeOrder: 1);
+                if (order == 1)
+                {
+                    return u1 / spot;
+                }
+
+                double u2 = EvaluateLogSpotApproximation(firstLog, x, settings, derivativeOrder: 2);
+                return (u2 - u1) / (spot * spot);
+            };
         }
         else
         {
@@ -266,6 +297,52 @@ internal static double EvaluateApproximation(
         return approximation.VectorizedEval([clamped], [derivativeOrder]);
     }
 
+    private static ChebyshevApproximation BuildLogSpotContinuationApproximation(
+        DynamicChebyshevSettings settings,
+        Func<double, double> continuation,
+        int maxDerivativeOrder)
+    {
+        // point[0] is x = log(S); the continuation function is already expressed in x.
+        double Function(double[] point, object? _) => continuation(point[0]);
+
+        var approximation = new ChebyshevApproximation(
+            Function,
+            numDimensions: 1,
+            domain: [[Math.Log(settings.SpotLower), Math.Log(settings.SpotUpper)]],
+            nNodes: [settings.SpotNodeCount],
+            maxDerivativeOrder: maxDerivativeOrder);
+        approximation.Build(verbose: false);
+        return approximation;
+    }
+
+    private static double LogSpotContinuationValue(
+        double x,
+        double drift,
+        double diffusion,
+        double discount,
+        Func<double, double> nextValue)
+    {
+        double sum = 0.0;
+        for (int i = 0; i < HermiteNodes8.Length; i++)
+        {
+            // Additive in x (bounded), unlike the multiplicative linear image spot * exp(...).
+            double nextX = x + drift + Math.Sqrt(2.0) * diffusion * HermiteNodes8[i];
+            sum += HermiteWeights8[i] * nextValue(nextX);
+        }
+
+        return discount * sum / Math.Sqrt(Math.PI);
+    }
+
+    internal static double EvaluateLogSpotApproximation(
+        ChebyshevApproximation approximation,
+        double x,
+        DynamicChebyshevSettings settings,
+        int derivativeOrder)
+    {
+        double clamped = Math.Clamp(x, Math.Log(settings.SpotLower), Math.Log(settings.SpotUpper));
+        return approximation.VectorizedEval([clamped], [derivativeOrder]);
+    }
+
     internal static double Payoff(AmericanOptionRequest request, double spot)
     {
         return request.Right switch
@@ -339,41 +416,6 @@ internal static double FindExerciseBoundary(
         return MathNet.Numerics.RootFinding.Brent.FindRoot(Gap, lo, hi, accuracy: 1e-8, maxIterations: 100);
     }
 
-    private static ChebyshevSpline BuildContinuationSpline(
-        DynamicChebyshevSettings settings, Func<double, double> continuation, double knot)
-    {
-        double Function(double[] point, object? _) => continuation(point[0]);
-
-        var spline = new ChebyshevSpline(
-            Function,
-            numDimensions: 1,
-            domain: [[settings.SpotLower, settings.SpotUpper]],
-            nNodes: [settings.SpotNodeCount],
-            knots: [[knot]],
-            maxDerivativeOrder: 2);
-        spline.Build(verbose: false);
-        return spline;
-    }
-
-    private static double EvaluateSpline(
-        ChebyshevSpline spline,
-        double spot,
-        DynamicChebyshevSettings settings,
-        int derivativeOrder,
-        double knot)
-    {
-        double clamped = Math.Clamp(spot, settings.SpotLower, settings.SpotUpper);
-
-        // Derivatives are undefined exactly at the knot (adjacent pieces disagree there); nudge onto
-        // the continuation side so Greeks remain queryable right at the boundary.
-        if (derivativeOrder > 0 && Math.Abs(clamped - knot) < 1e-9)
-        {
-            clamped = knot + 1e-7;
-        }
-
-        return spline.Eval([clamped], [derivativeOrder]);
-    }
-
     private static void Validate(
         AmericanOptionRequest request,
         DynamicChebyshevSettings settings)

tests/ChebyshevSharp.Tests/Finance/AmericanOptionBoundaryFinderTests.cs

@@ -0,0 +1,41 @@
+using System.Reflection;
+using AmericanOptionDynamicChebyshev;
+
+namespace ChebyshevSharp.Tests.Finance;
+
+/// <summary>
+/// The exercise-boundary finder: B (where payoff = continuation) located by Math.NET Brent
+/// root-finding. The boundary itself is a useful diagnostic and the basis for any future
+/// boundary-aware variant. (The earlier boundary-<i>split</i> experiment was retired as a documented
+/// negative result — see the case study — once the log-spot variant superseded it.)
+/// </summary>
+public sealed class AmericanOptionBoundaryFinderTests
+{
+    [Fact]
+    public void Exercise_boundary_is_located_near_the_diagnostic_value()
+    {
+        AmericanOptionRequest request = AmericanOptionScenarios.StandardPut();
+        var settings = new DynamicChebyshevSettings(80, 81, 5.0, 250.0, 8);
+        DynamicChebyshevAmericanOptionModel model =
+            new DynamicChebyshevAmericanOptionPricer().Build(request, settings);
+
+        double boundary = InvokeFindExerciseBoundary(model, lo: 5.0, hi: request.Strike);
+
+        // The diagnostics mode reports the first-step exercise boundary around spot 81.86.
+        Assert.InRange(boundary, 80.0, 84.0);
+    }
+
+    private static double InvokeFindExerciseBoundary(
+        DynamicChebyshevAmericanOptionModel model, double lo, double hi)
+    {
+        MethodInfo method = typeof(DynamicChebyshevAmericanOptionPricer)
+            .GetMethod(
+                "FindExerciseBoundary",
+                BindingFlags.Static | BindingFlags.NonPublic,
+                binder: null,
+                types: [typeof(DynamicChebyshevAmericanOptionModel), typeof(double), typeof(double)],
+                modifiers: null)
+            ?? throw new InvalidOperationException("FindExerciseBoundary is not implemented.");
+        return (double)method.Invoke(null, [model, lo, hi])!;
+    }
+}