"""Reproducible analysis for SHWRP-2026-013.

The script uses only public aggregate data.  It separates three questions:
1. monthly price co-movement with retail rice;
2. association with the CPI for prepared Chinese rice foods;
3. seasonally adjusted price changes in the month containing Ghost Festival.

The 2026 hog-market section is an event description, not a causal estimate.
"""

from __future__ import annotations

import json
import math
import xml.etree.ElementTree as ET
from pathlib import Path

import matplotlib.pyplot as plt
from matplotlib import font_manager
import numpy as np
import pandas as pd
import statsmodels.api as sm


ROOT = Path(__file__).resolve().parent
RAW = ROOT / "raw"
OUT = ROOT / "outputs"
FIG = OUT / "figures"
OUT.mkdir(exist_ok=True)
FIG.mkdir(exist_ok=True)

font_path = Path(r"C:\Windows\Fonts\NotoSansTC-VF.ttf")
if font_path.exists():
    font_manager.fontManager.addfont(font_path)
    chart_font = font_manager.FontProperties(fname=font_path).get_name()
else:
    chart_font = "sans-serif"

plt.rcParams.update(
    {
        "font.family": chart_font,
        "axes.unicode_minus": False,
        "figure.dpi": 160,
        "savefig.dpi": 220,
        "axes.spines.top": False,
        "axes.spines.right": False,
    }
)

ITEMS = {
    "總指數": "總指數(指數基期：民國110年=100)",
    "米": "1米(指數基期：民國110年=100)",
    "豬肉": "11豬肉(指數基期：民國110年=100)",
    "牛肉": "13牛肉、牛內臟(指數基期：民國110年=100)",
    "雞肉": "15雞肉(指數基期：民國110年=100)",
    "雞蛋": "20雞蛋(指數基期：民國110年=100)",
    "水產品": "5.水產品(指數基期：民國110年=100)",
    "蔬菜": "7.蔬菜(指數基期：民國110年=100)",
    "水果": "9.水果(指數基期：民國110年=100)",
    "食用油": "12.食用油(指數基期：民國110年=100)",
    "中式米食": "156中式米食(指數基期：民國110年=100)",
    "中式麵食": "157中式麵食(指數基期：民國110年=100)",
}

CANDIDATES = ["豬肉", "牛肉", "雞肉", "雞蛋", "水產品", "蔬菜", "水果", "食用油"]

# Gregorian dates of lunar 7/15.  The date determines the monthly indicator;
# it does not imply that all observances or purchases occur on one day.
GHOST_FESTIVAL = {
    2013: "2013-08-21",
    2014: "2014-08-10",
    2015: "2015-08-28",
    2016: "2016-08-17",
    2017: "2017-09-05",
    2018: "2018-08-25",
    2019: "2019-08-15",
    2020: "2020-09-02",
    2021: "2021-08-22",
    2022: "2022-08-12",
    2023: "2023-08-30",
    2024: "2024-08-18",
    2025: "2025-09-06",
    2026: "2026-08-27",
}


def parse_cpi() -> pd.DataFrame:
    rows: list[tuple[str, str, float]] = []
    wanted = set(ITEMS.values())
    for _, elem in ET.iterparse(RAW / "cpi-item-groups.xml", events=("end",)):
        if elem.tag != "Obs":
            continue
        row = {c.tag: c.text for c in elem}
        if row.get("TYPE") == "原始值" and row.get("Item") in wanted:
            rows.append((row["TIME_PERIOD"], row["Item"], float(row["Item_VALUE"])))
        elem.clear()
    long = pd.DataFrame(rows, columns=["period", "item", "index"])
    reverse = {v: k for k, v in ITEMS.items()}
    long["item"] = long["item"].map(reverse)
    long["date"] = pd.to_datetime(long["period"].str.replace("M", "-", regex=False) + "-01")
    wide = long.pivot(index="date", columns="item", values="index").sort_index()
    missing = sorted(set(ITEMS) - set(wide.columns))
    if missing:
        raise ValueError(f"CPI items missing: {missing}")
    return wide


