Calibrate a Heston Smile with Java Monte Carlo

import com.nablatensor.engine.ADouble;
import com.nablatensor.engine.AadRecorder;
import com.nablatensor.quant.Calibrator;
import com.nablatensor.quant.MultiOutput;
import java.util.LinkedHashMap;
import java.util.Map;

/**
 * nablatensor.com/learn, chapter 7.2: the same recording as the real {@code
 * HestonCalibrationTest} -- a full-truncation Euler path to maturity, three
 * strikes' payoffs returned as a named residual map, fit by Levenberg-
 * Marquardt with common random numbers between the target-generating run and
 * every calibration iteration -- rewritten on {@code .on("cpu")} instead of
 * the test's {@code cpu-jit}, at far fewer steps/scenarios so the scalar
 * interpreter stays interactive. {@code MultiOutput.Measures} and {@code
 * Calibrator.leastSquares} both take a plain {@code AadRecorder} lambda, so
 * there is no market-record reflection anywhere in this file.
 */
public final class HestonCalibrationRiskStudio {
  private static final double S0 = 100, R = 0.02, T = 1;
  private static final int STEPS = 8;
  private static final long SCEN = 500L;
  private static final long SEED = 42L;
  private static final double KAPPA = 1.5, THETA = 0.045;
  private static final double[] STRIKES = {
    85, 95, 100, 105, 115
  };
  private HestonCalibrationRiskStudio() {}

  public static void main(String[] args) {
    double trueV0 = 0.05, trueXi = 0.55, trueRho = -0.6;
    double[] target = new double[STRIKES.length];
    com.nablatensor.quant.transform.HestonCf truth = com.nablatensor.quant.transform.HestonCf.of()
        .rate(R)
        .v0(trueV0)
        .kappa(KAPPA)
        .theta(THETA)
        .xi(trueXi)
        .rho(trueRho)
        .build();
    for (int i = 0; i < STRIKES.length; i++) target[i] = com.nablatensor.quant.transform.CosMethod.price(truth,
        com.nablatensor.quant.OptionTypeEnum.CALL, S0, STRIKES[i], R, T);
    long t0 = System.nanoTime();
    Calibrator.Result res = Calibrator.leastSquares(rec -> {
      ADouble v0 = rec.input("v0", 0.03); ADouble kappa = rec.input("kappa", 1.0); ADouble theta = rec.input("theta",
          0.03); ADouble xi = rec.input("xi", 0.35); ADouble rho = rec.input("rho", -0.2); Map<String,
          ADouble> m = hestonMeasures(rec, v0, kappa, theta, xi, rho); Map<String, ADouble> residuals = new LinkedHashMap<>(); for (int i = 0; i < STRIKES.length; i++) {
        residuals.put("k" + i, m.get("k" + i)
            .sub(target[i]));
      }
      return residuals;
    })
        .parameter("v0", 0.03, 1e-3, 0.5)
        .parameter("kappa", 1.0, 1e-2, 20.0)
        .parameter("theta", 0.03, 1e-4, 1.0)
        .parameter("xi", 0.35, 1e-2, 3.0)
        .parameter("rho", -0.2, -0.98, 0.5)
        .scenarios(SCEN)
        .seed(SEED)
        .maxIterations(12)
        .tolerance(1e-10)
        .on("cpu")
        .solve();
    double ms = (System.nanoTime() - t0) / 1e6;
    double v0 = res.parameters()
        .get("v0"), kappa = res.parameters()
        .get("kappa"), theta = res.parameters()
        .get("theta"), xi = res.parameters()
        .get("xi"), rho = res.parameters()
        .get("rho");
    System.out.println("RESULT|" + v0 + "|" + kappa + "|" + theta + "|" + xi + "|" + rho + "|" + res.objective()
        + "|" + res.iterations() + "|" + res.converged() + "|" + ms);
    try (MultiOutput mo = MultiOutput.of(rec -> measures(rec, v0, kappa, theta, xi, rho))
        .threads(1)
        .on("cpu")
        .build()) {
      MultiOutput.Result fitted = mo.run(SCEN, SEED);
      for (int i = 0; i < STRIKES.length; i++) System.out.println("ROW|" + STRIKES[i] + "|" + target[i]
          + "|" + fitted.value("k" + i));
    }
  }

  private static Map<String, ADouble> measures(AadRecorder rec, double v0, double kappa, double theta,
      double xi, double rho) {
    return hestonMeasures(rec, rec.constant(v0), rec.constant(kappa), rec.constant(theta), rec.constant(xi),
        rec.constant(rho));
  }

  private static Map<String, ADouble> hestonMeasures(AadRecorder rec, ADouble v0, ADouble kappa,
      ADouble theta, ADouble xi, ADouble rho) {
    ADouble terminal = hestonTerminal(rec, v0, kappa, theta, xi, rho);
    ADouble disc = rec.constant(Math.exp(-R * T));
    Map<String, ADouble> m = new LinkedHashMap<>();
    for (int i = 0; i < STRIKES.length; i++) {
      m.put("k" + i, terminal.sub(STRIKES[i])
          .max(0.0)
          .mul(disc));
    }
    return m;
  }
  /** Full-truncation Euler Heston to T; returns the terminal spot. */
  private static ADouble hestonTerminal(AadRecorder rec, ADouble v0, ADouble kappa, ADouble theta,
      ADouble xi, ADouble rho) {
    double dt = T / STEPS, sqrtDt = Math.sqrt(dt);
    ADouble rhoBar = rho.mul(rho)
        .neg()
        .add(1.0)
        .sqrt();
    ADouble s = rec.constant(S0);
    ADouble v = v0;
    for (int t = 0; t < STEPS; t++) {
      ADouble z1 = rec.randn();
      ADouble z2 = rho.mul(z1)
          .add(rhoBar.mul(rec.randn()));
      ADouble vPlus = v.max(0.0);
      ADouble sqrtV = vPlus.sqrt();
      v = v.add(kappa.mul(theta.sub(vPlus))
          .mul(dt))
          .add(xi.mul(sqrtV)
          .mul(sqrtDt)
          .mul(z2));
      s = s.mul(rec.constant(R)
          .sub(vPlus.mul(0.5))
          .mul(dt)
          .add(sqrtV.mul(sqrtDt)
          .mul(z1))
          .exp());
    }
    return s;
  }
}