"""Stylized scenario simulation for SHWRP-2026-014. Outputs are conditional model results for a synthetic household population. They are not survey estimates, causal effects, or adoption forecasts. """ from __future__ import annotations import json from pathlib import Path import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap import numpy as np import pandas as pd ROOT = Path(__file__).resolve().parent OUT = ROOT / "outputs" FIG = OUT / "figures" OUT.mkdir(parents=True, exist_ok=True) FIG.mkdir(parents=True, exist_ok=True) SEED = 20260715 N = 100_000 RICE_KG_PERSON_YEAR = 42.42 HOUSEHOLD_SIZE = 2.35 REFERENCE_PRICE = 100.0 # NTD/kg; model normalization, not a market-price estimate REGIMES = ["人工按需購買", "固定商家訂閱", "彈性商家訂閱", "商家AI補貨", "開放式AI代理補貨"] COLORS = ["#8f8174", "#c07c48", "#ddb976", "#608b75", "#2f5f55"] SCENARIOS = { "2026基準": {"search": 0.15, "forecast": 0.25, "trust": 0.25, "interop": 0.15}, "2030過渡": {"search": 0.50, "forecast": 0.60, "trust": 0.55, "interop": 0.50}, "2035代理成熟": {"search": 0.80, "forecast": 0.82, "trust": 0.75, "interop": 0.78}, } def synthetic_households(n: int = N, seed: int = SEED) -> pd.DataFrame: rng = np.random.default_rng(seed) home_share = rng.triangular(0.40, 0.60, 0.80, n) annual_kg = RICE_KG_PERSON_YEAR * HOUSEHOLD_SIZE * home_share return pd.DataFrame({ "annual_kg": annual_kg, "home_share": home_share, "time_value_hour": np.clip(rng.lognormal(np.log(240), 0.35, n), 100, 650), "stockout_event_cost": rng.triangular(80, 180, 420, n), "control_value": np.clip(rng.gamma(2.0, 105.0, n), 0, 700), "data_risk_value": np.clip(rng.gamma(2.0, 100.0, n), 0, 700), "relationship_value": np.clip(rng.normal(240, 150, n), 0, 800), "demand_cv": rng.triangular(0.10, 0.25, 0.60, n), }) def generalized_costs( households: pd.DataFrame, levers: dict[str, float], reference_price: float = REFERENCE_PRICE, delivery_multiplier: float = 1.0, ) -> pd.DataFrame: q = households["annual_kg"].to_numpy() tv = households["time_value_hour"].to_numpy() stock = households["stockout_event_cost"].to_numpy() control = households["control_value"].to_numpy() data_risk = households["data_risk_value"].to_numpy() relationship = households["relationship_value"].to_numpy() cv = households["demand_cv"].to_numpy() search = levers["search"] forecast = levers["forecast"] trust = levers["trust"] interop = levers["interop"] def cost(price_index, shipments, shipping, minutes, stock_events, overstock_share, control_share, data_share, mismatch_scale, relationship_share=0.0, agent_fee=0.0): purchase = q * reference_price * price_index logistics = shipments * shipping * delivery_multiplier time_cost = tv * minutes / 60.0 stockout = stock * stock_events overstock = q * reference_price * overstock_share * 0.50 mismatch = mismatch_scale * np.square(cv) return (purchase + logistics + time_cost + stockout + overstock + control * control_share + data_risk * data_share + mismatch + agent_fee - relationship * relationship_share) manual = cost(1.00, 6, 30, 6 * 12, 0.50, 0.01, 0.00, 0.00, 50) fixed = cost(0.93, 12, 30, 60 + 12, 0.12, 0.08, 0.80, 0.10, 800, 0.55) flexible = cost(0.98, 10, 35, 50 + 10 * 2, 0.08, 0.03, 0.25, 0.15, 250, 0.55) seller_minutes = 65 + 9 * max(1.0, 5.0 * (1.0 - trust)) seller = cost( 0.99 - 0.055 * search * forecast, 9, 34, seller_minutes, 0.30 * (1.0 - 0.70 * forecast), 0.05 * (1.0 - 0.80 * forecast), max(0.12, 0.45 - 0.25 * trust), max(0.15, 0.40 - 0.20 * trust), 150 * (1.0 - 0.65 * forecast), 0.70, ) open_minutes = 50 + 8 * max(1.0, 6.0 * (1.0 - trust)) open_agent = cost( 1.00 - 0.10 * search * interop, 8, 35, open_minutes, 0.40 * (1.0 - 0.75 * forecast), 0.04 * (1.0 - 0.75 * forecast), max(0.05, 0.25 - 0.20 * trust), max(0.08, 0.50 - 0.42 * trust), 110 * (1.0 - 0.70 * forecast), 0.10, agent_fee=40 + 140 * (1.0 - interop), ) return pd.DataFrame(np.column_stack([manual, fixed, flexible, seller, open_agent]), columns=REGIMES) def chosen_shares(costs: pd.DataFrame) -> pd.Series: chosen = costs.idxmin(axis=1) return chosen.value_counts(normalize=True).reindex(REGIMES, fill_value=0.0) def scenario_analysis(households: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame]: shares = [] medians = [] for name, levers in SCENARIOS.items(): costs = generalized_costs(households, levers) share = chosen_shares(costs) shares.append(pd.DataFrame({"scenario": name, "regime": REGIMES, "share": share.values})) medians.append(pd.DataFrame({ "scenario": name, "regime": REGIMES, "median_generalized_cost_ntd_year": costs.median().values, })) return pd.concat(shares, ignore_index=True), pd.concat(medians, ignore_index=True) def phase_grid() -> tuple[np.ndarray, np.ndarray, np.ndarray]: representative = pd.DataFrame({ "annual_kg": [RICE_KG_PERSON_YEAR * HOUSEHOLD_SIZE * 0.60], "home_share": [0.60], "time_value_hour": [240.0], "stockout_event_cost": [180.0], "control_value": [210.0], "data_risk_value": [200.0], "relationship_value": [240.0], "demand_cv": [0.25], }) xs = np.linspace(0, 1, 51) ys = np.linspace(0, 1, 51) grid = np.zeros((len(ys), len(xs)), dtype=int) for iy, trust in enumerate(ys): for ix, search in enumerate(xs): levers = { "search": float(search), "trust": float(trust), "forecast": float(0.15 + 0.70 * (search + trust) / 2), "interop": float(0.10 + 0.80 * search), } costs = generalized_costs(representative, levers) grid[iy, ix] = int(np.argmin(costs.iloc[0].to_numpy())) return xs, ys, grid def industry_scores() -> pd.DataFrame: # Structured scenario assumptions, not measured industry outcomes. rows = [ ["衛生紙", .95, .85, .60, .70, .90, .85, .90, .75, .90], ["寵物食品", .90, .80, .80, .65, .75, .60, .85, .55, .80], ["白米", .88, .75, .45, .65, .72, .65, .82, .45, .78], ["咖啡豆", .78, .65, .40, .55, .65, .50, .78, .60, .45], ["生鮮蔬果", .90, .35, .75, .55, .35, .25, .65, .20, .65], ["服飾", .35, .30, .20, .35, .30, .20, .55, .70, .15], ["家電", .12, .15, .55, .25, .55, .40, .82, .60, .30], ["旅遊", .10, .10, .45, .40, .20, .15, .55, .35, .05], ] cols = ["industry", "recurrence", "predictability", "stockout_cost", "ordering_friction", "standardization", "mismatch_tolerance", "machine_readability", "reversibility", "low_experiential_choice"] df = pd.DataFrame(rows, columns=cols) weights = { "recurrence": .15, "predictability": .15, "stockout_cost": .10, "ordering_friction": .10, "standardization": .15, "mismatch_tolerance": .15, "machine_readability": .10, "reversibility": .05, "low_experiential_choice": .05, } df["automation_suitability_score"] = sum(df[k] * v for k, v in weights.items()) * 100 return df.sort_values("automation_suitability_score", ascending=False).reset_index(drop=True) def sensitivity(households: pd.DataFrame) -> pd.DataFrame: rows = [] base = SCENARIOS["2035代理成熟"] for price in [70.0, 100.0, 130.0]: for delivery in [0.65, 1.0, 1.40]: for trust in [0.55, 0.75, 0.90]: levers = {**base, "trust": trust} share = chosen_shares(generalized_costs(households, levers, price, delivery)) rows.append({ "reference_price_ntd_kg": price, "delivery_multiplier": delivery, "trust": trust, **{regime: share[regime] for regime in REGIMES}, }) return pd.DataFrame(rows) def charts(shares: pd.DataFrame, medians: pd.DataFrame, industries: pd.DataFrame, sens: pd.DataFrame) -> None: plt.rcParams.update({"font.family": ["Microsoft JhengHei", "DejaVu Sans"], "axes.unicode_minus": False}) fig, ax = plt.subplots(figsize=(12, 6.6)) stages = ["固定配送", "彈性訂閱", "資料輔助", "商家AI補貨", "開放式AI代理"] subtitles = ["固定量／固定日", "可跳過與改量", "依消耗修正", "單一商家內預測", "跨商家持續授權"] for i, (stage, subtitle) in enumerate(zip(stages, subtitles)): x = i * 2.15 ax.add_patch(plt.Rectangle((x, 0.7), 1.75, 1.2, facecolor=COLORS[i], edgecolor="#5b4637", linewidth=1.2)) ax.text(x + .875, 1.35, stage, ha="center", va="center", color="white", fontsize=13, weight="bold") ax.text(x + .875, 1.06, subtitle, ha="center", va="center", color="white", fontsize=9) if i < 4: ax.annotate("", xy=(x + 2.08, 1.3), xytext=(x + 1.78, 1.3), arrowprops=dict(arrowstyle="->", lw=1.8, color="#6f3519")) ax.text(4.95, 0.25, "制度核心由『定期扣款』轉向『可撤回、可驗證、有限範圍的持續授權』", ha="center", fontsize=12, color="#315c4c") ax.set_xlim(-.25, 10.5); ax.set_ylim(0, 2.25); ax.axis("off") ax.set_title("圖1　白米補貨制度的五階段演化假說", fontsize=15, weight="bold") fig.tight_layout(); fig.savefig(FIG / "fig1_evolution_stages.png", dpi=220); plt.close(fig) pivot = shares.pivot(index="scenario", columns="regime", values="share").reindex(SCENARIOS.keys()) fig, ax = plt.subplots(figsize=(11.5, 6.8)) bottom = np.zeros(len(pivot)) for regime, color in zip(REGIMES, COLORS): values = pivot[regime].to_numpy() * 100 ax.bar(pivot.index, values, bottom=bottom, label=regime, color=color) bottom += values ax.set_ylabel("合成家戶中廣義成本最低的比例（%）") ax.set_title("圖2　不同技術情境下的制度選擇（條件式模擬，非採用率預測）", weight="bold") ax.legend(ncol=2, frameon=False, loc="upper center", bbox_to_anchor=(.5, -0.12)) ax.set_ylim(0, 100); ax.grid(axis="y", alpha=.18) fig.tight_layout(); fig.savefig(FIG / "fig2_synthetic_choice_shares.png", dpi=220); plt.close(fig) xs, ys, grid = phase_grid() fig, ax = plt.subplots(figsize=(10.5, 7.2)) im = ax.imshow(grid, origin="lower", extent=[0, 1, 0, 1], aspect="auto", cmap=ListedColormap(COLORS), vmin=-.5, vmax=4.5) cbar = fig.colorbar(im, ax=ax, ticks=range(5)); cbar.ax.set_yticklabels(REGIMES) ax.set_xlabel("AI 搜尋／跨店比較效率"); ax.set_ylabel("授權與付款信任") ax.set_title("圖3　代表性家戶的制度優勢區域", weight="bold") fig.tight_layout(); fig.savefig(FIG / "fig3_regime_phase_map.png", dpi=220); plt.close(fig) fig, ax = plt.subplots(figsize=(10.8, 6.6)) plot = industries.sort_values("automation_suitability_score") ax.barh(plot["industry"], plot["automation_suitability_score"], color="#608b75") for i, value in enumerate(plot["automation_suitability_score"]): ax.text(value + 1, i, f"{value:.1f}", va="center", fontsize=9) ax.set_xlim(0, 100); ax.set_xlabel("補貨自動化適配分數（結構化假設）") ax.set_title("圖4　不是所有產業都同樣適合訂閱化", weight="bold") ax.grid(axis="x", alpha=.2) fig.tight_layout(); fig.savefig(FIG / "fig4_industry_suitability.png", dpi=220); plt.close(fig) mp = medians.pivot(index="regime", columns="scenario", values="median_generalized_cost_ntd_year").reindex(REGIMES) fig, ax = plt.subplots(figsize=(11.5, 6.6)) x = np.arange(len(REGIMES)); width = .24 for j, scenario in enumerate(SCENARIOS.keys()): ax.bar(x + (j - 1) * width, mp[scenario], width=width, label=scenario, color=["#c9b8a8", "#8aa796", "#315c4c"][j]) ax.set_xticks(x); ax.set_xticklabels(REGIMES, rotation=12) ax.set_ylabel("年度廣義成本中位數（模型化新臺幣）") ax.set_title("圖5　制度成本隨 AI 能力成熟而改變", weight="bold") ax.legend(frameon=False); ax.grid(axis="y", alpha=.2) fig.tight_layout(); fig.savefig(FIG / "fig5_generalized_costs.png", dpi=220); plt.close(fig) open_col = "開放式AI代理補貨" summary = sens.groupby("trust")[open_col].agg(["min", "median", "max"]).reset_index() fig, ax = plt.subplots(figsize=(10.5, 6.4)) ax.fill_between(summary["trust"], summary["min"] * 100, summary["max"] * 100, color="#a8c5b7", alpha=.55, label="價格×配送成本敏感度範圍") ax.plot(summary["trust"], summary["median"] * 100, marker="o", color="#315c4c", lw=2.5, label="中位情境") ax.set_xlabel("AI 授權與付款信任"); ax.set_ylabel("開放式 AI 代理為最低成本的合成家戶比例（%）") ax.set_title("圖6　開放式代理的擴張取決於信任，並非技術自動發生", weight="bold") ax.set_ylim(0, 100); ax.grid(alpha=.2); ax.legend(frameon=False) fig.tight_layout(); fig.savefig(FIG / "fig6_sensitivity_open_agent.png", dpi=220); plt.close(fig) def main() -> None: households = synthetic_households() shares, medians = scenario_analysis(households) industries = industry_scores() sens = sensitivity(households) shares.to_csv(OUT / "scenario_choice_shares.csv", index=False, encoding="utf-8-sig") medians.to_csv(OUT / "scenario_generalized_costs.csv", index=False, encoding="utf-8-sig") industries.to_csv(OUT / "industry_suitability_scores.csv", index=False, encoding="utf-8-sig") sens.to_csv(OUT / "sensitivity_results.csv", index=False, encoding="utf-8-sig") households.describe(percentiles=[.1, .5, .9]).to_csv(OUT / "synthetic_household_summary.csv", encoding="utf-8-sig") demand = { "official_rice_supply_kg_person_year_2024": RICE_KG_PERSON_YEAR, "official_average_household_size_may_2026": HOUSEHOLD_SIZE, "gross_household_rice_supply_kg_year": RICE_KG_PERSON_YEAR * HOUSEHOLD_SIZE, "home_share_range": [0.40, 0.80], "modeled_home_rice_kg_month_range": [ RICE_KG_PERSON_YEAR * HOUSEHOLD_SIZE * 0.40 / 12, RICE_KG_PERSON_YEAR * HOUSEHOLD_SIZE * 0.80 / 12, ], "baseline_home_rice_kg_month": RICE_KG_PERSON_YEAR * HOUSEHOLD_SIZE * 0.60 / 12, } analysis_summary = { "seed": SEED, "synthetic_households": N, "interpretation": "conditional scenario simulation; not an adoption forecast", "demand_anchor": demand, "scenario_levers": SCENARIOS, "choice_shares": { scenario: shares[shares.scenario == scenario].set_index("regime")["share"].to_dict() for scenario in SCENARIOS }, "open_agent_2035_sensitivity": { "min": float(sens["開放式AI代理補貨"].min()), "median": float(sens["開放式AI代理補貨"].median()), "max": float(sens["開放式AI代理補貨"].max()), }, } (OUT / "analysis_summary.json").write_text(json.dumps(analysis_summary, ensure_ascii=False, indent=2), encoding="utf-8") charts(shares, medians, industries, sens) print(json.dumps(analysis_summary, ensure_ascii=False, indent=2)) if __name__ == "__main__": main()