VaR Backtesting Calculator: Kupiec & Christoffersen Tests

import com.nablatensor.quant.analytic.Normal;
import java.util.Random;

public final class VarBacktestRiskStudio {
  // Java demo using nablatensor normal-distribution helpers: simulate P&L, then test exception frequency and timing.
  public static void main(String[] args) {
    // VaR is a positive loss threshold; alpha is its confidence as a probability.
    double forecast = 2000000, alpha = 99 / 100.0;
    // Multipliers control realised volatility and the extra volatility in the initial stress window.
    double volMult = 1, clusterMult = 2;
    int n = 500, stressDays = 0;
    // A fixed seed makes this constructed simulation reproducible.
    long seed = 7L;
    // For zero-mean Gaussian P&L, VaR = sigma * Phi^-1(alpha).
    // Normal.inverseCdf(0.99) is about 2.32635: the standard-normal 99th percentile.
    double sigma = forecast / Normal.inverseCdf(alpha);
    Random rng = new Random(seed);
    // hit[t] is true when day t's loss exceeds the forecast; x counts these exceptions.
    boolean[] hit = new boolean[n];
    int x = 0;
    for (int t = 0; t < n; t++) {
      // Stress affects the first stressDays observations; afterwards use ordinary realised volatility.
      double scale = sigma * volMult * (t < stressDays ? clusterMult : 1.0);
      // nextGaussian draws a standard normal; positive P&L is a gain, negative P&L a loss.
      double pnl = scale * rng.nextGaussian();
      hit[t] = -pnl > forecast;
      if (hit[t]) x++;
      // Return this run's daily loss, forecast, exception flag and expected running count for the graphs.
      // Display days from 1; array indices still start at 0.
      System.out.println("DAY|" + (t + 1) + "|" + (-pnl) + "|" + forecast + "|" + (hit[t] ? 1 : 0)
          + "|" + ((t + 1) * (1 - alpha)));
    }
    // Kupiec checks the exception rate against the forecast probability 1-alpha.
    double kupiec = pof(n, x, 1 - alpha);
    // Christoffersen checks whether exception probability depends on the previous day's outcome.
    double independence = independence(hit);
    // Conditional coverage combines the frequency and independence likelihood-ratio statistics.
    double cc = kupiec + independence;
    // Output: days, exceptions, expected exceptions, then each statistic and its p-value.
    // The combined test uses two chi-square degrees of freedom: survival = exp(-cc/2).
    System.out.println("RESULT|" + n + "|" + x + "|" + (n * (1 - alpha)) + "|" + kupiec + "|" + chi1(kupiec)
        + "|" + independence + "|" + chi1(independence) + "|" + cc + "|" + Math.exp(-0.5 * cc));
  }
  // Proportion-of-failures test: compare the forecast rate p with the fitted observed rate q.
  static double pof(int n, int x, double p) {
    // Boundary cases avoid evaluating log(0) multiplied by zero.
    if (x == 0) return -2.0 * n * Math.log(1 - p);
    if (x == n) return -2.0 * n * Math.log(p);
    double q = (double) x / n;
    // -2 * (forecast log-likelihood - fitted log-likelihood).
    return -2 * ((n - x) * Math.log(1 - p) + x * Math.log(p) - (n - x) * Math.log(1 - q) - x * Math.log(q));
  }
  // A count-weighted log probability; an unobserved event contributes zero even if p=0.
  static double xl(int n, double p) {
    return n == 0 ? 0 : n * Math.log(p);
  }
  // First-order independence test over adjacent days, using 0=no exception and 1=exception.
  static double independence(boolean[] h) {
    // nij counts transitions from yesterday's state i to today's state j.
    int n00 = 0, n01 = 0, n10 = 0, n11 = 0;
    // Start at day 1 because day 0 has no preceding observation in this sample.
    for (int t = 1; t < h.length; t++) {
      // 0 -> 0: two consecutive days without an exception.
      if (!h[t - 1] && !h[t]) n00++;
      // 0 -> 1: an exception following a day without one.
      else if (!h[t - 1]) n01++;
      // 1 -> 0: a non-exception following an exception.
      else if (!h[t]) n10++;
      // 1 -> 1: consecutive exceptions (one indication of clustering).
      else n11++;
    }
    // Row totals: how many transitions began in state 0 and state 1.
    int n0 = n00 + n01;
    int n1 = n10 + n11;
    // There are h.length-1 transitions; ones counts exceptions at their destinations.
    int total = n0 + n1;
    int ones = n01 + n11;
    // No transitions, or constant destinations: both fitted likelihoods agree, so LR=0.
    // This convention does not establish independence; a sparse sample offers little evidence.
    if (total == 0 || ones == 0 || ones == total) return 0;
    // Null model: one common exception probability, independent of yesterday's state.
    double pi = (double) ones / total;
    // Alternative model: separate exception probabilities after state 0 and state 1.
    // An empty row uses 0 as a placeholder; xl makes all its contributions zero.
    double pi0 = n0 == 0 ? 0 : (double) n01 / n0;
    double pi1 = n1 == 0 ? 0 : (double) n11 / n1;
    // Null log-likelihood: all destination exceptions/non-exceptions share the pooled rate.
    double pooled = xl(ones, pi) + xl(total - ones, 1 - pi);
    // Alternative log-likelihood: each transition uses its row's fitted probability.
    double split = xl(n01, pi0) + xl(n00, 1 - pi0) + xl(n11, pi1) + xl(n10, 1 - pi1);
    // Larger LR means stronger evidence against independence; clamp round-off below zero.
    return Math.max(0, -2 * (pooled - split));
  }
  // Approximate asymptotic p-value for a chi-square statistic with one degree of freedom.
  // If Z is standard normal, Z^2 is chi-square(1): P(Z^2>s)=2*(1-Phi(sqrt(s))).
  // Small samples and rare exceptions can make this asymptotic calibration unreliable.
  // Normal.cdf(x) returns P(Z <= x); nablatensor owns the numerical approximation.
  static double chi1(double s) {
    return s <= 0 ? 1 : Math.max(0, Math.min(1, 2 * (1 - Normal.cdf(Math.sqrt(s)))));
  }
}