Calibrate Fast, Without Recording a Tape

import com.nablatensor.quant.OptionTypeEnum;
import com.nablatensor.quant.analytic.GeneralizedBsm;
import com.nablatensor.quant.transform.BsmCf;
import com.nablatensor.quant.transform.CosMethod;
import com.nablatensor.quant.transform.HestonCf;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.Locale;
import java.util.function.ToDoubleFunction;

/**
 * nablatensor.com/learn, chapter 7.3: the real {@code CosCalibrationShowcase}'s
 * setup -- a Heston surface priced by the COS method, then fit back with
 * Nelder-Mead -- at a smaller quote grid and iteration cap. {@code
 * HestonCf}/{@code CosMethod} never touch {@code ADouble}, so there's no
 * reflection concern here at all; the reason for the smaller cap is
 * different: {@code HestonCosCalibrator.calibrate} hardcodes 2000 iterations
 * with no way to lower it, and the real class's own {@code NelderMead} is
 * package-private, so this file inlines the same reflect/expand/contract/
 * shrink simplex (copied verbatim from {@code NelderMead.minimise}) with a
 * cap this page's interpreter can clear quickly.
 */
public final class CosCalibrationRiskStudio {
  private CosCalibrationRiskStudio() {}

  public static void main(String[] args) {
    double s = 100;
    double r = 0.02;
    BsmCf bsm = BsmCf.of()
        .rate(r)
        .vol(0.2)
        .build();
    double maxErr = 0.0;
    for (double k : new double[] {
      70, 85, 100, 115, 130
    }) {
      double cos = CosMethod.price(bsm, OptionTypeEnum.CALL, s, k, r, 1.0);
      double closed = GeneralizedBsm.of()
          .type(OptionTypeEnum.CALL)
          .spot(s)
          .strike(k)
          .maturity(1.0)
          .rate(r)
          .dividend(0.0)
          .vol(0.2)
          .build()
          .price();
      maxErr = Math.max(maxErr, Math.abs(cos - closed));
    }
    HestonCf truth = heston(r, 0.045, 2, 0.05, 0.55, -0.65);
    double[] maturities = {
      0.5, 1.0
    };
    double[] strikes = {
      90, 100, 110
    };
    List<double[]> quotes = new ArrayList<>();
    // {strike, maturity, price}
    double bump = 1.0;
    for (double t : maturities) {
      for (double k : strikes) {
        double px = CosMethod.price(truth, OptionTypeEnum.CALL, s, k, r, t);
        quotes.add(new double[] {
          k, t, px * bump
        });
        bump = bump == 1.0 ? 1.004 : 1.0;
      }
    }
    ToDoubleFunction<double[]> sse = x -> {
      HestonCf cf = heston(r, x[0], x[1], x[2], x[3], x[4]);
      double sum = 0.0;
      for (double[] q : quotes) {
        double model = CosMethod.price(cf, OptionTypeEnum.CALL, s, q[0], r, q[1]);
        double d = model - q[2];
        sum += d * d;
      }
      return sum;
    };
    double[] start = {
      0.04, 1.0, 0.04, 0.3, -0.3
    };
    double[] step = {
      0.01, 0.5, 0.01, 0.05, 0.05
    };
    double[] lo = {
      1e-4, 1e-2, 1e-4, 1e-3, -0.999
    };
    double[] hi = {
      1.0, 20.0, 1.0, 5.0, 0.999
    };
    long t0 = System.nanoTime();
    Result fit = minimise(sse, start, step, lo, hi, 80, 1e-14);
    double ms = (System.nanoTime() - t0) / 1e6;
    double rmse = Math.sqrt(sse.applyAsDouble(fit.x) / quotes.size());
    System.out.println("RESULT|" + fit.x[0] + "|" + fit.x[1] + "|" + fit.x[2] + "|" + fit.x[3] + "|"
        + fit.x[4] + "|" + rmse + "|" + fit.iterations + "|" + fit.converged + "|" + ms);
    HestonCf fitted = heston(r, fit.x[0], fit.x[1], fit.x[2], fit.x[3], fit.x[4]);
    for (double[] q : quotes) System.out.println("ROW|" + q[0] + "|" + q[1] + "|" + q[2] + "|" + CosMethod.price(fitted,
        OptionTypeEnum.CALL, s, q[0], r, q[1]));
  }

