Add Merton jump-diffusion and fix variance/percentile estimators

Parity with Axis2/C commit 6e577eb + review fixes fb40242, 0b281fd.

Merton (1976) jump-diffusion model:
  - New request fields: model ("gbm"|"merton"), jumpIntensity,
    jumpMean, jumpVol — all with defaults matching C implementation.
  - Drift correction: (μ − σ²/2 − λk)·dt preserves E[S(T)].
  - Bernoulli jump process with lambda·dt validation (> 0.1 rejected).
  - Response includes "model" field echoing which model was used.

Numerical stability (from Gemini quant review):
  - Replaced one-pass variance (sumSq/N − mean²) with two-pass
    algorithm. The one-pass formula suffers catastrophic cancellation
    when stdDev << mean.
  - Fixed VaR percentile indexing: ceil(p·N) − 1 instead of
    floor(p·N). The floor estimator selects one observation too far
    from the tail, systematically understating VaR.

All existing tests pass.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: robertlazarski

modules/samples/userguide/src/userguide/springbootdemo-tomcat11/src/main/java/userguide/springboot/webservices/FinancialBenchmarkService.java

@@ -365,15 +365,43 @@ public MonteCarloResponse monteCarlo(MonteCarloRequest request)
             throw JsonRpcFaultException.validationError("volatility must be >= 0.");
         }
 
-        // ── Pre-computed GBM constants ────────────────────────────────────────
+        // ── Model selection ────────────────────────────────────────────────────
+        boolean isMerton = "merton".equalsIgnoreCase(request.getModel());
+        double jumpIntensity = request.getJumpIntensity();
+        double jumpMean = request.getJumpMean();
+        double jumpVol = request.getJumpVol();
+
+        if (isMerton && jumpVol < 0.0) {
+            throw JsonRpcFaultException.validationError("jumpVol must be >= 0.");
+        }
+
+        // ── Pre-computed constants ────────────────────────────────────────────
         // dt         — length of one time step in years (e.g., 1/252 for one
         //              trading day when nPeriodsPerYear = 252)
-        // drift      — (μ − σ²/2)·dt.  The (−σ²/2) term is the Itô correction
-        //              that keeps E[S(T)] = S(0)·exp(μ·T).  Do not drop it.
+        // drift      — (μ − σ²/2 − λk)·dt for Merton, (μ − σ²/2)·dt for GBM.
+        //              The −λk term compensates for the expected jump so that
+        //              E[S(T)] = S(0)·exp(μ·T) regardless of jump parameters.
         // volSqrtDt  — σ·√dt.  Scales each standard-normal shock by the
         //              one-step standard deviation.
         double dt = 1.0 / npy;
-        double drift = (mu - 0.5 * sigma * sigma) * dt;
+        double jumpCompensation = 0.0;
+        double jumpLambdaDt = 0.0;
+        if (isMerton) {
+            double k = Math.exp(jumpMean + 0.5 * jumpVol * jumpVol) - 1.0;
+            jumpCompensation = jumpIntensity * k;
+            jumpLambdaDt = jumpIntensity * dt;
+            // The Bernoulli approximation for a Poisson process is only valid
+            // when λ·dt << 1. At λ·dt > 0.1 the probability of ≥2 jumps per
+            // step becomes non-negligible; at λ·dt ≥ 1 the trial degenerates
+            // to a deterministic jump every step. Fail fast.
+            if (jumpLambdaDt > 0.1) {
+                throw JsonRpcFaultException.validationError(
+                    "jumpIntensity too high for time step: lambda*dt=" +
+                    String.format("%.4f", jumpLambdaDt) + " > 0.1. " +
+                    "Reduce jumpIntensity or increase nPeriodsPerYear.");
+            }
+        }
+        double drift = (mu - 0.5 * sigma * sigma - jumpCompensation) * dt;
         double volSqrtDt = sigma * Math.sqrt(dt);
 
         // ── PRNG: seeded for reproducibility, unseeded for production ─────────
@@ -386,7 +414,6 @@ public MonteCarloResponse monteCarlo(MonteCarloRequest request)
 
         double[] finalValues = new double[nSims];
         double sumFinal = 0.0;
-        double sumSqFinal = 0.0;
         int profitCount = 0;
         double maxDrawdown = 0.0;
 
