← back to Bertha

src/lib/model.js

140 lines

const { query } = require('./db');

// Calibrated ensemble probability estimation
// P(event) = fraction of "forecast scenarios" exceeding threshold
function ensembleExceedance(forecastFeatures, threshold, operator = '>=') {
  const values = forecastFeatures?.values || [];
  if (!values.length) return null;

  let count = 0;
  for (const v of values) {
    if (operator === '>=' && v >= threshold) count++;
    else if (operator === '<=' && v <= threshold) count++;
    else if (operator === '>' && v > threshold) count++;
    else if (operator === '<' && v < threshold) count++;
  }

  return count / values.length;
}

// Isotonic-style calibration (simple bin-based)
function calibrate(rawProb, calibrationBins) {
  if (!calibrationBins || !calibrationBins.length) return rawProb;

  // Find the closest bin
  let closestBin = calibrationBins[0];
  let minDist = Math.abs(rawProb - closestBin.center);

  for (const bin of calibrationBins) {
    const dist = Math.abs(rawProb - bin.center);
    if (dist < minDist) {
      minDist = dist;
      closestBin = bin;
    }
  }

  return closestBin.calibrated;
}

// Simple climatological prior (historical base rate)
function climatologicalPrior(variable, threshold, month, location) {
  // Simplified priors based on US weather patterns
  const priors = {
    temperature: {
      // P(temp >= threshold) by month range (rough US averages)
      summer: (t) => t > 100 ? 0.05 : t > 90 ? 0.35 : t > 80 ? 0.65 : 0.85,
      winter: (t) => t > 60 ? 0.15 : t > 40 ? 0.50 : t > 20 ? 0.75 : 0.90,
      spring: (t) => t > 80 ? 0.20 : t > 60 ? 0.55 : t > 40 ? 0.80 : 0.95,
      fall: (t) => t > 80 ? 0.25 : t > 60 ? 0.50 : t > 40 ? 0.75 : 0.90,
    },
    precipitation: {
      // P(rain >= threshold inches)
      any: (t) => t > 2 ? 0.10 : t > 1 ? 0.20 : t > 0.5 ? 0.35 : 0.50,
    },
    snow: {
      winter: (t) => t > 12 ? 0.05 : t > 6 ? 0.15 : t > 3 ? 0.25 : 0.40,
      any: () => 0.10,
    }
  };

  const season = month >= 6 && month <= 8 ? 'summer' :
                 month >= 12 || month <= 2 ? 'winter' :
                 month >= 3 && month <= 5 ? 'spring' : 'fall';

  const varPriors = priors[variable];
  if (!varPriors) return 0.50; // uninformative
  const fn = varPriors[season] || varPriors.any;
  if (!fn) return 0.50;
  return fn(threshold);
}

// Bayesian update: combine prior with forecast evidence
function bayesianUpdate(prior, forecastProb, forecastWeight = 0.7) {
  // Weighted combination (simplified)
  return prior * (1 - forecastWeight) + forecastProb * forecastWeight;
}

// Compute Brier score for a set of predictions
function brierScore(predictions) {
  if (!predictions.length) return null;
  const sum = predictions.reduce((acc, p) => {
    const outcome = p.resolved_yes ? 1 : 0;
    return acc + Math.pow(p.p_yes - outcome, 2);
  }, 0);
  return sum / predictions.length;
}

// Full prediction pipeline
async function generatePrediction({ conditionId, forecastFeatures, eventSpec, month }) {
  const { variable, threshold, operator } = eventSpec;

  // Step 1: Ensemble exceedance probability
  const rawProb = ensembleExceedance(forecastFeatures, threshold, operator);
  if (rawProb === null) {
    return { p_yes: null, confidence: 0, method: 'no_data' };
  }

  // Step 2: Climatological prior
  const prior = climatologicalPrior(variable, threshold, month);

  // Step 3: Bayesian update
  const combined = bayesianUpdate(prior, rawProb, 0.7);

  // Step 4: Clamp to [0.01, 0.99] to avoid extreme confidence
  const pYes = Math.max(0.01, Math.min(0.99, combined));

  // Step 5: Confidence based on data quality
  const dataPoints = forecastFeatures?.values?.length || 0;
  const confidence = Math.min(0.95, dataPoints / 20); // more data = more confidence

  return {
    p_yes: Math.round(pYes * 1000) / 1000,
    confidence: Math.round(confidence * 100) / 100,
    method: 'calibrated_ensemble_v1',
    components: {
      raw_exceedance: rawProb,
      climatological_prior: prior,
      bayesian_combined: combined,
      data_points: dataPoints,
    }
  };
}

// Store prediction in DB
async function storePrediction({ conditionId, asofTs, modelVersionId, pYes, confidence, edge, marketMid, recommendation, reasoning }) {
  const res = await query(`
    INSERT INTO predictions (condition_id, asof_ts, model_version_id, p_yes, confidence, edge, market_mid, recommendation, reasoning)
    VALUES ($1, $2, $3, $4, $5, $6, $7, $8, $9)
    ON CONFLICT (condition_id, asof_ts, model_version_id) DO UPDATE SET
      p_yes = EXCLUDED.p_yes, confidence = EXCLUDED.confidence, edge = EXCLUDED.edge,
      market_mid = EXCLUDED.market_mid, recommendation = EXCLUDED.recommendation, reasoning = EXCLUDED.reasoning
    RETURNING *
  `, [conditionId, asofTs || new Date(), modelVersionId || 1, pYes, confidence, edge, marketMid, recommendation, reasoning]);
  return res.rows[0];
}

module.exports = {
  ensembleExceedance, calibrate, climatologicalPrior, bayesianUpdate,
  brierScore, generatePrediction, storePrediction
};