def hac_regression(y: pd.Series, x: pd.Series, controls: pd.DataFrame | None = None) -> dict:
    frame = pd.concat([y.rename("y"), x.rename("x"), controls], axis=1).dropna()
    rhs = ["x"] + ([] if controls is None else list(controls.columns))
    X = sm.add_constant(frame[rhs])
    fit = sm.OLS(frame["y"], X).fit(cov_type="HAC", cov_kwds={"maxlags": 12})
    return {
        "n": int(fit.nobs),
        "beta": float(fit.params["x"]),
        "se_hac": float(fit.bse["x"]),
        "p_hac": float(fit.pvalues["x"]),
        "ci_low": float(fit.conf_int().loc["x", 0]),
        "ci_high": float(fit.conf_int().loc["x", 1]),
        "r2": float(fit.rsquared),
    }


def cpi_analysis(cpi: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame]:
    log_index = np.log(cpi)
    mom = log_index.diff() * 100
    yoy = log_index.diff(12) * 100

    rows = []
    for product in CANDIDATES:
        pair_m = mom[["米", product]].dropna()
        pair_y = yoy[["米", product]].dropna()
        rows.append(
            {
                "product": product,
                "n_mom": len(pair_m),
                "pearson_mom": pair_m["米"].corr(pair_m[product]),
                "spearman_mom": pair_m["米"].corr(pair_m[product], method="spearman"),
                "n_yoy": len(pair_y),
                "pearson_yoy": pair_y["米"].corr(pair_y[product]),
                "spearman_yoy": pair_y["米"].corr(pair_y[product], method="spearman"),
            }
        )
    corr = pd.DataFrame(rows).sort_values("pearson_yoy", ascending=False)

    # One candidate at a time: downstream prepared-rice-food inflation against
    # candidate inflation, controlling for all-items CPI and calendar-month FE.
    month_fe = pd.get_dummies(yoy.index.month, prefix="month", drop_first=True, dtype=float)
    month_fe.index = yoy.index
    prepared_rows = []
    for product in ["米"] + CANDIDATES:
        controls = pd.concat([yoy["總指數"].rename("all_cpi"), month_fe], axis=1)
        result = hac_regression(yoy["中式米食"], yoy[product], controls)
        prepared_rows.append({"product": product, **result})
    prepared = pd.DataFrame(prepared_rows).sort_values("beta", ascending=False)

    # Festival-month coefficient with month fixed effects and a linear trend.
    mom_model = mom.loc[:"2025-12-01"].copy()
    festival_months = {pd.Timestamp(v).to_period("M") for k, v in GHOST_FESTIVAL.items() if k <= 2025}
    festival = pd.Series(
        [int(d.to_period("M") in festival_months) for d in mom_model.index],
        index=mom_model.index,
        name="festival",
        dtype=float,
    )
    fe = pd.get_dummies(mom_model.index.month, prefix="month", drop_first=True, dtype=float)
    fe.index = mom_model.index
    trend = pd.Series(np.arange(len(mom_model), dtype=float), index=mom_model.index, name="trend")
    festival_rows = []
    for product in ["米", "中式米食"] + CANDIDATES:
        result = hac_regression(mom_model[product], festival, pd.concat([trend, fe], axis=1))
        festival_rows.append({"product": product, **result})
    festival_results = pd.DataFrame(festival_rows).sort_values("beta", ascending=False)

    cpi.to_csv(OUT / "cpi_selected_indices.csv", encoding="utf-8-sig")
    mom.to_csv(OUT / "cpi_log_changes_mom.csv", encoding="utf-8-sig")
    yoy.to_csv(OUT / "cpi_log_changes_yoy.csv", encoding="utf-8-sig")
    corr.to_csv(OUT / "rice_price_correlations.csv", index=False, encoding="utf-8-sig")
    prepared.to_csv(OUT / "prepared_rice_food_models.csv", index=False, encoding="utf-8-sig")
    festival_results.to_csv(OUT / "festival_month_models.csv", index=False, encoding="utf-8-sig")

    return corr, prepared, festival_results


def roc_to_timestamp(value: int | str) -> pd.Timestamp:
    s = str(value).split(".")[0].zfill(7)
    return pd.Timestamp(year=int(s[:3]) + 1911, month=int(s[3:5]), day=int(s[5:7]))


