// Difference-in-differences of Washington retail gasoline prices vs. comparison
// regions, using EIA weekly retail price series. Mirrors the Greene Economics
// approach in Ecology Pub. 26-14-043 (baseline 2015-2022, CPI-adjusted to Jan 2024 $)
// and extends it through 2025 and 2026 YTD.
//
// Usage: node did.js <eia_raw_dir> <out_dir>
const fs = require('fs');
const path = require('path');

const [rawDir, outDir] = process.argv.slice(2);
const MONTHS = { Jan: 1, Feb: 2, Mar: 3, Apr: 4, May: 5, Jun: 6, Jul: 7, Aug: 8, Sep: 9, Oct: 10, Nov: 11, Dec: 12 };

function parseEia(file) {
  const html = fs.readFileSync(file, 'utf8').replace(/\s+/g, ' ');
  const series = new Map();
  const rowRe = /<td class='B6'>(?:&nbsp;)*(\d{4})-(\w{3})<\/td>(.*?)<\/tr>/g;
  let m;
  while ((m = rowRe.exec(html))) {
    const year = +m[1];
    const rowMonth = MONTHS[m[2]];
    const cellRe = /<td class='B5'>(\d{2})\/(\d{2})&nbsp;<\/td> <td class='B3'>([\d.]+)/g;
    let c;
    while ((c = cellRe.exec(m[3]))) {
      const mm = +c[1], dd = +c[2];
      const y = mm < rowMonth ? year + 1 : year;
      series.set(`${y}-${String(mm).padStart(2, '0')}-${String(dd).padStart(2, '0')}`, +c[3]);
    }
  }
  return series;
}

// CPI-U (FRED CPIAUCSL), monthly
const cpi = new Map();
for (const line of fs.readFileSync(path.join(rawDir, 'CPIAUCSL.csv'), 'utf8').trim().split(/\r?\n/).slice(1)) {
  const [d, v] = line.split(',');
  cpi.set(d.slice(0, 7), v === '' ? null : +v);
}
// BLS did not publish October 2025 CPI (federal shutdown); interpolate linearly.
for (const [k, v] of cpi) {
  if (v != null) continue;
  const keys = [...cpi.keys()].sort(), i = keys.indexOf(k);
  cpi.set(k, (cpi.get(keys[i - 1]) + cpi.get(keys[i + 1])) / 2);
  console.log(`Note: CPI for ${k} missing; interpolated as ${cpi.get(k).toFixed(3)}`);
}
const CPI_BASE = cpi.get('2024-01');
const latestCpiMonth = [...cpi.keys()].sort().pop();
function realFactor(date) {
  const k = cpi.has(date.slice(0, 7)) ? date.slice(0, 7) : latestCpiMonth; // carry latest CPI forward
  return CPI_BASE / cpi.get(k);
}

// Washington state gasoline excise tax ($/gal). Normalize to the 2016-2025 rate of $0.494
// so that state gas-tax changes are not counted as climate-policy effects.
function waTax(date) {
  if (date < '2015-08-01') return 0.375;
  if (date < '2016-07-01') return 0.445;
  if (date < '2025-07-01') return 0.494;
  if (date < '2026-07-01') return 0.554;
  return 0.565;
}

const files = fs.readdirSync(rawDir).filter((f) => f.endsWith('.html'));
const S = {};
for (const f of files) S[f.replace('.html', '')] = parseEia(path.join(rawDir, f));

const NO_CARBON_STATES = ['SCO', 'STX', 'SFL', 'SMN', 'SOH', 'SNY', 'SMA'];
const CONTROLS = {
  'U.S. average': (d) => S.NUS.get(d),
  'Rocky Mountain (PADD 4, incl. Idaho)': (d) => S.R40.get(d),
  'Colorado': (d) => S.SCO.get(d),
  '7 EIA states w/o fuel carbon pricing': (d) => {
    const v = NO_CARBON_STATES.map((s) => S[s].get(d));
    return v.some((x) => x == null) ? undefined : v.reduce((a, b) => a + b, 0) / v.length;
  },
  'California (has cap-and-trade)': (d) => S.SCA.get(d),
};

const dates = [...S.SWA.keys()].filter((d) => d >= '2015-01-01').sort();
const lastDate = dates[dates.length - 1];

