# -*- coding: utf-8 -*-
"""
공시가격 현실화율의 역진성(수직 형평성) 통계 검정
=================================================
입력 : 결과/분석_실거래가_공시가격_결합.csv

(A) 회귀분석: 현실화율 ~ ln(거래가) + 통제변수 (강건표준오차 HC1)
    - ln(거래가) 계수가 음(-)이고 유의 → 가격이 높을수록 현실화율↓ = 역진적
    - 유형 상호작용으로 아파트/연립 역진성 차이 검정
(B) IAAO 수직형평성 지표 (유형별):
    - PRD (Price-Related Differential): >1.03 역진, <0.98 누진
    - PRB (Price-Related Bias): 회귀기반, 음(-)이면 역진 (|PRB|<0.05 양호)
    - COD (Coefficient of Dispersion): 균일성
"""
from __future__ import annotations

import sys
from pathlib import Path

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

try:
    sys.stdout.reconfigure(encoding="utf-8")
except Exception:
    pass

RESULT_DIR = Path(__file__).with_name("결과")
REGION = {"11680": "강남구", "11650": "서초구", "11710": "송파구"}
TYPE_SHORT = {"아파트_매매": "아파트", "연립다세대_매매": "연립다세대"}


def load() -> pd.DataFrame:
    df = pd.read_csv(RESULT_DIR / "분석_실거래가_공시가격_결합.csv",
                     encoding="utf-8-sig", low_memory=False)
    for c in ["현실화율", "거래금액_원", "공시가격_원", "전용면적_㎡", "거래연도", "buildYear"]:
        df[c] = pd.to_numeric(df[c], errors="coerce")
    df["유형"] = df["자료구분"].map(TYPE_SHORT).fillna(df["자료구분"])
    df["지역"] = df["sggCd"].astype(str).map(REGION).fillna(df["sggCd"].astype(str))
    # 유효표본: 매칭 + 이상치 제거 + 양수 면적/가격
    m = (df["현실화율"].between(0.1, 1.5) & (df["거래금액_원"] > 0)
         & (df["전용면적_㎡"] > 0) & df["공시가격_원"].notna())
    df = df[m].copy()
    df["lnprice"] = np.log(df["거래금액_원"])
    df["lnarea"] = np.log(df["전용면적_㎡"])
    df["age"] = (df["거래연도"] - df["buildYear"]).clip(lower=0)
    df["연도"] = df["거래연도"].astype("Int64").astype(str)
    df = df.dropna(subset=["lnprice", "lnarea", "age", "현실화율"])
    print(f"분석표본 N = {len(df):,}  (아파트 {sum(df['유형']=='아파트'):,}, "
          f"연립다세대 {sum(df['유형']=='연립다세대'):,})")
    return df


def run_ols(df: pd.DataFrame, formula: str, label: str) -> sm.regression.linear_model.RegressionResultsWrapper:
    res = smf.ols(formula, data=df).fit(cov_type="HC1")
    print(f"\n{'='*70}\n[{label}]  N={int(res.nobs):,}  R²={res.rsquared:.3f}")
    print(f"  formula: {formula}")
    # 핵심 계수만 추려서 표시
    keys = [k for k in res.params.index if ("lnprice" in k or "lnarea" in k or k == "age")]
    print(f"  {'변수':<28}{'계수':>12}{'표준오차':>12}{'t':>9}{'p':>10}")
    for k in keys:
        print(f"  {k:<28}{res.params[k]:>12.5f}{res.bse[k]:>12.5f}"
              f"{res.tvalues[k]:>9.2f}{res.pvalues[k]:>10.3g}")
    return res


