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

The script downloads Taiwan's official national daily rice-price series and
constructs a complete 2011-2025 monthly panel for japonica paddy, wholesale
milled rice, and retail milled rice. It estimates descriptive price gaps,
stationarity, error-correction models, asymmetric short-run responses, and
distributed-lag cumulative transmission. Reported stage gaps are not profits:
paddy and milled rice differ in physical form and the data contain no costs.
"""

from __future__ import annotations

import json
import re
from pathlib import Path
from urllib.request import Request, urlopen

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

ROOT = Path(__file__).resolve().parent
PRIVATE = ROOT / "private"
OUT = ROOT / "outputs"
FIG = OUT / "figures"
for directory in (PRIVATE, OUT, FIG):
    directory.mkdir(parents=True, exist_ok=True)

API_URL = "https://data.moa.gov.tw/Service/OpenData/Ricepriceavg.aspx?IsTransData=1&UnitId=E69"
START = pd.Timestamp("2011-01-01")
END = pd.Timestamp("2025-12-31")
COLORS = {"paddy": "#9a5b2d", "wholesale": "#d39a45", "retail": "#315c4c"}


def configure_plotting() -> None:
    candidates = [
        Path(__file__).parents[2] / "public" / "fonts" / "NotoSansTC-Regular.ttf",
        Path("C:/Windows/Fonts/msjh.ttc"),
    ]
    for candidate in candidates:
        if candidate.exists():
            font_manager.fontManager.addfont(str(candidate))
            plt.rcParams["font.family"] = font_manager.FontProperties(fname=str(candidate)).get_name()
            break
    plt.rcParams.update({"axes.unicode_minus": False, "figure.dpi": 150, "savefig.dpi": 180})
    plt.rcParams.update({"axes.grid": True, "grid.alpha": 0.18, "axes.facecolor": "#fffdf9"})


def download_daily() -> pd.DataFrame:
    target = PRIVATE / "national-daily-rice-prices-raw.json"
    request = Request(API_URL, headers={"User-Agent": "SHWRP-2026-019 research/1.0"})
    try:
        with urlopen(request, timeout=90) as response:
            raw = json.load(response)
        target.write_text(json.dumps(raw, ensure_ascii=False, indent=2), encoding="utf-8")
    except Exception:
        if not target.exists():
            raise
        raw = json.loads(target.read_text(encoding="utf-8-sig"))
    frame = pd.DataFrame(raw)

    def parse_roc(value: str) -> pd.Timestamp:
        digits = "".join(re.findall(r"\d", str(value).split(".")[0] if str(value).endswith(".0") else str(value)))
        value = digits.zfill(7)
        return pd.Timestamp(year=int(value[:3]) + 1911, month=int(value[3:5]), day=int(value[5:7]))

    frame["date"] = frame["pt_date"].map(parse_roc)
    selected = frame.loc[(frame["date"] >= START) & (frame["date"] <= END)].copy()
    rename = {
        "pt_4japt_price": "paddy_ntd_kg",
        "pt_2japt_price": "wholesale_ntd_kg",
        "pt_1japt_price": "retail_ntd_kg",
    }
    for source, destination in rename.items():
        selected[destination] = pd.to_numeric(selected[source], errors="coerce")
        if source.startswith(("pt_4", "pt_2")):
            selected[destination] /= 100.0
        selected.loc[selected[destination] <= 0, destination] = np.nan
    return selected[["date", *rename.values()]].sort_values("date")


def build_monthly(daily: pd.DataFrame) -> pd.DataFrame:
    monthly = daily.set_index("date").resample("MS").mean(numeric_only=True)
    expected = pd.date_range(START, END, freq="MS")
    monthly = monthly.reindex(expected)
    if monthly.isna().any().any():
        missing = monthly[monthly.isna().any(axis=1)].index.strftime("%Y-%m").tolist()
        raise ValueError(f"Incomplete monthly series: {missing}")
    monthly.index.name = "date"
    for stage in ("paddy", "wholesale", "retail"):
        price = f"{stage}_ntd_kg"
        monthly[f"log_{stage}"] = np.log(monthly[price])
        monthly[f"dlog_{stage}"] = monthly[f"log_{stage}"].diff()
        base = monthly.loc[monthly.index.year == 2011, price].mean()
        monthly[f"index_{stage}_2011"] = monthly[price] / base * 100
    monthly["reported_gap_wholesale_paddy"] = monthly["wholesale_ntd_kg"] - monthly["paddy_ntd_kg"]
    monthly["reported_gap_retail_wholesale"] = monthly["retail_ntd_kg"] - monthly["wholesale_ntd_kg"]
    monthly["reported_ratio_retail_paddy"] = monthly["retail_ntd_kg"] / monthly["paddy_ntd_kg"]
    monthly.to_csv(OUT / "monthly-stage-prices.csv", encoding="utf-8-sig")
    return monthly


def adf_summary(monthly: pd.DataFrame) -> pd.DataFrame:
    rows = []
    for stage in ("paddy", "wholesale", "retail"):
        for transform in ("log", "dlog"):
            series = monthly[f"{transform}_{stage}"].dropna()
            stat, pvalue, usedlag, nobs, *_ = adfuller(series, autolag="AIC")
            rows.append({"stage": stage, "transform": transform, "statistic": stat, "p_value": pvalue, "lags": usedlag, "n": nobs})
    result = pd.DataFrame(rows)
    result.to_csv(OUT / "stationarity-tests.csv", index=False, encoding="utf-8-sig")
    return result


def asymmetric_change_model(monthly: pd.DataFrame, upstream: str, downstream: str) -> tuple[dict, object]:
    long_run_x = sm.add_constant(monthly[[f"log_{upstream}"]])
    long_run = sm.OLS(monthly[f"log_{downstream}"], long_run_x).fit()
    residual_adf = adfuller(long_run.resid, autolag="AIC")
    data = pd.DataFrame(index=monthly.index)
    data["downstream_change"] = monthly[f"dlog_{downstream}"]
    change = monthly[f"dlog_{upstream}"]
    data["up_positive"] = change.clip(lower=0)
    data["up_negative"] = change.clip(upper=0)
    data["down_lag1"] = monthly[f"dlog_{downstream}"].shift(1)
    data = data.dropna()
    model = sm.OLS(data["downstream_change"], sm.add_constant(data.drop(columns="downstream_change"))).fit(cov_type="HAC", cov_kwds={"maxlags": 12})
    restriction = np.zeros((1, len(model.params)))
    restriction[0, list(model.params.index).index("up_positive")] = 1
    restriction[0, list(model.params.index).index("up_negative")] = -1
    wald = model.wald_test(restriction, scalar=True)
    record = {
        "link": f"{upstream}_to_{downstream}",
        "n": int(model.nobs),
        "level_association_slope_not_causal": float(long_run.params[f"log_{upstream}"]),
        "level_association_r2_not_causal": float(long_run.rsquared),
        "cointegration_residual_adf": float(residual_adf[0]),
        "cointegration_residual_p_approx": float(residual_adf[1]),
        "positive_short_run": float(model.params["up_positive"]),
        "positive_se": float(model.bse["up_positive"]),
        "positive_p": float(model.pvalues["up_positive"]),
        "negative_short_run": float(model.params["up_negative"]),
        "negative_se": float(model.bse["up_negative"]),
        "negative_p": float(model.pvalues["up_negative"]),
        "asymmetry_wald_p": float(wald.pvalue),
        "model_r2": float(model.rsquared),
    }
    return record, model


def distributed_lag(monthly: pd.DataFrame, upstream: str, downstream: str, max_lag: int = 6) -> tuple[pd.DataFrame, dict]:
    data = pd.DataFrame({"y": monthly[f"dlog_{downstream}"]})
    for lag in range(max_lag + 1):
        data[f"lag{lag}"] = monthly[f"dlog_{upstream}"].shift(lag)
    data["down_lag1"] = monthly[f"dlog_{downstream}"].shift(1)
    data = data.dropna()
    x_names = [f"lag{lag}" for lag in range(max_lag + 1)] + ["down_lag1"]
    model = sm.OLS(data["y"], sm.add_constant(data[x_names])).fit(cov_type="HAC", cov_kwds={"maxlags": 12})
    lag_names = [f"lag{lag}" for lag in range(max_lag + 1)]
    cov = model.cov_params().loc[lag_names, lag_names]
    coefficients = model.params[lag_names]
    rows = []
    for horizon in range(max_lag + 1):
        names = lag_names[: horizon + 1]
        estimate = coefficients.loc[names].sum()
        variance = cov.loc[names, names].to_numpy().sum()
        se = np.sqrt(max(variance, 0))
        rows.append({"link": f"{upstream}_to_{downstream}", "horizon_months": horizon, "cumulative_elasticity": estimate, "ci_low": estimate - 1.96 * se, "ci_high": estimate + 1.96 * se})
    summary = {"link": f"{upstream}_to_{downstream}", "n": int(model.nobs), "r2": float(model.rsquared), "six_month_cumulative": float(rows[-1]["cumulative_elasticity"]), "six_month_ci_low": float(rows[-1]["ci_low"]), "six_month_ci_high": float(rows[-1]["ci_high"])}
    return pd.DataFrame(rows), summary


def load_governance() -> pd.DataFrame:
    frame = pd.read_csv(ROOT / "governance-coding.csv")
    dimensions = [column for column in frame.columns if column.startswith("d_")]
    frame[dimensions] = frame[dimensions].apply(pd.to_numeric)
    frame.to_csv(OUT / "governance-coding-verified.csv", index=False, encoding="utf-8-sig")
    summary = pd.DataFrame({"dimension": [column.removeprefix("d_") for column in dimensions], "weighted_coding_score": frame[dimensions].sum(axis=0).astype(int).values, "documents_with_signal": (frame[dimensions] > 0).sum(axis=0).astype(int).values})
    summary.to_csv(OUT / "governance-dimension-summary.csv", index=False, encoding="utf-8-sig")
    return frame


def savefig(fig: plt.Figure, filename: str) -> None:
    fig.tight_layout()
    fig.savefig(FIG / filename, bbox_inches="tight", facecolor="white")
    plt.close(fig)


def figures(monthly: pd.DataFrame, ecm: list[dict], lags: pd.DataFrame, governance: pd.DataFrame) -> None:
    fig, ax = plt.subplots(figsize=(11, 5.8))
    for stage, label in (("paddy", "粳種稻穀"), ("wholesale", "粳種白米躉售"), ("retail", "粳種白米零售")):
        ax.plot(monthly.index, monthly[f"index_{stage}_2011"], label=label, color=COLORS[stage], linewidth=2)
    ax.axhline(100, color="#999", linewidth=0.8)
    ax.set(title="圖1　稻米三階段價格指數（2011平均＝100）", ylabel="價格指數")
    ax.legend(ncol=3, frameon=False)
    savefig(fig, "figure-1-stage-price-indices.png")

    fig, axes = plt.subplots(2, 1, figsize=(11, 8), sharex=True)
    for stage, label in (("paddy", "稻穀"), ("wholesale", "白米躉售"), ("retail", "白米零售")):
        axes[0].plot(monthly.index, monthly[f"{stage}_ntd_kg"], label=label, color=COLORS[stage], linewidth=1.8)
    axes[0].set(ylabel="官方報價（元／公斤）", title="圖2　官方三階段價格與表面價差")
    axes[0].legend(ncol=3, frameon=False)
    axes[1].plot(monthly.index, monthly["reported_gap_wholesale_paddy"], label="躉售－稻穀", color="#b9732e")
    axes[1].plot(monthly.index, monthly["reported_gap_retail_wholesale"], label="零售－躉售", color="#315c4c")
    axes[1].set(ylabel="表面價差（元／公斤）", xlabel="年月")
    axes[1].legend(frameon=False)
    axes[1].text(0.01, -0.32, "註：稻穀與白米物理形態不同；價差包含碾製損耗、包裝、物流、庫存、品質與服務，不是利潤。", transform=axes[1].transAxes, fontsize=9, color="#625b55")
    savefig(fig, "figure-2-reported-stage-gaps.png")

    rolling = pd.DataFrame(index=monthly.index)
    rolling["稻穀→躉售"] = monthly["dlog_paddy"].rolling(36).corr(monthly["dlog_wholesale"])
    rolling["躉售→零售"] = monthly["dlog_wholesale"].rolling(36).corr(monthly["dlog_retail"])
    rolling.to_csv(OUT / "rolling-correlations.csv", encoding="utf-8-sig")
    fig, ax = plt.subplots(figsize=(11, 5.8))
    ax.plot(rolling.index, rolling["稻穀→躉售"], label="稻穀→躉售", color="#9a5b2d", linewidth=1.8)
    ax.plot(rolling.index, rolling["躉售→零售"], label="躉售→零售", color="#315c4c", linewidth=1.8)
    ax.axhline(0, color="#777", linewidth=0.8)
    ax.set(title="圖3　三年滾動月變動相關", ylabel="Pearson相關係數", ylim=(-0.6, 1.0))
    ax.legend(frameon=False)
    savefig(fig, "figure-3-rolling-transmission.png")

    fig, axes = plt.subplots(1, 2, figsize=(11, 4.8), sharey=True)
    for ax, (link, group) in zip(axes, lags.groupby("link", sort=False)):
        ax.plot(group["horizon_months"], group["cumulative_elasticity"], marker="o", color="#315c4c")
        ax.fill_between(group["horizon_months"], group["ci_low"], group["ci_high"], color="#315c4c", alpha=0.16)
        ax.axhline(0, color="#777", linewidth=0.8)
        ax.set(title="稻穀→躉售" if link.startswith("paddy") else "躉售→零售", xlabel="累積月份")
    axes[0].set_ylabel("累積價格傳遞彈性（95%信賴區間）")
    fig.suptitle("圖4　上游價格變動的累積傳遞", y=1.02)
    savefig(fig, "figure-4-cumulative-pass-through.png")

    plot_rows = []
    for record in ecm:
        for direction in ("positive", "negative"):
            plot_rows.append({"link": "稻穀→躉售" if record["link"].startswith("paddy") else "躉售→零售", "direction": "上游上漲" if direction == "positive" else "上游下跌", "estimate": record[f"{direction}_short_run"], "se": record[f"{direction}_se"]})
    plot = pd.DataFrame(plot_rows)
    fig, ax = plt.subplots(figsize=(9, 5.5))
    x = np.arange(2)
    for offset, direction, color in ((-0.18, "上游上漲", "#9a5b2d"), (0.18, "上游下跌", "#315c4c")):
        subset = plot[plot["direction"] == direction]
        ax.bar(x + offset, subset["estimate"], width=0.34, label=direction, color=color, alpha=0.9)
        ax.errorbar(x + offset, subset["estimate"], yerr=1.96 * subset["se"], fmt="none", ecolor="#211a15", capsize=4)
    ax.axhline(0, color="#777", linewidth=0.8)
    ax.set_xticks(x, ["稻穀→躉售", "躉售→零售"])
    ax.set(ylabel="當月價格傳遞係數", title="圖5　上游上漲與下跌的短期反應")
    ax.legend(frameon=False)
    savefig(fig, "figure-5-asymmetric-responses.png")

    dimensions = [column for column in governance.columns if column.startswith("d_")]
    labels = {"d_price_rule": "價格形成", "d_quality": "品質規格", "d_purchase": "收購承諾", "d_market_access": "市場進入", "d_inventory": "庫存風險", "d_yield": "生產風險", "d_payment": "付款", "d_dispute": "爭議處理", "d_transparency": "資訊透明", "d_traceability": "可追溯性"}
    matrix = governance.set_index("short_name")[dimensions].rename(columns=labels)
    fig, ax = plt.subplots(figsize=(11, 6.4))
    from matplotlib.colors import ListedColormap
    image = ax.imshow(matrix.to_numpy(), cmap=ListedColormap(["#f4eee7", "#d39a45", "#315c4c"]), vmin=0, vmax=2, aspect="auto")
    ax.set_xticks(range(len(matrix.columns)), matrix.columns)
    ax.set_yticks(range(len(matrix.index)), matrix.index)
    for row in range(matrix.shape[0]):
        for column in range(matrix.shape[1]):
            value = int(matrix.iloc[row, column])
            ax.text(column, row, str(value), ha="center", va="center", color="white" if value == 2 else "#211a15", fontsize=9)
    colorbar = fig.colorbar(image, ax=ax, ticks=[0, 1, 2], fraction=0.03, pad=0.02)
    colorbar.set_label("0未聚焦／1提及／2明確規範或實證")
    ax.set(title="圖6　制度與文獻的關係治理資訊矩陣", xlabel="治理面向", ylabel="文件")
    ax.tick_params(axis="x", rotation=35)
    ax.tick_params(axis="y", rotation=0)
    savefig(fig, "figure-6-governance-matrix.png")


def main() -> None:
    configure_plotting()
    daily = download_daily()
    monthly = build_monthly(daily)
    tests = adf_summary(monthly)
    ecm_records = []
    for upstream, downstream in (("paddy", "wholesale"), ("wholesale", "retail")):
        record, _ = asymmetric_change_model(monthly, upstream, downstream)
        ecm_records.append(record)
    pd.DataFrame(ecm_records).to_csv(OUT / "asymmetric-change-results.csv", index=False, encoding="utf-8-sig")
    lag_frames, lag_summaries = [], []
    for upstream, downstream in (("paddy", "wholesale"), ("wholesale", "retail")):
        frame, summary = distributed_lag(monthly, upstream, downstream)
        lag_frames.append(frame)
        lag_summaries.append(summary)
    lags = pd.concat(lag_frames, ignore_index=True)
    lags.to_csv(OUT / "distributed-lag-results.csv", index=False, encoding="utf-8-sig")
    governance = load_governance()
    figures(monthly, ecm_records, lags, governance)

    summary = {
        "sample_start": str(monthly.index.min().date()),
        "sample_end": str(monthly.index.max().date()),
        "months": int(len(monthly)),
        "daily_observations": int(len(daily)),
        "mean_prices_ntd_kg": {stage: float(monthly[f"{stage}_ntd_kg"].mean()) for stage in ("paddy", "wholesale", "retail")},
        "end_prices_ntd_kg": {stage: float(monthly.iloc[-1][f"{stage}_ntd_kg"]) for stage in ("paddy", "wholesale", "retail")},
        "mean_reported_gaps_ntd_kg": {
            "wholesale_minus_paddy": float(monthly["reported_gap_wholesale_paddy"].mean()),
            "retail_minus_wholesale": float(monthly["reported_gap_retail_wholesale"].mean()),
        },
        "adf": tests.to_dict(orient="records"),
        "asymmetric_change_models": ecm_records,
        "distributed_lag": lag_summaries,
        "governance_documents": int(len(governance)),
        "warning": "Reported stage price gaps are not processor or distributor profits and do not allocate value added." 
    }
    (ROOT / "analysis-results.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()