#!/usr/bin/env python3
# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca
"""Phase 6 — comparison with external indices (Teranet–National Bank, CREA MLS HPI), harmonised to 2021 average = 100.
SHARP metropolitan aggregates are stock-weighted bottom-up combinations of municipal cells inside the CMM (Montréal) and
CMQ (Québec) metropolitan communities (SDA `comet_s`), the closest available match to the CMA used by Teranet/CREA.
Writes results/validation/external_comparison.csv and results/validation/external_series.csv."""
from __future__ import annotations

import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
import numpy as np
import polars as pl

from sharp_hpi import config as C

OUT = C.RESULTS / "validation"
EXT = C.EXTERNAL / "external_indices"


def rebase(df: pl.DataFrame, col: str = "value") -> pl.DataFrame:
    b = df.filter(pl.col("period").str.starts_with("2021")).select(pl.col(col).log().mean().alias("b"))["b"][0]
    return df.with_columns((pl.col(col).log() - b).alias("li"))


def main() -> None:
    panel = pl.read_parquet(C.PROCESSED / "index_panel.parquet")
    geo = pl.read_parquet(C.GEO_TX, columns=["mun_code", "comet_code"]).unique().filter(pl.col("comet_code").is_not_null())
    stock = pl.read_csv(OUT / "stock_weights_by_municipality_class.csv")
    series = []
    for cls in ["house", "condo"]:
        mun = panel.filter((pl.col("geography_level") == "municipality") & (pl.col("property_type") == cls) & ~pl.col("partial_month_flag")).with_columns(np.log(pl.col("index_value") / 100).alias("li"))
        mun = mun.join(geo, left_on="geography_id", right_on="mun_code", how="inner").join(stock.filter(pl.col("cls") == cls).select(["geography_id", "stock_units"]), on="geography_id", how="inner")
        for com, nm in [("663", "Montréal (CMM)"), ("235", "Québec (CMQ)")]:
            a = mun.filter(pl.col("comet_code") == com).group_by("month").agg(((pl.col("li") * pl.col("stock_units")).sum() / pl.col("stock_units").sum()).alias("li")).sort("month")
            series.append(a.select(pl.lit("SHARP-HPI").alias("source"), pl.lit(f"SHARP {cls} stock-weighted").alias("series_name"), pl.lit(nm).alias("geography"), pl.col("month").alias("period"), pl.col("li")))
        p = panel.filter((pl.col("geography_level") == "province") & (pl.col("property_type") == cls) & ~pl.col("partial_month_flag")).sort("month")
        series.append(p.select(pl.lit("SHARP-HPI").alias("source"), pl.lit(f"SHARP {cls} province").alias("series_name"), pl.lit("Quebec (province)").alias("geography"), pl.col("month").alias("period"), np.log(pl.col("index_value") / 100).alias("li")))
    ext = []
    t = pl.read_csv(EXT / "teranet_nb_hpi_monthly.csv").filter(pl.col("series_name").str.contains("not seasonally")).filter(pl.col("geography").is_in(["Montréal CMA", "Québec CMA"]))
    for g, sub in t.group_by("geography"):
        ext.append(rebase(sub.sort("period")).select(pl.lit("Teranet-NB").alias("source"), pl.lit("Teranet NSA (repeat sales)").alias("series_name"), pl.col("geography"), "period", "li"))
    c = pl.read_csv(EXT / "crea_mls_hpi_monthly.csv").filter(pl.col("series_name").str.contains("NSA") & pl.col("series_name").str.ends_with("_HPI"))
    for g, sn in [("Montréal CMA", "Single_Family"), ("Montréal CMA", "Apartment"), ("Montréal CMA", "Composite"), ("Québec CMA", "Single_Family"), ("Québec CMA", "Apartment"), ("Québec CMA", "Composite")]:
        sub = c.filter((pl.col("geography") == g) & pl.col("series_name").str.contains(sn))
        if sub.height:
            ext.append(rebase(sub.sort("period")).select(pl.lit("CREA MLS HPI").alias("source"), pl.lit(f"CREA {sn.replace('_', ' ')} NSA").alias("series_name"), pl.col("geography"), "period", "li"))
    qc = c.filter(pl.col("geography").str.contains("Quebec") & ~pl.col("geography").str.contains("CMA"))
    for sn in ["Single_Family", "Apartment", "Composite"]:
        sub = qc.filter(pl.col("series_name").str.contains(sn))
        if sub.height:
            ext.append(rebase(sub.sort("period")).select(pl.lit("CREA MLS HPI").alias("source"), pl.lit(f"CREA {sn.replace('_', ' ')} NSA").alias("series_name"), pl.lit("Quebec (province)").alias("geography"), "period", "li"))
    allser = pl.concat(series + ext); allser.write_csv(OUT / "external_series.csv")
    # pairwise comparison: SHARP vs each external in the same metro / province
    rows = []
    pairs = [("SHARP house stock-weighted", "Montréal (CMM)", "Montréal CMA"), ("SHARP condo stock-weighted", "Montréal (CMM)", "Montréal CMA"),
             ("SHARP house stock-weighted", "Québec (CMQ)", "Québec CMA"), ("SHARP condo stock-weighted", "Québec (CMQ)", "Québec CMA"),
             ("SHARP house province", "Quebec (province)", "Quebec (province)"), ("SHARP condo province", "Quebec (province)", "Quebec (province)")]
    for sname, sgeo, egeo in pairs:
        s = allser.filter((pl.col("series_name") == sname) & (pl.col("geography") == sgeo)).select("period", pl.col("li").alias("li_s"))
        for (src, ename), e in allser.filter((pl.col("geography") == egeo) & (pl.col("source") != "SHARP-HPI")).group_by(["source", "series_name"]):
            j = s.join(e.select("period", pl.col("li").alias("li_e")), on="period").sort("period")
            if j.height < 24:
                continue
            a, b = j["li_s"].to_numpy(), j["li_e"].to_numpy()
            d12a, d12b = a[12:] - a[:-12], b[12:] - b[:-12]; d1a, d1b = np.diff(a), np.diff(b)
            rows.append({"sharp_series": sname, "sharp_geo": sgeo, "external_source": src, "external_series": ename, "external_geo": egeo, "months": j.height, "last_period": j["period"][-1],
                         "rmse_level": float(np.sqrt(np.mean((a - b) ** 2))), "max_abs_gap": float(np.max(np.abs(a - b))), "total_change_sharp": float(a[-1] - a[0]), "total_change_external": float(b[-1] - b[0]),
                         "corr_12m": float(np.corrcoef(d12a, d12b)[0, 1]), "corr_1m": float(np.corrcoef(d1a, d1b)[0, 1]), "sd_1m_sharp": float(np.std(d1a)), "sd_1m_external": float(np.std(d1b))})
    out = pl.DataFrame(rows); out.write_csv(OUT / "external_comparison.csv"); print(out)