def iaao_metrics(g: pd.DataFrame) -> dict:
    """PRD, PRB, COD, 중앙값 현실화율."""
    ratio = g["현실화율"].values
    av = g["공시가격_원"].values      # 공시가격 (assessed)
    sp = g["거래금액_원"].values      # 실거래가 (sale price)
    med = np.median(ratio)

    cod = 100 * np.mean(np.abs(ratio - med)) / med
    prd = np.mean(ratio) / (av.sum() / sp.sum())

    # PRB (IAAO): Y=(ratio-med)/med, X=ln(0.5*(AV/med + SP))/ln(2)
    value_proxy = 0.5 * (av / med + sp)
    X = np.log(value_proxy) / np.log(2)
    Y = (ratio - med) / med
    Xc = sm.add_constant(X)
    prb_res = sm.OLS(Y, Xc).fit(cov_type="HC1")
    prb = prb_res.params[1]
    prb_p = prb_res.pvalues[1]
    prb_ci = prb_res.conf_int()[1]

    return {"N": len(g), "중앙현실화율": round(med, 4),
            "COD": round(cod, 2), "PRD": round(prd, 4),
            "PRB": round(prb, 4), "PRB_p": prb_p,
            "PRB_CI": (round(prb_ci[0], 4), round(prb_ci[1], 4))}


def interpret_prb(prb: float) -> str:
    if prb < -0.10: return "강한 역진"
    if prb < -0.05: return "역진 경향"
    if prb > 0.10:  return "강한 누진"
    if prb > 0.05:  return "누진 경향"
    return "형평(양호)"


def main() -> None:
    df = load()

    print("\n" + "#"*70 + "\n(A) 회귀분석 — ln(거래가) 계수가 음(-)이면 역진적\n" + "#"*70)
    base_ctrl = "lnarea + age + C(연도) + C(지역)"
    run_ols(df, f"현실화율 ~ lnprice + {base_ctrl} + C(유형)", "전체(pooled)")
    run_ols(df[df["유형"] == "아파트"], f"현실화율 ~ lnprice + {base_ctrl}", "아파트")
    run_ols(df[df["유형"] == "연립다세대"], f"현실화율 ~ lnprice + {base_ctrl}", "연립다세대")
    # 유형별 역진성 차이(상호작용)
    inter = run_ols(df, f"현실화율 ~ lnprice * C(유형) + {base_ctrl}",
                    "상호작용 (lnprice × 유형)")

    print("\n" + "#"*70 + "\n(B) IAAO 수직형평성 지표 (유형별·전체)\n" + "#"*70)
    rows = []
    for name, g in list(df.groupby("유형")) + [("전체", df)]:
        m = iaao_metrics(g)
        m["유형"] = name
        m["PRB판정"] = interpret_prb(m["PRB"])
        m["PRB_유의"] = "유의(p<.05)" if m["PRB_p"] < 0.05 else "비유의"
        rows.append(m)
    tbl = pd.DataFrame(rows).set_index("유형")[
        ["N", "중앙현실화율", "COD", "PRD", "PRB", "PRB_p", "PRB판정", "PRB_유의"]]
    print("\n" + tbl.to_string())
    print("\n※ 기준: PRD>1.03 역진 / PRB<-0.05 역진(음수·유의) / COD 낮을수록 균일")

    # 연도×유형 PRB 추이
    print("\n[연도별 PRB (역진성 추이)]")
    yr_rows = []
    for (typ, yr), g in df.groupby(["유형", "연도"]):
        if len(g) < 30:
            continue
        m = iaao_metrics(g)
        yr_rows.append({"유형": typ, "연도": yr, "PRB": m["PRB"],
                        "PRD": m["PRD"], "N": m["N"]})
    yr = pd.DataFrame(yr_rows)
    piv = yr.pivot(index="연도", columns="유형", values="PRB")
    print(piv.round(4).to_string())

    # 저장
    out = RESULT_DIR / "분석_역진성_회귀.xlsx"
    with pd.ExcelWriter(out) as xw:
        tbl.to_excel(xw, sheet_name="IAAO지표")
        piv.to_excel(xw, sheet_name="연도별PRB")
        # 회귀 계수표(상호작용 모델)
        coef = pd.DataFrame({"계수": inter.params, "표준오차": inter.bse,
                             "t": inter.tvalues, "p": inter.pvalues}).round(5)
        coef.to_excel(xw, sheet_name="회귀_상호작용")
    print(f"\n저장: {out}")


if __name__ == "__main__":
    main()