import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np

# 设置学术风格
plt.style.use('seaborn-v0_8-whitegrid')
sns.set_palette("husl")

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# ============ 图1: 数字普惠金融指数时空演变热力图 ============
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 子图1: 各省份数字金融指数热力图
provinces = ['北京', '上海', '广东', '浙江', '江苏', '四川']
years = [2015, 2018, 2021, 2023]
data_heatmap = np.array([
    [245, 312, 398, 456],
    [238, 298, 385, 442],
    [198, 256, 334, 398],
    [205, 268, 345, 412],
    [192, 248, 328, 389],
    [165, 218, 298, 356]
])

im = axes[0].imshow(data_heatmap, cmap='RdYlGn', aspect='auto')
axes[0].set_xticks(range(len(years)))
axes[0].set_xticklabels(years)
axes[0].set_yticks(range(len(provinces)))
axes[0].set_yticklabels(provinces)
axes[0].set_title('各省份数字普惠金融指数演变', fontweight='bold')

# 添加数值标签
for i in range(len(provinces)):
    for j in range(len(years)):
        text = axes[0].text(j, i, f'{data_heatmap[i,j]}',
                          ha="center", va="center", fontsize=9)

plt.colorbar(im, ax=axes[0], label='指数值')

# 子图2: 东中西部地区差异趋势
axes[1].plot(years, data_heatmap[0], marker='o', label='北京')
axes[1].plot(years, data_heatmap[4], marker='s', label='江苏')
axes[1].plot(years, data_heatmap[5], marker='^', label='四川')
axes[1].set_title('区域数字金融发展对比', fontweight='bold')
axes[1].set_xlabel('年份')
axes[1].set_ylabel('数字普惠金融指数')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('figure1_spatial_temporal.png', dpi=300)
print("✅ 图1保存成功!")

# ============ 图2: 回归系数森林图 ============
fig, ax = plt.subplots(figsize=(10, 6))

models = ['基准模型', '稳健性检验1', '稳健性检验2', '稳健性检验3']
coefficients = [0.152, 0.135, 0.168, 0.141]
ci_lower = [0.110, 0.095, 0.128, 0.105]
ci_upper = [0.194, 0.175, 0.208, 0.177]

y_pos = np.arange(len(models))
ax.errorbar(coefficients, y_pos,
            xerr=[[(c-l) for c,l in zip(coefficients, ci_lower)],
                  [(u-c) for c,u in zip(coefficients, ci_upper)]],
            fmt='o', markersize=10, capsize=5, linewidth=2.5,
            color='#007bff', ecolor='#6f42c1')

ax.set_yticks(y_pos)
ax.set_yticklabels(models)
ax.set_xlabel('系数值')
ax.set_title('数字普惠金融对家庭收入的影响：稳健性检验', fontweight='bold')
ax.axvline(x=0, color='black', linestyle='--', linewidth=1.5)
ax.grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.savefig('figure2_forest_plot.png', dpi=300)
print("✅ 图2保存成功!")

\documentclass[12pt,a4paper]{article}
\usepackage[UTF8]{ctex}  % 中文支持
\usepackage{amsmath,amssymb}  % 数学公式
\usepackage{booktabs}  % 三线表
\usepackage{graphicx}  % 图片
\usepackage{cite}  % 参考文献
\usepackage{hyperref}  % 超链接

\title{数字普惠金融对农村家庭消费的影响研究}
\author{张三 \\ 某某大学金融学院}
\date{\today}

\begin{document}

\maketitle

\section{引言}
数字普惠金融作为金融创新的重要形式，对促进农村消费升级具有重要意义\cite{zhang2023}。

\section{模型设定}
基准回归模型设定如下：

\begin{equation}
    \ln(Consumption_{it}) = \alpha_0 + \alpha_1 Digital_{it} + \alpha_2 X_{it} + \mu_i + \lambda_t + \varepsilon_{it}
    \label{eq:baseline}
\end{equation}

其中，下标$i$和$t$分别表示省份和年份；$Digital_{it}$为数字普惠金融指数；$X_{it}$为控制变量向量；$\mu_i$和$\lambda_t$分别为省份和年份固定效应。

\section{实证结果}
表\ref{tab:baseline}报告了基准回归结果。

\begin{table}[htbp]
\centering
\caption{基准回归结果}
\label{tab:baseline}
\begin{tabular}{lcccc}
\toprule
                    & \multicolumn{2}{c}{固定效应模型}        & \multicolumn{2}{c}{随机效应模型}    \\
\cmidrule(lr){2-3} \cmidrule(lr){4-5}
变量                 & 系数          & 标准误         & 系数          & 标准误         \\
\midrule
数字普惠金融指数     & 0.152***      & (0.021)        & 0.148***      & (0.019)        \\
人均GDP（对数）      & 0.089**       & (0.035)        & 0.092**       & (0.033)        \\
受教育年限           & 0.045***      & (0.012)        & 0.043***      & (0.011)        \\
\midrule
省份固定效应         & Yes           &                & No            &                \\
年份固定效应         & Yes           &                & Yes           &                \\
观测数               & 12,450        &                & 12,450        &                \\
$R^2$                & 0.312         &                & 0.298         &                \\
\bottomrule
\end{tabular}
\medskip
\\ \footnotesize 注：***、**、*分别表示在1%、5%、10%水平上显著。
\end{table}

从表\ref{tab:baseline}可以看出，数字普惠金融指数的系数为0.152，在1\%水平上显著，表明数字普惠金融每提高1个单位，农村家庭消费支出增加约15.2\%。

\bibliographystyle{plain}
\bibliography{references}

\end{document}

