// Builds the saved analyses for the FactHub dataset "homelessness-rate-vs-rent".
// Inputs: ../../data/coc_homelessness_rent_panel.csv, ../../data/peers/peers.csv
// Outputs (in ../out/): king_rent_residuals.csv, analyses.json
// Usage: node build_rate_rent.js
// King County is left out of each fit (leave-one-out), so its 95% interval is a
// true out-of-sample prediction interval for a CoC with King's rent.
const fs = require('fs');
const path = require('path');

const DATA = path.join(__dirname, '..', '..', 'data');
const OUT = path.join(__dirname, '..', 'out');
fs.mkdirSync(OUT, { recursive: true });

function parseCsv(text) {
  const rows = [];
  for (const line of text.split(/\r?\n/)) {
    if (!line) continue;
    const cells = [];
    let cur = '', q = false;
    for (let i = 0; i < line.length; i++) {
      const c = line[i];
      if (q) { if (c === '"') { if (line[i + 1] === '"') { cur += '"'; i++; } else q = false; } else cur += c; }
      else if (c === '"') q = true;
      else if (c === ',') { cells.push(cur); cur = ''; }
      else cur += c;
    }
    cells.push(cur);
    rows.push(cells);
  }
  const [h, ...body] = rows;
  return body.map(r => Object.fromEntries(h.map((k, i) => [k, r[i]])));
}

// Two-sided 95% t quantile (Cornish-Fisher expansion; accurate to ~0.001 for df >= 20).
function t975(df) {
  const z = 1.959964;
  return z + (z ** 3 + z) / (4 * df) + (5 * z ** 5 + 16 * z ** 3 + 3 * z) / (96 * df * df);
}

function ols(xs, ys) {
  const n = xs.length;
  const xbar = xs.reduce((a, b) => a + b, 0) / n;
  const ybar = ys.reduce((a, b) => a + b, 0) / n;
  let sxx = 0, sxy = 0, syy = 0;
  for (let i = 0; i < n; i++) { sxx += (xs[i] - xbar) ** 2; sxy += (xs[i] - xbar) * (ys[i] - ybar); syy += (ys[i] - ybar) ** 2; }
  const b = sxy / sxx, a = ybar - b * xbar;
  let sse = 0;
  for (let i = 0; i < n; i++) sse += (ys[i] - a - b * xs[i]) ** 2;
  const s = Math.sqrt(sse / (n - 2));
  return { n, a, b, s, xbar, sxx, r2: 1 - sse / syy, seB: s / Math.sqrt(sxx) };
}

const r2 = v => Math.round(v * 100) / 100;
const panel = parseCsv(fs.readFileSync(path.join(DATA, 'coc_homelessness_rent_panel.csv'), 'utf8'));
const YEARS = [2015, 2016, 2017, 2018, 2019, 2022, 2024];

// sample: 'all' | 'majorCity' | 'exNYC'; measure: 'rate_per_10k' | 'unsheltered_rate'
function kingFit(year, sample, measure, kingAdjustPeople = 0) {
  const rows = panel.filter(r => +r.year === year && r.sheltered_only_count === '0' && r.rate_per_10k !== '' && r.median_rent !== '');
  const king = rows.find(r => r.coc_code === 'WA-500');
  let fitRows = rows.filter(r => r.coc_code !== 'WA-500');
  if (sample === 'majorCity') fitRows = fitRows.filter(r => r.hud_category_2025 === 'Major City CoC');
  if (sample === 'exNYC') fitRows = fitRows.filter(r => r.coc_code !== 'NY-600');
  const m = ols(fitRows.map(r => +r.median_rent / 100), fitRows.map(r => +r[measure]));
  const x0 = +king.median_rent / 100;
  const pred = m.a + m.b * x0;
  const actual = +king[measure] - kingAdjustPeople / +king.population * 10000;
  const sePred = m.s * Math.sqrt(1 + 1 / m.n + (x0 - m.xbar) ** 2 / m.sxx);
  const t = t975(m.n - 2);
  const resid = actual - pred;
  return {
    year, label: `${year} King County ${measure === 'rate_per_10k' ? 'total' : 'unsheltered'} rate minus rent-predicted rate`,
    estimate: r2(resid), ciLow: r2(resid - t * sePred), ciHigh: r2(resid + t * sePred), ciLevel: 0.95,
    standardError: r2(sePred), unit: 'homeless people per 10,000 residents',
    method: 'Leave-one-out OLS on median rent; 95% prediction interval',
    n: m.n, actual: r2(actual), predicted: r2(pred), slopePer100: r2(m.b), slopeSE: r2(m.seB), r2: Math.round(m.r2 * 1000) / 1000,
    kingMedianRent: +king.median_rent,
  };
}

const analyses = {};
const series = (sample, measure) => YEARS.map(y => kingFit(y, sample, measure));
analyses.main = series('all', 'rate_per_10k');
analyses.majorCity = series('majorCity', 'rate_per_10k');
analyses.exNYC = series('exNYC', 'rate_per_10k');
analyses.unsheltered = series('all', 'unsheltered_rate');
analyses.unshelteredMajorCity = series('majorCity', 'unsheltered_rate');
analyses.methodAdj2024 = [
  { ...kingFit(2024, 'all', 'rate_per_10k', 1100), label: '2024 total, central method adjustment (1,100 people removed)' },
  { ...kingFit(2024, 'all', 'rate_per_10k', 3600), label: '2024 total, high method adjustment (3,600 people removed)' },
  { ...kingFit(2024, 'majorCity', 'rate_per_10k', 1100), label: '2024 total, Major City sample, central adjustment' },
  { ...kingFit(2024, 'majorCity', 'rate_per_10k', 3600), label: '2024 total, Major City sample, high adjustment' },
];