  private record Result(double[] x, int iterations, boolean converged) {}

  private static HestonCf heston(double r, double v0, double kappa, double theta, double xi, double rho) {
    return HestonCf.of()
        .rate(r)
        .v0(v0)
        .kappa(kappa)
        .theta(theta)
        .xi(xi)
        .rho(rho)
        .build();
  }
  /** {@code NelderMead.minimise}, inlined -- the real class is package-private. */
  private static Result minimise(ToDoubleFunction<double[]> f, double[] start, double[] step, double[] lo,
      double[] hi, int maxIter, double tol) {
    int n = start.length;
    double[][] p = new double[n + 1][];
    double[] fv = new double[n + 1];
    p[0] = clamp(start.clone(), lo, hi);
    fv[0] = f.applyAsDouble(p[0]);
    for (int i = 0; i < n; i++) {
      double[] v = start.clone();
      v[i] += step[i];
      p[i + 1] = clamp(v, lo, hi);
      fv[i + 1] = f.applyAsDouble(p[i + 1]);
    }
    int iter = 0;
    for (; iter < maxIter; iter++) {
      order(p, fv);
      if (Math.abs(fv[n] - fv[0]) <= tol * (Math.abs(fv[0]) + tol)) {
        break;
      }
      double[] centroid = new double[n];
      for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
          centroid[j] += p[i][j] / n;
        }
      }
      double[] worst = p[n];
      double[] reflected = combine(centroid, worst, 1.0, lo, hi);
      double fr = f.applyAsDouble(reflected);
      if (fr < fv[0]) {
        double[] expanded = combine(centroid, worst, 2.0, lo, hi);
        double fe = f.applyAsDouble(expanded);
        replaceWorst(p, fv, fe < fr ? expanded : reflected, Math.min(fe, fr));
      } else if (fr < fv[n - 1]) {
        replaceWorst(p, fv, reflected, fr);
      } else {
        double[] contracted = combine(centroid, worst, -0.5, lo, hi);
        double fc = f.applyAsDouble(contracted);
        if (fc < fv[n]) {
          replaceWorst(p, fv, contracted, fc);
        } else {
          for (int i = 1; i <= n; i++) {
            for (int j = 0; j < n; j++) {
              p[i][j] = p[0][j] + 0.5 * (p[i][j] - p[0][j]);
            }
            p[i] = clamp(p[i], lo, hi);
            fv[i] = f.applyAsDouble(p[i]);
          }
        }
      }
    }
    order(p, fv);
    return new Result(p[0].clone(), iter, iter < maxIter);
  }

  private static double[] combine(double[] centroid, double[] worst, double coef, double[] lo, double[] hi) {
    double[] out = new double[centroid.length];
    for (int i = 0; i < out.length; i++) {
      out[i] = centroid[i] + coef * (centroid[i] - worst[i]);
    }
    return clamp(out, lo, hi);
  }

  private static void replaceWorst(double[][] p, double[] fv, double[] point, double value) {
    p[p.length - 1] = point;
    fv[fv.length - 1] = value;
  }

  private static void order(double[][] p, double[] fv) {
    Integer[] idx = new Integer[fv.length];
    for (int i = 0; i < idx.length; i++) {
      idx[i] = i;
    }
    Arrays.sort(idx, (x, y) -> Double.compare(fv[x], fv[y]));
    double[][] np = new double[p.length][];
    double[] nfv = new double[fv.length];
    for (int i = 0; i < idx.length; i++) {
      np[i] = p[idx[i]];
      nfv[i] = fv[idx[i]];
    }
    System.arraycopy(np, 0, p, 0, p.length);
    System.arraycopy(nfv, 0, fv, 0, fv.length);
  }

  private static double[] clamp(double[] x, double[] lo, double[] hi) {
    for (int i = 0; i < x.length; i++) {
      x[i] = Math.max(lo[i], Math.min(hi[i], x[i]));
    }
    return x;
  }
}