```mermaid
graph TD
    A["数据来源"] --> A1["CFPS 中国家庭追踪调查"]
    A --> A2["CHFS 中国家庭金融调查"]
    A --> A3["北大数字普惠金融指数"]
    A --> A4["统计年鉴宏观数据"]
    A1 --> B["数据清洗与处理"]
    A2 --> B
    A3 --> B
    A4 --> B
    B --> B1["缺失值处理（删除/插值）"]
    B --> B2["异常值处理（IQR方法）"]
    B --> B3["数据合并（面板数据构建）"]
    B1 --> C["变量选择"]
    B2 --> C
    B3 --> C
    C --> C1["核心变量：数字普惠金融指数"]
    C --> C2["被解释变量：家庭收入/消费"]
    C --> C3["控制变量：GDP/教育/年龄"]
    C --> C4["机制变量：创业/信贷可得性"]
    C1 --> D["模型构建"]
    C2 --> D
    C3 --> D
    C4 --> D
    D --> D1["基准回归（OLS/固定效应）"]
    D --> D2["工具变量法（IV/2SLS）"]
    D --> D3["双重差分（DID）"]
    D --> D4["中介效应模型"]
    D1 --> E["实证分析"]
    D2 --> E
    D3 --> E
    D4 --> E
    E --> E1["描述性统计（均值/标准差/相关系数）"]
    E --> E2["基准回归（系数显著性）"]
    E --> E3["平行趋势检验（DID前）"]
    E1 --> F{"稳健性检验"}
    E2 --> F
    E3 --> F
    F --> F1["替换变量（滞后一期）"]
    F --> F2["改变样本（剔除直辖市）"]
    F --> F3["PSM 倾向得分匹配"]
    F --> F4["安慰剂检验（随机化处理）"]
    F1 --> G["异质性分析"]
    F2 --> G
    F3 --> G
    F4 --> G
    G --> G1["城乡差异"]
    G --> G2["区域差异（东中西部）"]
    G --> G3["收入分组（低中高）"]
    G --> G4["教育水平"]
    G1 --> H["机制检验"]
    G2 --> H
    G3 --> H
    G4 --> H
    H --> H1["中介效应（Bootstrap检验）"]
    H --> H2["调节效应（交互项分析）"]
    H --> H3["渠道分析（分步回归）"]
    H1 --> I["研究结论与政策建议"]
    H2 --> I
    H3 --> I
    I --> I1["核心结论"]
    I --> I2["政策建议"]
    I --> I3["研究局限"]
    I --> I4["未来研究方向"]
    classDef dataSource fill:#e3f2fd,stroke:#1976d2,stroke-width:2px
    classDef dataProcess fill:#bbdefb,stroke:#1565c0,stroke-width:2px
    classDef variable fill:#f3e5f5,stroke:#7b1fa2,stroke-width:2px
    classDef model fill:#e8f5e9,stroke:#388e3c,stroke-width:2px
    classDef analysis fill:#fff3e0,stroke:#f57c00,stroke-width:2px
    classDef robust fill:#fce4ec,stroke:#c2185b,stroke-width:2px
    classDef hetero fill:#e0f2f1,stroke:#00796b,stroke-width:2px
    classDef mechanism fill:#f1f8e9,stroke:#558b2f,stroke-width:2px
    classDef conclusion fill:#e8eaf6,stroke:#4527a0,stroke-width:2px
    class A,A1,A2,A3,A4 dataSource
    class B,B1,B2,B3 dataProcess
    class C,C1,C2,C3,C4 variable
    class D,D1,D2,D3,D4 model
    class E,E1,E2,E3 analysis
    class F,F1,F2,F3,F4 robust
    class G,G1,G2,G3,G4 hetero
    class H,H1,H2,H3 mechanism
    class I,I1,I2,I3,I4 conclusion
```

import pandas as pd
import stata_setup

# ============ 配置Stata路径 ============
# 首次使用需要配置Stata安装路径
stata_setup.config("C:/Program Files/Stata17", "se")

# 导入Stata API
from pystata import stata

# ============ 加载数据到Stata ============
print("📂 加载数据到Stata...")
df = pd.read_csv('data_cleaned.csv')

# 将Pandas DataFrame导入Stata
stata.df_store(df, "mydata")

# ============ 在Stata中执行命令 ============
print("\n🔬 执行固定效应回归...")

# 设置面板数据
stata.run("tsset province_id year")

# 固定效应回归
stata.run("""
    xtreg household_income digital_index gdp_per_capita education, fe
        cluster(province)
""", quietly=False)

# 保存回归结果
stata.run("estimates store fe_model")

# ============ 生成三线表（esttab） ============
print("\n📊 生成LaTeX格式结果表...")
stata.run("""
    esttab fe_model using "regression_results.tex", ///
        replace b(3) se(3) star(* 0.1 ** 0.05 *** 0.01) ///
        scalars("r2_a R-squared N N") ///
        title("数字普惠金融对家庭收入的影响")
""")

# ============ 导出结果到Python ============
# 获取回归系数
coef_df = stata.get_df()
print("\n✅ 回归完成！")
print(coef_df.head())

# ============ 批量回归示例 ============
print("\n🔄 批量回归不同模型...")
models = [
    "household_income digital_index",
    "household_income digital_index gdp_per_capita",
    "household_income digital_index gdp_per_capita education",
    "household_income digital_index gdp_per_capita education i.province i.year"
]

for i, model in enumerate(models):
    print(f"模型 {i+1}: {model}")
    stata.run(f"reg {model}")
    stata.run(f"estimates store model{i+1}")

# 比较所有模型
stata.run("esttab model* using "all_models.tex", replace star(* 0.1 ** 0.05 *** 0.01)")
print("\n✅ 批量回归完成！结果已保存。")