def hog_analysis() -> tuple[pd.DataFrame, dict]:
    raw = pd.read_csv(RAW / "pig-transactions.csv")
    raw["date"] = raw["交易日期"].map(roc_to_timestamp)
    heads = pd.to_numeric(raw["成交總數(不含冷凍廠)-頭數"], errors="coerce").fillna(0)
    price = pd.to_numeric(raw["成交總數(不含冷凍廠)-平均價格"], errors="coerce")
    raw["heads"] = heads
    raw["price"] = price
    raw["value_proxy"] = raw["heads"] * raw["price"]
    daily = raw.groupby("date", as_index=False).agg(heads=("heads", "sum"), value_proxy=("value_proxy", "sum"))
    daily["weighted_price"] = daily["value_proxy"] / daily["heads"].replace(0, np.nan)
    daily = daily.drop(columns="value_proxy").sort_values("date")
    daily["price_7d"] = daily["weighted_price"].rolling(7, min_periods=3).mean()
    daily["heads_7d"] = daily["heads"].rolling(7, min_periods=3).mean()

    def window(start: str, end: str) -> dict:
        x = daily[(daily.date >= start) & (daily.date <= end)]
        return {
            "start": start,
            "end": end,
            "trading_days": int(len(x)),
            "weighted_price": float(np.average(x.weighted_price, weights=x.heads)),
            "mean_daily_heads": float(x.heads.mean()),
            "total_heads": int(x.heads.sum()),
        }

    current = window("2026-07-01", "2026-07-14")
    prior_year = window("2025-07-01", "2025-07-14")
    baseline = window("2026-01-01", "2026-05-31")
    latest = daily.loc[daily.date.idxmax()]
    summary = {
        "data_start": daily.date.min().date().isoformat(),
        "data_end": daily.date.max().date().isoformat(),
        "current_july": current,
        "same_period_2025": prior_year,
        "jan_may_2026": baseline,
        "current_vs_2025_price_pct": (current["weighted_price"] / prior_year["weighted_price"] - 1) * 100,
        "current_vs_2025_heads_pct": (current["mean_daily_heads"] / prior_year["mean_daily_heads"] - 1) * 100,
        "latest_date": latest.date.date().isoformat(),
        "latest_weighted_price": float(latest.weighted_price),
        "latest_heads": int(latest.heads),
    }
    daily.to_csv(OUT / "hog_market_daily_aggregate.csv", index=False, encoding="utf-8-sig")
    (OUT / "hog_event_summary.json").write_text(json.dumps(summary, ensure_ascii=False, indent=2), encoding="utf-8")
    return daily, summary