// OLS of y on [1, dummies...] with Newey-West HAC standard errors.
function olsHac(y, X, lag) {
  const n = y.length, k = X[0].length;
  const XtX = Array.from({ length: k }, () => Array(k).fill(0));
  const Xty = Array(k).fill(0);
  for (let i = 0; i < n; i++) for (let a = 0; a < k; a++) { Xty[a] += X[i][a] * y[i]; for (let b = 0; b < k; b++) XtX[a][b] += X[i][a] * X[i][b]; }
  const inv = invert(XtX);
  const beta = inv.map((r) => r.reduce((s, v, j) => s + v * Xty[j], 0));
  const e = y.map((yi, i) => yi - X[i].reduce((s, v, j) => s + v * beta[j], 0));
  const meat = Array.from({ length: k }, () => Array(k).fill(0));
  for (let l = 0; l <= lag; l++) {
    const w = l === 0 ? 1 : 1 - l / (lag + 1);
    for (let t = l; t < n; t++) for (let a = 0; a < k; a++) for (let b = 0; b < k; b++) {
      const v = w * e[t] * e[t - l] * X[t][a] * X[t - l][b];
      meat[a][b] += v; if (l > 0) meat[b][a] += v;
    }
  }
  const V = mul(mul(inv, meat), inv);
  return { beta, se: beta.map((_, i) => Math.sqrt(V[i][i])) };
}
function mul(A, B) { return A.map((r) => B[0].map((_, j) => r.reduce((s, v, i) => s + v * B[i][j], 0))); }
function invert(M) {
  const n = M.length, A = M.map((r, i) => [...r, ...Array.from({ length: n }, (_, j) => (i === j ? 1 : 0))]);
  for (let c = 0; c < n; c++) {
    let p = c; for (let r = c + 1; r < n; r++) if (Math.abs(A[r][c]) > Math.abs(A[p][c])) p = r;
    [A[c], A[p]] = [A[p], A[c]];
    const d = A[c][c]; for (let j = 0; j < 2 * n; j++) A[c][j] /= d;
    for (let r = 0; r < n; r++) if (r !== c) { const f = A[r][c]; for (let j = 0; j < 2 * n; j++) A[r][j] -= f * A[c][j]; }
  }
  return A.map((r) => r.slice(n));
}

// Builds weekly gap rows: WA (tax-normalized, real) minus control (real).
function gaps(ctrlFn, { taxAdjust = true, real = true } = {}) {
  const rows = [];
  for (const d of dates) {
    const wa = S.SWA.get(d), c = ctrlFn(d);
    if (wa == null || c == null) continue;
    const f = real ? realFactor(d) : 1;
    const waAdj = taxAdjust ? wa - (waTax(d) - 0.494) : wa;
    rows.push({ d, gap: waAdj * f - c * f });
  }
  return rows;
}

// Effect for each period = mean gap in period - mean gap in baseline (HAC SE).
function did(rows, periodOf, baseline = (d) => d < '2023-01-01', lag = 12) {
  const rs = rows.filter((r) => baseline(r.d) || periodOf(r.d));
  const periods = [...new Set(rs.map((r) => periodOf(r.d)).filter(Boolean))].sort();
  const X = rs.map((r) => [1, ...periods.map((p) => (periodOf(r.d) === p ? 1 : 0))]);
  const { beta, se } = olsHac(rs.map((r) => r.gap), X, lag);
  return { baselineGap: beta[0], effects: periods.map((p, i) => ({ period: p, effect: beta[i + 1], se: se[i + 1], weeks: rs.filter((r) => periodOf(r.d) === p).length })) };
}

const byYear = (d) => (d >= '2023-01-01' ? d.slice(0, 4) : null);
const byQuarter = (d) => (d >= '2023-01-01' ? `${d.slice(0, 4)}Q${Math.ceil(+d.slice(5, 7) / 3)}` : null);
const fmt = (x) => (x * 100).toFixed(1).padStart(6);