// Rent growth 2019-24 vs change in homelessness, across the other 66 CoCs (from peers.csv).
const peers = parseCsv(fs.readFileSync(path.join(DATA, 'peers', 'peers.csv'), 'utf8'));
const others = peers.filter(r => r.coc_code !== 'WA-500' && r.rent_growth_2019_24_pct !== '' && r.chg_rate_pct !== '');
function corr(xs, ys) {
  const n = xs.length, mx = xs.reduce((a, b) => a + b) / n, my = ys.reduce((a, b) => a + b) / n;
  let sxy = 0, sxx = 0, syy = 0;
  for (let i = 0; i < n; i++) { sxy += (xs[i] - mx) * (ys[i] - my); sxx += (xs[i] - mx) ** 2; syy += (ys[i] - my) ** 2; }
  const r = sxy / Math.sqrt(sxx * syy);
  const z = Math.atanh(r), se = 1 / Math.sqrt(n - 3);
  return { r, lo: Math.tanh(z - 1.96 * se), hi: Math.tanh(z + 1.96 * se), n };
}
const rg = others.map(r => +r.rent_growth_2019_24_pct);
const mk = (field, label) => {
  const c = corr(rg, others.map(r => +r[field]));
  return { label, estimate: Math.round(c.r * 1000) / 1000, ciLow: Math.round(c.lo * 1000) / 1000, ciHigh: Math.round(c.hi * 1000) / 1000, ciLevel: 0.95, unit: 'Pearson correlation coefficient', method: 'Pearson r across large CoCs excluding King County; 95% CI by Fisher z-transform', n: c.n };
};
analyses.rentGrowthCorr = [
  mk('chg_rate_pct', 'Rent growth 2019-24 vs % change in total homelessness rate'),
  mk('chg_unsheltered_pct', 'Rent growth 2019-24 vs % change in unsheltered rate'),
  mk('chg_rate', 'Rent growth 2019-24 vs change in total rate per 10,000'),
];
const all = peers.filter(r => r.rent_growth_2019_24_pct !== '');
const rank = (field, desc) => [...all].sort((a, b) => desc ? +b[field] - +a[field] : +a[field] - +b[field]).findIndex(r => r.coc_code === 'WA-500') + 1;
analyses.kingRanks = {
  rows: [
    { measure: 'Rent growth 2019-24 (%)', kingValue: 21.5, rankOf: all.length, rank: rank('rent_growth_2019_24_pct', true), rankDirection: '1 = fastest rent growth' },
    { measure: 'Change in unsheltered rate per 10,000, 2019 to latest full count', kingValue: 18.71, rankOf: all.length, rank: rank('chg_unsheltered_rate', true), rankDirection: '1 = largest increase' },
    { measure: 'Change in total rate per 10,000, 2019 to latest full count', kingValue: 22.37, rankOf: all.length, rank: rank('chg_rate', true), rankDirection: '1 = largest increase' },
  ],
};

// Derived table: one row per year x sample x measure for King County.
const lines = ['year,sample,measure,n,king_median_rent,actual_rate,predicted_rate,residual,pi95_low,pi95_high,slope_per_100,slope_se,r2'];
for (const [key, sample, measure] of [['main', 'all_full_count', 'total'], ['majorCity', 'major_city', 'total'], ['exNYC', 'excl_nyc', 'total'], ['unsheltered', 'all_full_count', 'unsheltered'], ['unshelteredMajorCity', 'major_city', 'unsheltered']]) {
  for (const r of analyses[key]) lines.push([r.year, sample, measure, r.n, r.kingMedianRent, r.actual, r.predicted, r.estimate, r.ciLow, r.ciHigh, r.slopePer100, r.slopeSE, r.r2].join(','));
}
fs.writeFileSync(path.join(OUT, 'king_rent_residuals.csv'), lines.join('\n') + '\n');
fs.writeFileSync(path.join(OUT, 'analyses.json'), JSON.stringify(analyses, null, 1));
for (const k of ['main', 'majorCity', 'unsheltered', 'unshelteredMajorCity']) {
  console.log(k);
  for (const r of analyses[k]) console.log(`  ${r.year} n=${r.n} actual=${r.actual} pred=${r.predicted} resid=${r.estimate} [${r.ciLow}, ${r.ciHigh}] R2=${r.r2}`);
}
console.log('methodAdj2024'); for (const r of analyses.methodAdj2024) console.log(`  ${r.label}: ${r.estimate} [${r.ciLow}, ${r.ciHigh}]`);
console.log('rentGrowthCorr'); for (const r of analyses.rentGrowthCorr) console.log(`  ${r.label}: r=${r.estimate} [${r.ciLow}, ${r.ciHigh}] n=${r.n}`);
console.log('kingRanks', JSON.stringify(analyses.kingRanks.rows));