def make_figures(cpi: pd.DataFrame, corr: pd.DataFrame, prepared: pd.DataFrame, festival: pd.DataFrame, hog: pd.DataFrame) -> None:
    colors = {"米": "#7a3f1f", "豬肉": "#c96c45", "雞肉": "#d6a848", "雞蛋": "#7b8f4e", "水產品": "#417a89", "蔬菜": "#4d7851"}

    fig, ax = plt.subplots(figsize=(10.5, 5.4))
    for col in ["米", "豬肉", "雞肉", "雞蛋", "水產品", "蔬菜"]:
        ax.plot(cpi.index, cpi[col], label=col, lw=1.7, color=colors[col])
    ax.axhline(100, color="#999", lw=0.8, ls="--")
    ax.set(title="主要農產品消費者物價指數（2021=100）", ylabel="指數")
    ax.legend(ncol=3, frameon=False)
    ax.grid(axis="y", alpha=.2)
    fig.tight_layout()
    fig.savefig(FIG / "fig1_cpi_indices.png", bbox_inches="tight")
    plt.close(fig)

    plot = corr.sort_values("pearson_yoy")
    fig, ax = plt.subplots(figsize=(8.4, 5.1))
    bars = ax.barh(plot["product"], plot.pearson_yoy, color=["#7a3f1f" if v >= 0 else "#8aa0a6" for v in plot.pearson_yoy])
    ax.axvline(0, color="#333", lw=.8)
    ax.bar_label(bars, labels=[f"{v:.2f}" for v in plot.pearson_yoy], padding=3, fontsize=9)
    ax.set(xlabel="與米價年增率之 Pearson 相關係數", title="價格共變排名：年增率，不使用價格水準")
    ax.grid(axis="x", alpha=.2)
    fig.tight_layout()
    fig.savefig(FIG / "fig2_rice_correlations.png", bbox_inches="tight")
    plt.close(fig)

    plot = prepared.sort_values("beta")
    fig, ax = plt.subplots(figsize=(8.4, 5.2))
    err = np.vstack([plot.beta - plot.ci_low, plot.ci_high - plot.beta])
    ax.errorbar(plot.beta, plot["product"], xerr=err, fmt="o", color="#7a3f1f", ecolor="#bda48f", capsize=3)
    ax.axvline(0, color="#333", lw=.8)
    ax.set(xlabel="中式米食年增率對該品項年增率之係數（95% HAC CI）", title="餐食端關聯：逐項模型，控制總 CPI 與月份")
    ax.grid(axis="x", alpha=.2)
    fig.tight_layout()
    fig.savefig(FIG / "fig3_prepared_rice_models.png", bbox_inches="tight")
    plt.close(fig)

    plot = festival.sort_values("beta")
    fig, ax = plt.subplots(figsize=(8.5, 5.5))
    err = np.vstack([plot.beta - plot.ci_low, plot.ci_high - plot.beta])
    ax.errorbar(plot.beta, plot["product"], xerr=err, fmt="o", color="#7a3f1f", ecolor="#bda48f", capsize=3)
    ax.axvline(0, color="#333", lw=.8)
    ax.set(xlabel="中元節所在月的額外月增率（百分點，95% HAC CI）", title="節慶月份訊號：控制月份固定效果與趨勢")
    ax.grid(axis="x", alpha=.2)
    fig.tight_layout()
    fig.savefig(FIG / "fig4_festival_effects.png", bbox_inches="tight")
    plt.close(fig)

    x = hog[hog.date >= "2025-01-01"].copy()
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10.4, 6.7), sharex=True)
    ax1.plot(x.date, x.weighted_price, color="#d4a48a", alpha=.4, lw=.7)
    ax1.plot(x.date, x.price_7d, color="#7a3f1f", lw=2, label="7交易日移動平均")
    ax1.set(ylabel="加權平均價格（元／公斤）", title="臺灣毛豬交易市場：價格與成交頭數")
    ax1.legend(frameon=False)
    ax1.grid(axis="y", alpha=.2)
    ax2.plot(x.date, x.heads, color="#9ab2a0", alpha=.35, lw=.7)
    ax2.plot(x.date, x.heads_7d, color="#315d48", lw=2)
    ax2.set(ylabel="成交頭數", xlabel="交易日期")
    ax2.grid(axis="y", alpha=.2)
    event_date = pd.Timestamp("2026-07-15")
    for ax in [ax1, ax2]:
        ax.axvline(event_date, color="#a33", ls="--", lw=1)
        ax.text(event_date, ax.get_ylim()[1], " 7/15臨時休市（樣本外）", va="top", fontsize=9, color="#a33")
    fig.tight_layout()
    fig.savefig(FIG / "fig5_hog_event.png", bbox_inches="tight")
    plt.close(fig)


def main() -> None:
    cpi = parse_cpi()
    corr, prepared, festival = cpi_analysis(cpi)
    hog, hog_summary = hog_analysis()
    make_figures(cpi, corr, prepared, festival, hog)
    summary = {
        "cpi_start": cpi.index.min().date().isoformat(),
        "cpi_end": cpi.index.max().date().isoformat(),
        "cpi_months": int(len(cpi)),
        "top_rice_yoy_correlation": corr.iloc[0].to_dict(),
        "top_prepared_rice_beta": prepared.iloc[0].to_dict(),
        "largest_festival_beta": festival.iloc[0].to_dict(),
        "hog": hog_summary,
    }
    (OUT / "analysis_summary.json").write_text(json.dumps(summary, ensure_ascii=False, indent=2), encoding="utf-8")
    print(json.dumps(summary, ensure_ascii=False, indent=2))


if __name__ == "__main__":
    main()