@@ -399,7 +426,18 @@ public MonteCarloResponse monteCarlo(MonteCarloRequest request)
 
             for (int period = 0; period < nPeriods; period++) {
                 double z = rng.nextGaussian();
-                value *= Math.exp(drift + volSqrtDt * z);
+                double exponent = drift + volSqrtDt * z;
+
+                // Merton jump component: compound Poisson process.
+                // At each step, a jump occurs with probability λ·dt.
+                // The jump magnitude is log-normal: J = exp(μ_J + σ_J · W)
+                // where W ~ N(0,1) is independent of the diffusion Z.
+                if (isMerton && rng.nextDouble() < jumpLambdaDt) {
+                    double w = rng.nextGaussian();
+                    exponent += jumpMean + jumpVol * w;
+                }
+
+                value *= Math.exp(exponent);
 
                 if (value > peak) {
                     peak = value;
@@ -411,26 +449,35 @@ public MonteCarloResponse monteCarlo(MonteCarloRequest request)
 
             finalValues[sim] = value;
             sumFinal += value;
-            sumSqFinal += value * value;
             if (value > initialValue) profitCount++;
             if (simMaxDrawdown > maxDrawdown) maxDrawdown = simMaxDrawdown;
         }
 
         long elapsedUs = (System.nanoTime() - startNs) / 1_000;
 
         // ── Statistics ────────────────────────────────────────────────────────
+        // Two-pass algorithm for variance: the one-pass formula
+        // (sumSq/N − mean²) suffers catastrophic cancellation when
+        // stdDev << mean (common for low-vol strategies). Two-pass:
+        // compute mean first, then sum squared deviations.
         double mean = sumFinal / nSims;
-        double variance = (sumSqFinal / nSims) - (mean * mean);
-        if (variance < 0.0) variance = 0.0;
+        double sumSqDiff = 0.0;
+        for (int i = 0; i < nSims; i++) {
+            double d = finalValues[i] - mean;
+            sumSqDiff += d * d;
+        }
+        double variance = sumSqDiff / nSims;
 
         // Sort ascending so finalValues[0] is the worst outcome and
         // finalValues[nSims-1] is the best.  VaR and CVaR then become
-        // simple index lookups: the p-th percentile loss corresponds to
-        // finalValues[floor(p * nSims)].
+        // simple index lookups.
         Arrays.sort(finalValues);
 
-        int idx5  = (int)(0.05 * nSims);
-        int idx1  = (int)(0.01 * nSims);
+        // Percentile indexing: ceil(p·N) − 1 selects the k-th order
+        // statistic such that exactly floor(p·N) observations are strictly
+        // below the VaR level. Matches the standard quantile definition.
+        int idx5  = Math.max(0, (int) Math.ceil(0.05 * nSims) - 1);
+        int idx1  = Math.max(0, (int) Math.ceil(0.01 * nSims) - 1);
 
         // Sample median of a sorted array: for odd N take the middle element,
         // for even N average the two central elements.  Using a single index
@@ -467,14 +514,15 @@ public MonteCarloResponse monteCarlo(MonteCarloRequest request)
             for (int pi = 0; pi < maxPct; pi++) {
                 double p = request.getPercentiles()[pi];
                 if (p <= 0.0 || p >= 1.0) continue;
-                int idx = (int)(p * nSims);
+                int idx = Math.max(0, (int) Math.ceil(p * nSims) - 1);
                 if (idx >= nSims) idx = nSims - 1;
                 percentileVars.add(new MonteCarloResponse.PercentileVar(
                     p, initialValue - finalValues[idx]));
             }
         }
 
         MonteCarloResponse response = new MonteCarloResponse();
+        response.setModel(isMerton ? "merton" : "gbm");
         response.setStatus("SUCCESS");
         response.setMeanFinalValue(mean);
         response.setMedianFinalValue(median);

modules/samples/userguide/src/userguide/springbootdemo-tomcat11/src/main/java/userguide/springboot/webservices/MonteCarloRequest.java

