"""Reproducible analysis for SHWRP-2026-021. The study separates long-run employment-stock change from contemporary vacancy, recruitment-duration, and turnover signals. It is descriptive screening rather than a causal estimate or a forecast of firm-level labor shortages. """ from __future__ import annotations import json from pathlib import Path import matplotlib matplotlib.use("Agg") import matplotlib.font_manager as fm import matplotlib.pyplot as plt import numpy as np import pandas as pd ROOT = Path(__file__).resolve().parent OUT = ROOT / "outputs" FIG = OUT / "figures" OUT.mkdir(exist_ok=True) FIG.mkdir(exist_ok=True) FONT_PATH = ROOT.parents[1] / "public" / "fonts" / "NotoSansTC-Regular.ttf" if FONT_PATH.exists(): fm.fontManager.addfont(str(FONT_PATH)) plt.rcParams["font.family"] = "Noto Sans TC" plt.rcParams.update({"axes.unicode_minus": False, "figure.dpi": 170, "savefig.dpi": 190}) LABELS = { "agriculture": "農林漁牧", "mining": "礦業土石", "manufacturing": "製造", "electricity_gas": "電力燃氣", "water_remediation": "用水污染整治", "construction": "營建工程", "wholesale_retail": "批發零售", "transport_storage": "運輸倉儲", "accommodation_food": "住宿餐飲", "information_communication": "出版影音資通訊", "finance_insurance": "金融保險", "real_estate": "不動產", "professional_scientific": "專業科學技術", "support_services": "支援服務", "public_administration": "公共行政國防", "education": "教育", "health_social": "醫療社會工作", "arts_recreation": "藝術休閒", "other_services": "其他服務", } emp = pd.read_csv(ROOT / "employment-by-industry-2014-2025.csv") vac = pd.read_csv(ROOT / "vacancy-march-2025.csv") agr = pd.read_csv(ROOT / "agriculture-succession-demand.csv") records = [] for col, label in LABELS.items(): start = float(emp.loc[emp.year == 2014, col].iloc[0]) end = float(emp.loc[emp.year == 2025, col].iloc[0]) pct = (end / start - 1) * 100 cagr = (end / start) ** (1 / 11) - 1 slope = float(np.polyfit(emp.year, emp[col], 1)[0]) records.append({ "industry_code": col, "industry": label, "employment_2014_thousand": start, "employment_2025_thousand": end, "change_thousand": end - start, "percent_change": pct, "cagr_percent": cagr * 100, "ols_slope_thousand_per_year": slope, }) trend = pd.DataFrame(records).merge( vac[["employment_column", "vacancies", "vacancy_rate", "entry_rate", "exit_rate", "recruitment_months"]], how="left", left_on="industry_code", right_on="employment_column" ).drop(columns=["employment_column"]) trend["turnover_gap_pp"] = trend["exit_rate"] - trend["entry_rate"] trend["employment_state"] = np.select( [trend.percent_change <= -5, trend.percent_change > 5], ["明顯下降", "明顯增加"], default="大致持平" ) trend["vacancy_pressure"] = (trend.vacancy_rate >= 3.1) | (trend.recruitment_months >= 3.5) def classify(row: pd.Series) -> str: code = row.industry_code if code == "agriculture": return "人力退場與接班風險（農業專案證據）" if code == "public_administration": return "資料範圍外／不判定" if row.percent_change <= -5: return "小基數專業替補風險" if row.recruitment_months >= 3.5 else "需求仍待確認的收縮" if row.percent_change > 5 and row.vacancy_pressure: if row.turnover_gap_pp >= 0.3: return "擴張伴隨留任壓力" return "需求擴張快於補實" if row.percent_change <= 5 and row.vacancy_pressure: return "存量持平但技能／補實壓力" if row.turnover_gap_pp >= 0.3: return "高流動留任壓力" return "本次未見強烈背離" trend["screening_type"] = trend.apply(classify, axis=1) trend = trend.sort_values("percent_change") trend.to_csv(OUT / "industry-screening.csv", index=False, encoding="utf-8-sig") selected = ["agriculture", "manufacturing", "construction", "transport_storage", "accommodation_food", "health_social"] indexed = emp[["year"] + selected].copy() for col in selected: indexed[col] = indexed[col] / indexed.loc[indexed.year == 2014, col].iloc[0] * 100 indexed.to_csv(OUT / "employment-index-2014-2025.csv", index=False, encoding="utf-8-sig") colors = {"brown": "#7a3f22", "orange": "#c97832", "green": "#2f6b4f", "blue": "#376a8a", "gray": "#8c8178", "cream": "#f6efe5", "red": "#a74c43"} # Figure 1: indexed employment paths. fig, ax = plt.subplots(figsize=(10.8, 6.2)) line_colors = [colors["brown"], colors["gray"], colors["orange"], colors["blue"], colors["red"], colors["green"]] for col, color in zip(selected, line_colors): ax.plot(indexed.year, indexed[col], marker="o", lw=2.1, ms=3.8, label=LABELS[col], color=color) ax.axhline(100, color="#bdb3aa", lw=1, ls="--") ax.set(title="圖1 2014–2025年六類產業就業人數指數", xlabel="年", ylabel="2014年＝100") ax.legend(ncol=3, frameon=False, loc="upper left") ax.grid(axis="y", color="#e5ddd5", lw=.7) fig.tight_layout(); fig.savefig(FIG / "figure-1-employment-index.png", bbox_inches="tight"); plt.close(fig) # Figure 2: endpoint changes across all broad industries. plot = trend.sort_values("percent_change") fig, ax = plt.subplots(figsize=(10.8, 7.2)) bar_colors = [colors["red"] if v <= -5 else colors["green"] if v > 5 else colors["orange"] for v in plot.percent_change] ax.barh(plot.industry, plot.percent_change, color=bar_colors) ax.axvline(-5, color="#8c8178", ls="--", lw=1); ax.axvline(5, color="#8c8178", ls="--", lw=1) ax.axvline(0, color="#33271f", lw=.8) for y, v in enumerate(plot.percent_change): ax.text(v + (0.7 if v >= 0 else -0.7), y, f"{v:+.1f}%", ha="left" if v >= 0 else "right", va="center", fontsize=8) ax.set(title="圖2 19類產業就業人數端點變化", xlabel="2014至2025年變化（%）", ylabel="") ax.grid(axis="x", color="#e5ddd5", lw=.7) fig.tight_layout(); fig.savefig(FIG / "figure-2-industry-change.png", bbox_inches="tight"); plt.close(fig) # Figure 3: vacancy rate and recruitment duration. joined = vac.merge(trend[["industry_code", "percent_change"]], left_on="employment_column", right_on="industry_code") fig, ax = plt.subplots(figsize=(10.8, 6.6)) sc = ax.scatter(joined.vacancy_rate, joined.recruitment_months, s=np.sqrt(joined.vacancies) * 5.2, c=joined.percent_change, cmap="RdYlGn", vmin=-10, vmax=30, alpha=.8, edgecolor="white", linewidth=.8) ax.axvline(3.1, color="#8c8178", ls="--", lw=1); ax.axhline(3.5, color="#8c8178", ls="--", lw=1) for _, r in joined.iterrows(): if r.vacancy_rate >= 4 or r.recruitment_months >= 4.1 or r.vacancies >= 16000: ax.annotate(LABELS[r.employment_column], (r.vacancy_rate, r.recruitment_months), xytext=(4, 4), textcoords="offset points", fontsize=8) ax.set(title="圖3 職缺率、招募時間與長期就業變化", xlabel="2025年3月底職缺率（%）", ylabel="全時職缺平均招募時間（月）") fig.colorbar(sc, ax=ax, label="2014–2025就業變化（%）") ax.grid(color="#e5ddd5", lw=.7) fig.tight_layout(); fig.savefig(FIG / "figure-3-vacancy-recruitment.png", bbox_inches="tight"); plt.close(fig) # Figure 4: employment change vs combined demand-pressure index. joined["demand_pressure_index"] = ((joined.vacancy_rate - 3.1) / joined.vacancy_rate.std(ddof=0) + (joined.recruitment_months - 3.5) / joined.recruitment_months.std(ddof=0)) / 2 fig, ax = plt.subplots(figsize=(10.8, 6.5)) ax.scatter(joined.percent_change, joined.demand_pressure_index, s=75, color=colors["brown"], alpha=.8) ax.axvline(-5, color="#8c8178", ls="--"); ax.axvline(5, color="#8c8178", ls="--"); ax.axhline(0, color="#8c8178", ls="--") for _, r in joined.iterrows(): ax.annotate(LABELS[r.employment_column], (r.percent_change, r.demand_pressure_index), xytext=(3, 3), textcoords="offset points", fontsize=7.5) ax.text(-24, 2.0, "人力下降＋補實壓力", color=colors["red"], fontsize=10) ax.text(12, 2.0, "需求擴張＋補實壓力", color=colors["green"], fontsize=10) ax.set(title="圖4 非農產業需求—勞動供給背離篩選", xlabel="2014–2025就業變化（%）", ylabel="需求壓力指數（職缺率與招募時間標準化平均）") ax.grid(color="#e5ddd5", lw=.7) fig.tight_layout(); fig.savefig(FIG / "figure-4-divergence-matrix.png", bbox_inches="tight"); plt.close(fig) # Figure 5: entry and exit rates. fig, ax = plt.subplots(figsize=(9.8, 6.5)) ax.scatter(vac.entry_rate, vac.exit_rate, s=np.sqrt(vac.vacancies) * 5, color=colors["orange"], alpha=.8, edgecolor="white") lim = max(vac.entry_rate.max(), vac.exit_rate.max()) + .3 ax.plot([0, lim], [0, lim], color="#8c8178", ls="--", lw=1) for _, r in vac.iterrows(): gap = r.exit_rate - r.entry_rate if abs(gap) >= .3 or r.vacancies >= 16000: ax.annotate(LABELS[r.employment_column], (r.entry_rate, r.exit_rate), xytext=(4, 3), textcoords="offset points", fontsize=8) ax.set(title="圖5 進入率與退出率：把留任問題與存量下降分開", xlabel="2025年3月進入率（%）", ylabel="2025年3月退出率（%）", xlim=(0, lim), ylim=(0, lim)) ax.grid(color="#e5ddd5", lw=.7) fig.tight_layout(); fig.savefig(FIG / "figure-5-entry-exit.png", bbox_inches="tight"); plt.close(fig) # Figure 6: agricultural succession risk. age = agr[agr.indicator == "manager_age"].copy() fig, ax1 = plt.subplots(figsize=(10.2, 6.1)) ax2 = ax1.twinx() ax1.plot(age.year, age.manager_mean_age, color=colors["brown"], marker="o", lw=2.5, label="平均年齡") ax2.bar(age.year, age.share_age_65_plus, width=2.6, color="#d9c7b5", alpha=.75, label="65歲以上占比") ax1.set(title="圖6 農業經營管理者高齡化與接班風險", xlabel="農林漁牧業普查年", ylabel="平均年齡（歲）") ax2.set_ylabel("65歲以上占比（%）") ax1.set_xticks(age.year) lines, labels = ax1.get_legend_handles_labels(); bars, labels2 = ax2.get_legend_handles_labels() ax1.legend(lines + bars, labels + labels2, frameon=False, loc="upper left") ax1.grid(axis="y", color="#e5ddd5", lw=.7) fig.tight_layout(); fig.savefig(FIG / "figure-6-agriculture-succession.png", bbox_inches="tight"); plt.close(fig) ag_start = float(emp.loc[emp.year == 2014, "agriculture"].iloc[0]) ag_end = float(emp.loc[emp.year == 2025, "agriculture"].iloc[0]) mining = trend.loc[trend.industry_code == "mining"].iloc[0] results = { "report_number": "SHWRP-2026-021", "period": [2014, 2025], "industry_count": len(LABELS), "vacancy_industry_count": len(vac), "employment_total_change_percent": float((emp.total.iloc[-1] / emp.total.iloc[0] - 1) * 100), "agriculture": { "employment_2014_thousand": ag_start, "employment_2025_thousand": ag_end, "change_thousand": ag_end - ag_start, "percent_change": float((ag_end / ag_start - 1) * 100), "manager_mean_age_2020": 64.13, "share_age_65_plus_2020": 48.76, "farm_sales_change_2015_2020_percent": 7.0, "rice_area_change_2015_2020_percent": 7.6, "screening_conclusion": "high-confidence labor-exit and succession-risk candidate; demand persistence is supported by separate agricultural indicators, not by the vacancy survey" }, "mining": { "employment_change_percent": float(mining.percent_change), "employment_2025_thousand": float(mining.employment_2025_thousand), "vacancy_rate": float(mining.vacancy_rate), "recruitment_months": float(mining.recruitment_months), "screening_conclusion": "small-base specialist replacement watchlist; broad demand persistence not established" }, "selected_expansion_pressure": trend[trend.screening_type.str.contains("需求擴張", na=False)][["industry", "percent_change", "vacancy_rate", "recruitment_months", "screening_type"]].to_dict("records"), "method_thresholds": {"employment_decline_percent": -5, "employment_growth_percent": 5, "vacancy_rate_percent": 3.1, "recruitment_months": 3.5}, "limitations": [ "Vacancy data are a March 2025 cross-section and exclude agriculture and public administration.", "Industry classifications changed across the period; the official linked series is used but discontinuity risk remains.", "Agricultural output value and sales are nominal and do not alone identify real demand or profitability.", "Broad-industry averages can hide occupation, region, contract, wage, and skill heterogeneity." ] } (ROOT / "analysis-results.json").write_text(json.dumps(results, ensure_ascii=False, indent=2), encoding="utf-8") print(json.dumps(results, ensure_ascii=False, indent=2))