if __name__ == "__main__" and len(sys.argv) == 1:
    main()


def yearly_and_repeat_sales() -> None:
    """Year-by-year gap decomposition and our own Case–Shiller repeat sales on the CMM / island of Montréal, for comparison with Teranet."""
    from sharp_hpi import estimators as E
    from sharp_hpi.hedonic import build_features
    f = build_features(pl.read_parquet(C.ANALYTIC_TX))
    geo = pl.read_parquet(C.GEO_TX, columns=["id", "comet_code"]).rename({"id": "tx_id"})
    f = f.join(geo, on="tx_id", how="left")
    grid = E.month_grid("2021-01", "2026-10")
    rows = []
    for nm, sub in [("Montréal (CMM) houses", f.filter((pl.col("comet_code") == "663") & (pl.col("property_class") == "house"))),
                    ("Montréal (CMM) condos", f.filter((pl.col("comet_code") == "663") & (pl.col("property_class") == "condo"))),
                    ("Montréal (CMM) all three classes", f.filter(pl.col("comet_code") == "663")),
                    ("Island of Montréal (region 06) all", f.filter(pl.col("region_code") == "06")),
                    ("Québec (CMQ) all three classes", f.filter(pl.col("comet_code") == "235")),
                    ("Province houses", f.filter(pl.col("property_class") == "house")), ("Province all three classes", f)]:
        for hold in (6, 12):
            rs = E.repeat_sales(sub, grid, weighted=True, min_hold_months=hold); li = rs.log_index[:-1]; okm = ~np.isnan(li)
            med = E.median_index(sub, grid).log_index[:-1]
            yr = {}
            for y in range(2021, 2027):
                idx = [i for i, m in enumerate(grid[:-1]) if m.startswith(str(y))]
                yr[y] = (float(np.nanmean(li[idx])), float(np.nanmean(med[idx])))
            rows.append({"sample": nm, "min_hold_months": hold, "pairs": rs.meta.get("pairs", 0), "rs_total_change": float(li[okm][-1] - li[okm][0]), "median_total_change": float(med[-1] - med[0]),
                         **{f"rs_{y}_vs_2021": yr[y][0] - yr[2021][0] for y in range(2022, 2027)}, **{f"median_{y}_vs_2021": yr[y][1] - yr[2021][1] for y in range(2022, 2027)}})
    pl.DataFrame(rows).write_csv(OUT / "external_repeat_sales_check.csv"); print(pl.DataFrame(rows))
    # yearly averages of every series for the gap decomposition
    s = pl.read_csv(OUT / "external_series.csv").with_columns(pl.col("period").str.slice(0, 4).alias("year")).filter(pl.col("period") <= "2026-09")
    y = s.group_by(["source", "series_name", "geography", "year"]).agg(pl.col("li").mean()).sort(["source", "series_name", "geography", "year"])
    y.pivot(on="year", index=["source", "series_name", "geography"], values="li").write_csv(OUT / "external_yearly.csv")


if __name__ == "__main__" and len(sys.argv) > 1 and sys.argv[1] == "yearly":
    yearly_and_repeat_sales()