@@ -21,10 +21,23 @@
 /**
  * Request for Monte Carlo Value-at-Risk simulation.
  *
- * <p>Simulates portfolio value paths using Geometric Brownian Motion:
+ * <p>Simulates portfolio value paths using one of two stochastic models:
+ *
+ * <p><b>GBM</b> (model="gbm", default): Geometric Brownian Motion.
  * <pre>S(t+dt) = S(t) × exp((μ − σ²/2)·dt + σ·√dt·Z)</pre>
  * where dt = 1/nPeriodsPerYear and Z ~ N(0,1).
  *
+ * <p><b>Merton jump-diffusion</b> (model="merton"): GBM plus a compound
+ * Poisson jump process that captures fat tails and discontinuities
+ * (crashes, gaps) that GBM systematically understates.
+ * <pre>S(t+dt) = S(t) × exp((μ − σ²/2 − λk)·dt + σ·√dt·Z) × J</pre>
+ * where J = exp(μ_J + σ_J·W) with probability λ·dt per step, else J = 1.
+ * Z, W ~ N(0,1) independent. k = exp(μ_J + σ_J²/2) − 1.
+ * The −λk drift correction preserves E[S(T)] = S(0)·exp(μ·T).
+ *
+ * <p>Reference: Merton, R. (1976). "Option pricing when underlying stock
+ * returns are discontinuous." J. Financial Economics, 3(1-2): 125-144.
+ *
  * <p>All fields have defaults so a minimal {@code {}} request body is valid.
  *
  * <h3>Example</h3>
@@ -196,6 +209,38 @@ public class MonteCarloRequest {
      */
     private double[] percentiles = {0.01, 0.05};
 
+    /**
+     * Simulation model: "gbm" (default) or "merton".
+     *
+     * <p>GBM is the textbook baseline — constant vol, continuous paths.
+     * Merton adds a compound Poisson jump process for fat tails.
+     * See class javadoc for the full model equations.
+     */
+    private String model = "gbm";
+
+    /**
+     * Jump intensity (Merton only): expected number of jumps per YEAR
+     * (Poisson λ). Must be &gt;= 0. Default: 1.0. At λ=1, approximately
+     * one jump per simulated year. Higher values increase kurtosis.
+     * Ignored when model="gbm".
+     */
+    private double jumpIntensity = 1.0;
+
+    /**
+     * Mean of the log-normal jump size distribution (Merton only).
+     * Negative values produce downward jumps on average (crashes).
+     * Default: -0.03 (average jump ≈ 3% drop). Ignored when model="gbm".
+     */
+    private double jumpMean = -0.03;
+
+    /**
+     * Volatility of the log-normal jump size (Merton only). Controls how
+     * variable individual jump magnitudes are. Must be &gt;= 0.
+     * Default: 0.05. At jumpVol=0 every jump has exactly size
+     * exp(jumpMean). Ignored when model="gbm".
+     */
+    private double jumpVol = 0.05;
+
     /** Optional identifier echoed in the response for request tracing. */
     private String requestId;
 
@@ -218,6 +263,10 @@ public class MonteCarloRequest {
     public int getNPeriodsPerYear() { return nPeriodsPerYear > 0 ? nPeriodsPerYear : 252; }
     public double[] getPercentiles() { return percentiles != null ? percentiles : new double[]{0.01, 0.05}; }
     public String getRequestId() { return requestId; }
+    public String getModel() { return model != null ? model : "gbm"; }
+    public double getJumpIntensity() { return jumpIntensity >= 0 ? jumpIntensity : 1.0; }
+    public double getJumpMean() { return jumpMean; }
+    public double getJumpVol() { return jumpVol >= 0 ? jumpVol : 0.05; }
 
     // ── setters ──────────────────────────────────────────────────────────────
 
@@ -230,4 +279,8 @@ public class MonteCarloRequest {
     public void setNPeriodsPerYear(int nPeriodsPerYear) { this.nPeriodsPerYear = nPeriodsPerYear; }
     public void setPercentiles(double[] percentiles) { this.percentiles = percentiles; }
     public void setRequestId(String requestId) { this.requestId = requestId; }
+    public void setModel(String model) { this.model = model; }
+    public void setJumpIntensity(double jumpIntensity) { this.jumpIntensity = jumpIntensity; }
+    public void setJumpMean(double jumpMean) { this.jumpMean = jumpMean; }
+    public void setJumpVol(double jumpVol) { this.jumpVol = jumpVol; }
 }

modules/samples/userguide/src/userguide/springbootdemo-tomcat11/src/main/java/userguide/springboot/webservices/MonteCarloResponse.java

@@ -81,6 +81,9 @@ public class MonteCarloResponse {
     /** JVM heap used at response time in MB */
     private long memoryUsedMb;
 
+    /** Model used for this simulation: "gbm" or "merton" */
+    private String model;
+
     /** Echoed from request */
     private String requestId;
 
@@ -161,4 +164,7 @@ public static MonteCarloResponse failed(String errorMessage) {
 
     public String getRequestId() { return requestId; }
     public void setRequestId(String requestId) { this.requestId = requestId; }
+
+    public String getModel() { return model; }
+    public void setModel(String model) { this.model = model; }
 }