const results = { lastDate, cpiBase: '2024-01', latestCpiMonth, annual: {}, quarterly: {}, checks: {} };
console.log(`EIA weekly data through ${lastDate}; CPI through ${latestCpiMonth}; cents/gal, Jan-2024 dollars, WA gas tax normalized\n`);
console.log('Annual Washington-specific price effect vs. 2015-2022 baseline (cents/gal, ±HAC SE)');
console.log('Control'.padEnd(40) + ['2023', '2024', '2025', '2026 YTD'].map((s) => s.padStart(16)).join(''));
for (const [name, fn] of Object.entries(CONTROLS)) {
  const r = did(gaps(fn), byYear);
  results.annual[name] = r;
  console.log(name.padEnd(40) + r.effects.map((e) => `${fmt(e.effect)} ±${(e.se * 100).toFixed(1).padStart(4)}   `).join(''));
}

console.log('\nQuarterly effect (cents/gal)');
const qControls = ['U.S. average', 'Rocky Mountain (PADD 4, incl. Idaho)', '7 EIA states w/o fuel carbon pricing'];
const q = {};
for (const name of qControls) q[name] = did(gaps(CONTROLS[name]), byQuarter);
results.quarterly = q;
console.log('Quarter'.padEnd(10) + qControls.map((n) => n.slice(0, 22).padStart(24)).join(''));
for (let i = 0; i < q[qControls[0]].effects.length; i++) {
  console.log(q[qControls[0]].effects[i].period.padEnd(10) + qControls.map((n) => fmt(q[n].effects[i].effect).padStart(24)).join(''));
}

// Checks
// 1. Placebo: pretend a policy began in 2019 with a 2015-2018 baseline.
const placebo = {};
for (const name of qControls) {
  placebo[name] = did(gaps(CONTROLS[name]).filter((r) => r.d < '2023-01-01'),
    (d) => (d >= '2019-01-01' ? d.slice(0, 4) : null), (d) => d < '2019-01-01');
}
results.checks.placebo2019 = placebo;
console.log('\nPlacebo (fake policy start 2019, baseline 2015-2018), cents/gal');
for (const [n, r] of Object.entries(placebo)) console.log(n.padEnd(40) + r.effects.map((e) => `${e.period}:${fmt(e.effect)}`).join('  '));

// 2. Sensitivity: baseline 2019-2022 only.
const shortBase = {};
for (const name of qControls) shortBase[name] = did(gaps(CONTROLS[name]), byYear, (d) => d >= '2019-01-01' && d < '2023-01-01');
results.checks.baseline2019to2022 = shortBase;
console.log('\nSensitivity: baseline 2019-2022 only');
for (const [n, r] of Object.entries(shortBase)) console.log(n.padEnd(40) + r.effects.map((e) => `${e.period}:${fmt(e.effect)}`).join('  '));

// 3. Without tax normalization and in nominal dollars (comparability).
const raw = {};
for (const name of qControls) raw[name] = did(gaps(CONTROLS[name], { taxAdjust: false, real: false }), byYear);
results.checks.nominalNoTaxAdjust = raw;
console.log('\nNominal $, no gas-tax normalization');
for (const [n, r] of Object.entries(raw)) console.log(n.padEnd(40) + r.effects.map((e) => `${e.period}:${fmt(e.effect)}`).join('  '));

// 4. Nominal $, gas-tax normalized (used for dollar cost estimates).
const nominal = {};
for (const name of qControls) nominal[name] = did(gaps(CONTROLS[name], { taxAdjust: true, real: false }), byYear);
results.checks.nominalTaxNormalized = nominal;
console.log('\nNominal $, gas-tax normalized');
for (const [n, r] of Object.entries(nominal)) console.log(n.padEnd(40) + r.effects.map((e) => `${e.period}:${fmt(e.effect)}`).join('  '));

fs.mkdirSync(outDir, { recursive: true });
fs.writeFileSync(path.join(outDir, 'did_results.json'), JSON.stringify(results, null, 2));
// Weekly table for the dataset
const weekly = dates.map((d) => ({ weekEnding: d, WA: S.SWA.get(d), US: S.NUS.get(d), PADD4: S.R40.get(d), CA: S.SCA.get(d),
  CO: S.SCO.get(d), TX: S.STX.get(d), FL: S.SFL.get(d), MN: S.SMN.get(d), OH: S.SOH.get(d), NY: S.SNY.get(d), MA: S.SMA.get(d),
  waStateTax: waTax(d), cpiFactorJan2024: +realFactor(d).toFixed(6) }));
fs.writeFileSync(path.join(outDir, 'weekly_prices.json'), JSON.stringify(weekly));