"""Reproducible analysis for SHWRP-2026-020. This study describes market-capacity dilution in Taiwan's real-estate brokerage sector. Buy/sale transfer registrations are a market-volume proxy, not broker-mediated closings. Ratios that divide transfers by businesses or personnel therefore measure available market throughput, not firm or worker productivity. Policy dates are descriptive event markers and are not treated as exogenous causal instruments. """ from __future__ import annotations import json from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd import statsmodels.api as sm from matplotlib import font_manager from scipy import stats ROOT = Path(__file__).resolve().parent OUT = ROOT / "outputs" FIG = OUT / "figures" OUT.mkdir(exist_ok=True) FIG.mkdir(exist_ok=True) MARKET_FILE = ROOT / "brokerage-market-2006-2025.csv" PIPELINE_FILE = ROOT / "housing-pipeline-2016-2025.csv" COLORS = { "transactions": "#7b3f21", "businesses": "#d08b3f", "personnel": "#315c4c", "permit": "#a05a2c", "start": "#cf8a3a", "use": "#55715f", "initial": "#335f77", } def configure_plotting() -> None: candidates = [ ROOT.parents[1] / "public" / "fonts" / "NotoSansTC-Regular.ttf", Path("C:/Windows/Fonts/msjh.ttc"), Path("C:/Windows/Fonts/mingliu.ttc"), ] for candidate in candidates: if candidate.exists(): font_manager.fontManager.addfont(str(candidate)) family = font_manager.FontProperties(fname=str(candidate)).get_name() plt.rcParams["font.family"] = family break plt.rcParams.update( { "axes.unicode_minus": False, "figure.dpi": 150, "savefig.dpi": 180, "axes.facecolor": "#fffdf8", "figure.facecolor": "#fffdf8", "axes.grid": True, "grid.alpha": 0.18, "axes.spines.top": False, "axes.spines.right": False, } ) def load_market() -> pd.DataFrame: frame = pd.read_csv(MARKET_FILE) frame["transfers_per_business"] = frame["buy_sale_transfers"] / frame["operating_brokerage_businesses"] frame["transfers_per_person"] = frame["buy_sale_transfers"] / frame["employed_brokerage_personnel"] for column in ["buy_sale_transfers", "operating_brokerage_businesses", "employed_brokerage_personnel"]: frame[f"{column}_index_2006"] = frame[column] / frame.loc[0, column] * 100 frame[f"{column}_growth_pct"] = frame[column].pct_change() * 100 return frame def log_trend(frame: pd.DataFrame, column: str) -> dict[str, float]: x = sm.add_constant(frame["year"] - frame["year"].min()) fit = sm.OLS(np.log(frame[column]), x).fit(cov_type="HC3") slope = float(fit.params["year"]) return { "annual_log_slope": slope, "annual_percent_trend": float((np.exp(slope) - 1) * 100), "p_value_hc3": float(fit.pvalues["year"]), "r_squared": float(fit.rsquared), } def best_piecewise_break(frame: pd.DataFrame, column: str, min_segment: int = 5) -> dict[str, float | int]: years = frame["year"].to_numpy() y = np.log(frame[column].to_numpy()) candidates: list[dict[str, float | int]] = [] for split in range(min_segment, len(frame) - min_segment + 1): sse = 0.0 slopes: list[float] = [] for indices in (slice(0, split), slice(split, len(frame))): xx = sm.add_constant(years[indices] - years[indices][0]) fit = sm.OLS(y[indices], xx).fit() sse += float(np.sum(fit.resid**2)) slopes.append(float(fit.params[1])) candidates.append( { "break_year": int(years[split]), "sse": sse, "pre_annual_pct": (np.exp(slopes[0]) - 1) * 100, "post_annual_pct": (np.exp(slopes[1]) - 1) * 100, } ) return min(candidates, key=lambda item: float(item["sse"])) def chow_test(frame: pd.DataFrame, column: str, break_year: int) -> dict[str, float | int]: before = frame.loc[frame["year"] < break_year] after = frame.loc[frame["year"] >= break_year] if min(len(before), len(after)) < 3: return {"break_year": break_year, "f_statistic": float("nan"), "p_value": float("nan")} def sse(part: pd.DataFrame) -> float: x = sm.add_constant(part["year"] - frame["year"].min()) return float(np.sum(sm.OLS(np.log(part[column]), x).fit().resid**2)) pooled = sse(frame) separate = sse(before) + sse(after) k = 2 denominator_df = len(frame) - 2 * k f_value = ((pooled - separate) / k) / (separate / denominator_df) return { "break_year": break_year, "f_statistic": float(f_value), "p_value": float(stats.f.sf(f_value, k, denominator_df)), } def pipeline_lags(pipeline: pd.DataFrame) -> pd.DataFrame: rows = [] changes = pipeline.set_index("year").pct_change().dropna() * 100 for upstream in ["h2_building_permit_households", "h2_construction_start_households", "h2_use_permit_households"]: for lag in range(0, 4): left = changes[upstream].iloc[: len(changes) - lag] right = changes["initial_building_registrations"].iloc[lag:] rows.append( { "upstream_indicator": upstream, "lag_years_to_initial_registration": lag, "pearson_correlation_of_annual_changes": float(np.corrcoef(left, right)[0, 1]), "n_pairs": len(left), } ) return pd.DataFrame(rows) def scenario_paths(market: pd.DataFrame) -> tuple[pd.DataFrame, dict[str, float]]: recent_growth = market.loc[market["year"].between(2017, 2025), "buy_sale_transfers_growth_pct"].dropna() contraction = float(recent_growth.quantile(0.25)) / 100 normalization = float(recent_growth.median()) / 100 assumptions = { "contraction_annual_rate": contraction, "stabilization_annual_rate": 0.0, "normalization_annual_rate": normalization, } base = float(market.loc[market["year"] == 2025, "buy_sale_transfers"].iloc[0]) rows = [] for year in range(2025, 2031): horizon = year - 2025 rows.append( { "year": year, "structural_contraction": round(base * ((1 + contraction) ** horizon)), "stabilization": round(base), "transaction_normalization": round(base * ((1 + normalization) ** horizon)), } ) return pd.DataFrame(rows), assumptions def make_figures(market: pd.DataFrame, pipeline: pd.DataFrame, lags: pd.DataFrame, scenarios: pd.DataFrame) -> None: # Figure 1: indexed divergence fig, ax = plt.subplots(figsize=(9.2, 5.2)) ax.plot(market["year"], market["buy_sale_transfers_index_2006"], marker="o", lw=2.4, color=COLORS["transactions"], label="買賣移轉棟數") ax.plot(market["year"], market["operating_brokerage_businesses_index_2006"], marker="o", lw=2.2, color=COLORS["businesses"], label="經紀業營業家數") ax.plot(market["year"], market["employed_brokerage_personnel_index_2006"], marker="o", lw=2.2, color=COLORS["personnel"], label="受僱經紀人員") ax.axhline(100, color="#777", lw=0.8) ax.set(title="交易市場與仲介供給的長期分歧", xlabel="年", ylabel="指數（2006=100）") ax.legend(frameon=False, ncol=3, loc="upper left") fig.tight_layout() fig.savefig(FIG / "figure-1-indexed-divergence.png", bbox_inches="tight") plt.close(fig) # Figure 2: opportunity-density proxies fig, ax = plt.subplots(figsize=(9.2, 5.2)) ax.plot(market["year"], market["transfers_per_business"], marker="o", lw=2.4, color=COLORS["businesses"], label="每家營業業者市場移轉量") ax.set_ylabel("棟／營業業者", color=COLORS["businesses"]) ax2 = ax.twinx() ax2.plot(market["year"], market["transfers_per_person"], marker="s", lw=2.2, color=COLORS["personnel"], label="每名登記人員市場移轉量") ax2.set_ylabel("棟／登記從業人員", color=COLORS["personnel"]) ax.set(title="案件機會密度代理指標（不是實際成交生產力）", xlabel="年") lines = ax.get_lines() + ax2.get_lines() ax.legend(lines, [line.get_label() for line in lines], frameon=False, loc="upper right") fig.tight_layout() fig.savefig(FIG / "figure-2-market-opportunity-density.png", bbox_inches="tight") plt.close(fig) # Figure 3: annual changes with policy markers fig, ax = plt.subplots(figsize=(9.2, 5.2)) growth = market.set_index("year")["buy_sale_transfers_growth_pct"] colors = ["#b34f3d" if value < 0 else "#4f7a5d" for value in growth.dropna()] ax.bar(growth.dropna().index, growth.dropna().values, color=colors, width=0.75) for year, label in [(2016, "房地合一1.0"), (2021, "房地合一2.0／實登2.0"), (2023, "平均地權修法／新青安"), (2024, "第七度信用管制")]: ax.axvline(year, color="#4b3628", lw=0.8, ls="--") ax.text(year + 0.08, ax.get_ylim()[1] * 0.88, label, rotation=90, va="top", fontsize=8, color="#4b3628") ax.axhline(0, color="#444", lw=0.8) ax.set(title="買賣移轉量年增減與政策事件", xlabel="年", ylabel="年增減（%）") fig.tight_layout() fig.savefig(FIG / "figure-3-policy-events-and-growth.png", bbox_inches="tight") plt.close(fig) # Figure 4: indexed supply pipeline fig, ax = plt.subplots(figsize=(9.2, 5.2)) mapping = { "h2_building_permit_households": ("住宅類建造執照", COLORS["permit"]), "h2_construction_start_households": ("住宅類開工", COLORS["start"]), "h2_use_permit_households": ("住宅類使用執照", COLORS["use"]), "initial_building_registrations": ("建物第一次登記", COLORS["initial"]), "buy_sale_transfers": ("建物買賣移轉", COLORS["transactions"]), } for column, (label, color) in mapping.items(): ax.plot(pipeline["year"], pipeline[column] / pipeline[column].iloc[0] * 100, marker="o", lw=2, label=label, color=color) ax.axhline(100, color="#777", lw=0.8) ax.set(title="住宅供給管線與交易量（各序列2016=100）", xlabel="年", ylabel="指數") ax.legend(frameon=False, ncol=2) fig.tight_layout() fig.savefig(FIG / "figure-4-housing-pipeline.png", bbox_inches="tight") plt.close(fig) # Figure 5: lag correlations fig, ax = plt.subplots(figsize=(9.2, 5.2)) lag_mapping = { "h2_building_permit_households": ("建照→第一次登記", COLORS["permit"]), "h2_construction_start_households": ("開工→第一次登記", COLORS["start"]), "h2_use_permit_households": ("使照→第一次登記", COLORS["use"]), } for column, (label, color) in lag_mapping.items(): part = lags.loc[lags["upstream_indicator"] == column] ax.plot(part["lag_years_to_initial_registration"], part["pearson_correlation_of_annual_changes"], marker="o", lw=2.2, label=label, color=color) ax.axhline(0, color="#555", lw=0.8) ax.set(title="供給管線與第一次登記年增率的探索性領先落後相關", xlabel="上游指標年增率領先年數", ylabel="Pearson相關係數", xticks=range(4), ylim=(-1, 1)) ax.legend(frameon=False) fig.tight_layout() fig.savefig(FIG / "figure-5-pipeline-lag-correlations.png", bbox_inches="tight") plt.close(fig) # Figure 6: scenarios fig, ax = plt.subplots(figsize=(9.2, 5.2)) actual = market.loc[market["year"] >= 2016, ["year", "buy_sale_transfers"]] ax.plot(actual["year"], actual["buy_sale_transfers"], marker="o", lw=2.4, color=COLORS["transactions"], label="歷史值") ax.plot(scenarios["year"], scenarios["structural_contraction"], marker="o", lw=2, color="#8f5260", label="結構收縮情境") ax.plot(scenarios["year"], scenarios["stabilization"], marker="o", lw=2, color="#82764f", label="低量穩定情境") ax.plot(scenarios["year"], scenarios["transaction_normalization"], marker="o", lw=2, color="#3c6c59", label="交易正常化情境") ax.axvline(2025, color="#555", lw=0.8, ls="--") ax.text(2025.1, ax.get_ylim()[1] * 0.96, "情境起點", va="top", fontsize=9) ax.set(title="2026–2030年買賣移轉量情境展望（非單點預測）", xlabel="年", ylabel="建物買賣移轉棟數") ax.legend(frameon=False) fig.tight_layout() fig.savefig(FIG / "figure-6-scenario-outlook.png", bbox_inches="tight") plt.close(fig) def main() -> None: configure_plotting() market = load_market() pipeline = pd.read_csv(PIPELINE_FILE) lags = pipeline_lags(pipeline) scenarios, assumptions = scenario_paths(market) trends = { column: log_trend(market, column) for column in ["buy_sale_transfers", "operating_brokerage_businesses", "employed_brokerage_personnel", "transfers_per_business", "transfers_per_person"] } breaks = { column: best_piecewise_break(market, column) for column in ["buy_sale_transfers", "operating_brokerage_businesses", "employed_brokerage_personnel"] } chow = { column: [chow_test(market, column, year) for year in (2016, 2021, 2024)] for column in ["buy_sale_transfers", "operating_brokerage_businesses", "employed_brokerage_personnel"] } endpoints = {} for column in ["buy_sale_transfers", "operating_brokerage_businesses", "employed_brokerage_personnel", "transfers_per_business", "transfers_per_person"]: start = float(market.loc[market["year"] == 2006, column].iloc[0]) end = float(market.loc[market["year"] == 2025, column].iloc[0]) endpoints[column] = {"2006": start, "2025": end, "percent_change": (end / start - 1) * 100} counterfactual = { "businesses_at_2019_opportunity_density": float(market.loc[market["year"] == 2025, "buy_sale_transfers"].iloc[0] / market.loc[market["year"] == 2019, "transfers_per_business"].iloc[0]), "personnel_at_2019_opportunity_density": float(market.loc[market["year"] == 2025, "buy_sale_transfers"].iloc[0] / market.loc[market["year"] == 2019, "transfers_per_person"].iloc[0]), } correlation = { "level_transfers_businesses": float(market["buy_sale_transfers"].corr(market["operating_brokerage_businesses"])), "level_transfers_personnel": float(market["buy_sale_transfers"].corr(market["employed_brokerage_personnel"])), "growth_transfers_businesses": float(market["buy_sale_transfers_growth_pct"].corr(market["operating_brokerage_businesses_growth_pct"])), "growth_transfers_personnel": float(market["buy_sale_transfers_growth_pct"].corr(market["employed_brokerage_personnel_growth_pct"])), } first_registration = pipeline.set_index("year")["initial_building_registrations"] transaction = pipeline.set_index("year")["buy_sale_transfers"] results = { "report_number": "SHWRP-2026-020", "analysis_date": "2026-07-21", "observations": {"market_years": len(market), "pipeline_years": len(pipeline)}, "endpoint_comparisons": endpoints, "robust_log_trends": trends, "exploratory_best_breaks": breaks, "predetermined_chow_tests": chow, "correlations": correlation, "2024_to_2025": { "transfers_pct": float(market.loc[market["year"] == 2025, "buy_sale_transfers_growth_pct"].iloc[0]), "businesses_pct": float(market.loc[market["year"] == 2025, "operating_brokerage_businesses_growth_pct"].iloc[0]), "personnel_pct": float(market.loc[market["year"] == 2025, "employed_brokerage_personnel_growth_pct"].iloc[0]), "initial_registration_pct": float((first_registration.loc[2025] / first_registration.loc[2024] - 1) * 100), }, "pipeline": { "initial_registration_2016": int(first_registration.loc[2016]), "initial_registration_2025": int(first_registration.loc[2025]), "percent_change": float((first_registration.loc[2025] / first_registration.loc[2016] - 1) * 100), "initial_registration_to_transactions_2016": float(first_registration.loc[2016] / transaction.loc[2016]), "initial_registration_to_transactions_2025": float(first_registration.loc[2025] / transaction.loc[2025]), }, "scenario_assumptions": assumptions, "scenario_2030": scenarios.loc[scenarios["year"] == 2030].iloc[0].to_dict(), "counterfactual_not_recommendation": counterfactual, "interpretation_guardrail": "Transfer-per-business and transfer-per-person ratios are market-capacity proxies, not actual broker closings or productivity. Policy-event timing is descriptive, not causal identification.", } market.to_csv(OUT / "market-series-with-derived-indicators.csv", index=False) pipeline.to_csv(OUT / "housing-pipeline-series.csv", index=False) lags.to_csv(OUT / "pipeline-lag-correlations.csv", index=False) scenarios.to_csv(OUT / "scenario-paths-2025-2030.csv", index=False) pd.DataFrame( [ {"series": column, **values} for column, values in trends.items() ] ).to_csv(OUT / "robust-trend-models.csv", index=False) pd.DataFrame( [ {"series": column, **values} for column, values in breaks.items() ] ).to_csv(OUT / "exploratory-change-points.csv", index=False) (ROOT / "analysis-results.json").write_text(json.dumps(results, ensure_ascii=False, indent=2), encoding="utf-8") make_figures(market, pipeline, lags, scenarios) if __name__ == "__main__": main()