名师讲堂｜使用 Python 计算上市公司供应链地理加权距离

4033 字 约 13 分钟 名师讲堂 距离与边界上市公司Python

插图

由于借助 AI 工具学习编程已经变得非常容易了，因此之后的课程就不再默认进行视频讲解了，如果特别需要视频讲解也可以联系李老师预约讲解～讲义材料学习过程中遇到的问题也可以及时与李老师联系。

购买 RStata 名师讲堂会员即可参加该课程啦（之前的和未来的都可以参加）！

价格：2800/年 或者 4800/长期

购买会员可以从这里下单：https://rstata.duanshu.com/#/card/list/

名师讲堂会员权益：

  1. 参加每个月 3～4 次的名师讲堂课程；
  2. 参加平台上的其他 R 语言和 Stata 的课程；
  3. 以会员折扣价购买我们分享的数据资料（10 元/份）；
  4. 课程内外的提问解答服务（课程外的尽量帮忙解决）。

* 如果发票可添加小编微信 r_stata2 （RStata 李老师）开具。如需数据资料，购买后可添加小编微信免费领取数据折扣卡。

更多关于 RStata 会员的更多信息可添加微信号 r_stata2 咨询：

课程主页（点击文末的阅读原文即可跳转）：<>

指标来源

供应链地理加权距离指标来自邹颖、石福安、祁亚发表在《世界经济》2026 年第 1 期的论文《以数促联：公共数据开放与企业供应链地理布局》。

该论文采用堆叠双重差分模型考察公共数据开放对企业供应链地理布局的影响，其中核心被解释变量即为供应链地理加权距离，包括两个维度：

  • 供应商加权距离（Disws）：衡量上市企业与主要供应商之间的地理距离（加权）；
  • 客户加权距离（Diswc）：衡量上市企业与主要客户之间的地理距离（加权）。

论文发现，公共数据开放能打破地理距离约束，拓宽企业供应链分布范围。

指标定义与计算公式

数据来源

该指标的测算需要三类数据：

  1. 上市公司注册地址与办公地址数据：包含 2000～2024 年所有沪深 A 股上市公司的注册地址和办公地址经纬度信息，以及所处省市区县。
  2. 上市公司前 5 大供应商数据：包含 2001～2024 年上市公司前 5 大供应商的工商注册信息匹配结果，含供应商经纬度、采购额及采购额占比。
  3. 上市公司前 5 大客户数据：包含 2001～2024 年上市公司前 5 大客户的工商注册信息匹配结果，含客户经纬度、销售额及销售额占比。

计算公式

地理距离计算方法

地理距离采用 Haversine 大圆距离公式计算，与 R 语言 sf::st_distance(crs=4326) 使用完全相同的公式，保证了结果的可复现性。Python 实现如下：

import numpy as np

def haversine_km(lat1, lon1, lat2, lon2):
    """
    计算两点之间的 Haversine 大圆距离（千米）。
    参数为 pandas Series，支持向量化运算。
    """
    R = 6371.0  # 地球平均半径（千米）
    lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = np.sin(dlat / 2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2
    c = 2 * np.arcsin(np.sqrt(a))
    return R * c

其中：

  • R = 6371.0 为地球平均半径（千米）；
  • 使用 numpy 的向量化运算，支持对整列数据批量计算；
  • 最终结果单位为千米（km）。

上市公司地址使用的是办公地址（而非注册地址），这与论文中基于实际经营地址的选择一致。

加权方式说明

权重采用的是供应商采购额占前五大供应商采购总额的比例（而非供应商采购额占企业总采购额的比例）。这是因为上市公司年报中通常只披露前五大供应商/客户的名称和交易金额，无法获取完整的采购/销售总额数据。

计算过程

第 1 步：参数设置与数据读取

首先设置数据路径并读取三类数据：

import pandas as pd
import os

dir_path = "/Users/ac/Desktop/使用 Python 计算上市公司供应链地理加权距离"

gys_file  = os.path.join(dir_path,
    "2001~2024年上市公司前5大供应商工商注册信息匹配结果（含经纬度及所处的省市区县）.dta")
kh_file   = os.path.join(dir_path,
    "2001~2024年上市公司前5大客户工商注册信息匹配结果（含经纬度及所处的省市区县）.dta")
addr_file = os.path.join(dir_path,
    "2000～2024年上市公司注册地址与办公地址（含经纬度、所处的省市区县及搬迁距离）.dta")
outfile   = os.path.join(dir_path,
    "2001～2024年上市公司供应链地理加权距离_Python.dta")

# 读取数据
gys = pd.read_stata(gys_file)
kh  = pd.read_stata(kh_file)
addr = pd.read_stata(addr_file)
print(f"供应商记录数: {len(gys)}, 客户记录数: {len(kh)}, 地址记录数: {len(addr)}")

# 准备上市公司办公地址
addr_sub = addr[["股票代码", "年份", "办公地址_经度", "办公地址_纬度"]].copy()
addr_sub = addr_sub.rename(columns={"年份": "统计年份"})
addr_sub = addr_sub.dropna(subset=["统计年份", "办公地址_经度", "办公地址_纬度"])
print(f"有效地址记录: {len(addr_sub)}")

说明：使用 pandas.read_stata() 直接读取 Stata 的 .dta 格式文件，无需格式转换。dropna() 用于过滤经纬度缺失的记录。

第 2 步：计算供应商加权距离

首先读取供应商数据，筛选有效记录（年报、有经纬度、有采购额），然后合并办公地址并计算距离：

# 此处代码需下载讲义材料查看～

# 按公司-年度汇总
gys_weighted = gys_matched.groupby(["股票代码", "统计年份"]).agg(
    Disw_s=("weighted_dist", lambda x: np.log(1 + x.sum())),
    供应商有效数=("ratio_s", "count")
).reset_index()
print(f"Disw_s 计算完成: {len(gys_weighted)} 个公司-年度观测值")

关键逻辑说明：haversine_km 的参数顺序是 (纬度, 经度, 纬度, 经度)，注意不要写反经纬度的位置。groupby().transform("sum") 计算组内总额，再除以总额得到权重比例。lambda x: np.log(1 + x.sum()) 在 agg() 中直接计算加权距离的对数值。

第 3 步：计算客户加权距离

客户加权距离的计算逻辑与供应商完全对称，只是将”供应商采购额”替换为”客户销售额”：

# 筛选有效记录
kh_work = kh[
    (kh["报表类型"] == 1) &
    kh["经度"].notna() & kh["纬度"].notna() &
    kh["客户销售额"].notna() & kh["客户销售额占比"].notna()
].copy()
print(f"有效客户记录（年报+有经纬度+有销售额）: {len(kh_work)} / {len(kh)}")

kh_work = kh_work.merge(addr_sub, on=["股票代码", "统计年份"], how="left")
kh_matched = kh_work.dropna(subset=["办公地址_经度", "办公地址_纬度"])
print(f"成功匹配公司地址的记录: {len(kh_matched)}")

# 计算距离
kh_matched = kh_matched.copy()
kh_matched["dist_km"] = haversine_km(
    kh_matched["办公地址_纬度"], kh_matched["办公地址_经度"],
    kh_matched["纬度"], kh_matched["经度"]
)

# 计算加权距离
此处代码需下载讲义材料查看～

第 4 步：合并结果并输出

将供应商加权距离和客户加权距离合并到上市公司信息表中：

# 合并公司信息
info = addr[["股票代码", "年份", "股票简称", "行业代码C",
             "办公地址_经度", "办公地址_纬度", "办公地址_省", "办公地址_市"]].copy()
info = info.rename(columns={"年份": "统计年份"})

result = info.merge(gys_weighted, on=["股票代码", "统计年份"], how="left")
result = result.merge(kh_weighted, on=["股票代码", "统计年份"], how="left")
result = result.sort_values(["股票代码", "统计年份"]).reset_index(drop=True)
print(f"最终样本: {len(result)} 个公司-年度观测值")

# 去除 2000 年数据
result = result[result["统计年份"] >= 2001].copy()

# 保存（英文列名版本，pandas.to_stata 对中文列名有限制）
result_en = result.rename(columns={
    "股票代码": "stkcd", "统计年份": "year", "股票简称": "stkname",
    "行业代码C": "indC", "办公地址_经度": "off_lon", "办公地址_纬度": "off_lat",
    "办公地址_省": "off_prov", "办公地址_市": "off_city",
    "供应商有效数": "n_supplier", "客户有效数": "n_customer",
})
result_en.to_stata(outfile, write_index=False, version=118)

# 同时保存 CSV（保留中文列名）
csv_file = outfile.replace(".dta", ".csv")
result.to_csv(csv_file, index=False, encoding="utf-8-sig")
print(f"结果已保存至: {outfile}")
print(f"CSV 版本: {csv_file}")

注意：pandas.to_stata() 对中文列名有严格限制（仅支持 ASCII 字符），因此保存 .dta 文件时使用英文列名。如需保留中文列名，可同时保存 CSV 版本。

描述性统计

# 描述性统计
print("--- Disw_s（供应商加权距离）---")
ds = result["Disw_s"].dropna()
print(f"  N = {len(ds)}")
print(f"  均值 = {ds.mean():.3f}, 中位数 = {ds.median():.3f}, 标准差 = {ds.std():.3f}")
print(f"  最小值 = {ds.min():.3f}, 最大值 = {ds.max():.3f}")

print("--- Disw_c（客户加权距离）---")
dc = result["Disw_c"].dropna()
print(f"  N = {len(dc)}")
print(f"  均值 = {dc.mean():.3f}, 中位数 = {dc.median():.3f}, 标准差 = {dc.std():.3f}")
print(f"  最小值 = {dc.min():.3f}, 最大值 = {dc.max():.3f}")

可以看到：

  • 供应商加权距离（Disws）和客户加权距离（Diswc）的中位数分别为 6.197 和 6.408，均值分别为 5.840 和 5.915；
  • 由于取了对数变换，中位数大于均值说明存在明显的左偏分布；
  • 客户加权距离的标准差（1.504）略大于供应商加权距离（1.361），说明不同企业与客户的地理距离差异更大。

说明：上述描述性统计基于全样本（2001～2024 年）计算，与论文中 2007～2022 年的子样本统计量存在差异属于正常现象。

使用 reticulate 创建与管理 Python 虚拟环境

在 R 中通过 reticulate 包来调用 Python，最好的实践是为项目创建一个专属的 Python 虚拟环境，将所需依赖隔离到独立空间，避免与系统 Python（如 Anaconda）发生版本冲突。

重要说明（避免”已初始化”报错）：reticulate 在 R 会话中只能绑定一次 Python——一旦某个 {python} 代码块运行，Python 解释器就被锁定，之后再调用 use_virtualenv() 会报错：

ERROR: The requested version of Python cannot be used, as another version has already been initialized.

因此，虚拟环境的激活必须在所有 {python} 代码块之前完成。本文档的解决方案是在 setup chunk 中通过 Sys.setenv(RETICULATE_PYTHON = ...) 提前锁定 Python 路径，这是 reticulate 选取 Python 的最高优先级入口。

安装 reticulate（仅首次）

# 设置 CRAN 镜像（knit 时 R 处于非交互模式，不会自动选择镜像）
options(repos = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))

# 仅在尚未安装时才安装，避免每次 knit 都重装
if (!requireNamespace("reticulate", quietly = TRUE)) {
  install.packages("reticulate")
  message("reticulate 安装完成！")
} else {
  message("reticulate 已安装，版本：", packageVersion("reticulate"))
}

虚拟环境初始化原理（已在 setup chunk 中完成）

本文档的 setup chunk（隐藏运行）包含如下逻辑：

library(reticulate)

.venv_name   <- ".venv"
.venv_python <- virtualenv_python(.venv_name)

# 虚拟环境不存在时自动创建
if (!file.exists(.venv_python)) {
  virtualenv_create(.venv_name)
  .venv_python <- virtualenv_python(.venv_name)
}

# 通过环境变量抢先锁定 Python（优先级最高，早于任何 {python} chunk）
Sys.setenv(RETICULATE_PYTHON = .venv_python)
use_virtualenv(.venv_name, required = TRUE)

这样做的关键在于：knitr 在处理第一个 {python} chunk 时，reticulate 已经通过 RETICULATE_PYTHON 环境变量知道要使用 .venv，不会再去碰 Anaconda。

在虚拟环境中安装 Python 包（仅首次）

# 检查关键包是否已安装，缺失的才安装
py_pkgs <- c("numpy", "pandas")
installed <- py_list_packages(".venv")$package
need_install <- setdiff(py_pkgs, installed)

if (length(need_install) > 0) {
  virtualenv_install(".venv", packages = need_install)
  message("已安装缺失的包：", paste(need_install, collapse = ", "))
} else {
  message("所有 Python 包已就绪，无需安装")
}

验证激活状态

# 验证当前绑定的 Python 路径（应指向 .venv 目录）
py_config()

虚拟环境管理常用命令

# 查看所有已创建的虚拟环境
virtualenv_list()

# 删除虚拟环境（当不再需要时）
# virtualenv_remove(".venv")

# 升级某个包
# virtualenv_install(".venv", packages = "pandas", ignore_installed = TRUE)

完整 Python 脚本

# -*- coding: utf-8 -*-
# ============================================================
# 计算上市公司与供应商/客户的加权地理距离（Python 版本）
# 参考论文：以数促联：公共数据开放与企业供应链地理布局（邹颖等, 2026）
#
# 公式：
#   Disw_s_it = ln(1 + Σ Dis_s_ipt × Ratio_s_ipt)
#   Disw_c_it = ln(1 + Σ Dis_c_iqt × Ratio_c_iqt)
#
# 距离计算：Haversine 大圆距离（与 R sf::st_distance 一致）
# 地址选择：上市公司使用办公地址
# 数据计算：微信公众号 RStata
# ============================================================

import pandas as pd
import numpy as np
from math import radians, sin, cos, sqrt, atan2, log
import os

# ---- 1. 参数设置 ----
dir_path = "/Users/ac/Desktop/使用 Python 计算上市公司供应链地理加权距离"

gys_file  = os.path.join(dir_path,
    "2001~2024年上市公司前5大供应商工商注册信息匹配结果（含经纬度及所处的省市区县）.dta")
kh_file   = os.path.join(dir_path,
    "2001~2024年上市公司前5大客户工商注册信息匹配结果（含经纬度及所处的省市区县）.dta")
addr_file = os.path.join(dir_path,
    "2000～2024年上市公司注册地址与办公地址（含经纬度、所处的省市区县及搬迁距离）.dta")
outfile   = os.path.join(dir_path,
    "2001～2024年上市公司供应链地理加权距离_Python.dta")

# ---- 2. Haversine 大圆距离计算函数 ----
def haversine_km(lat1, lon1, lat2, lon2):
    """
    计算两点之间的 Haversine 大圆距离（千米）。
    参数为 pandas Series，支持向量化运算。
    """
    R = 6371.0  # 地球平均半径（千米）
    lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = np.sin(dlat / 2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2
    c = 2 * np.arcsin(np.sqrt(a))
    return R * c

# ---- 3. 读取数据 ----
print("读取供应商数据...")
gys = pd.read_stata(gys_file)
print(f"  供应商记录数: {len(gys)}, 唯一公司数: {gys['股票代码'].nunique()}")

print("读取客户数据...")
kh = pd.read_stata(kh_file)
print(f"  客户记录数: {len(kh)}, 唯一公司数: {kh['股票代码'].nunique()}")

print("读取上市公司地址数据...")
addr = pd.read_stata(addr_file)
print(f"  地址记录数: {len(addr)}, 唯一公司数: {addr['股票代码'].nunique()}")

# ---- 4. 准备上市公司办公地址 ----
print("准备上市公司办公地址...")
addr_sub = addr[["股票代码", "年份", "办公地址_经度", "办公地址_纬度"]].copy()
addr_sub = addr_sub.rename(columns={"年份": "统计年份"})
addr_sub = addr_sub.dropna(subset=["统计年份", "办公地址_经度", "办公地址_纬度"])
print(f"  有效地址记录: {len(addr_sub)}")

# ---- 5. 计算供应商加权距离 ----
# 此处代码需下载讲义材料查看～

# ---- 6. 计算客户加权距离 ----
此处代码需下载讲义材料查看～

# ---- 7. 合并结果 ----
print("\n========== 合并结果 ==========")

info = addr[["股票代码", "年份", "股票简称", "行业代码C",
             "办公地址_经度", "办公地址_纬度", "办公地址_省", "办公地址_市"]].copy()
info = info.rename(columns={"年份": "统计年份"})

result = info.merge(gys_weighted, on=["股票代码", "统计年份"], how="left")
result = result.merge(kh_weighted, on=["股票代码", "统计年份"], how="left")
result = result.sort_values(["股票代码", "统计年份"]).reset_index(drop=True)

print(f"  最终样本: {len(result)} 个公司-年度观测值")
print(f"  唯一公司数: {result['股票代码'].nunique()}")
print(f"  年份范围: {int(result['统计年份'].min())} - {int(result['统计年份'].max())}")

# ---- 8. 描述性统计 ----
print("\n========== 描述性统计 ==========")
print("--- Disw_s（供应商加权距离）---")
ds = result["Disw_s"].dropna()
print(f"  N = {len(ds)}")
print(f"  均值 = {ds.mean():.3f}, 中位数 = {ds.median():.3f}, 标准差 = {ds.std():.3f}")
print(f"  最小值 = {ds.min():.3f}, 最大值 = {ds.max():.3f}")

print("--- Disw_c（客户加权距离）---")
dc = result["Disw_c"].dropna()
print(f"  N = {len(dc)}")
print(f"  均值 = {dc.mean():.3f}, 中位数 = {dc.median():.3f}, 标准差 = {dc.std():.3f}")
print(f"  最小值 = {dc.min():.3f}, 最大值 = {dc.max():.3f}")

# ---- 9. 去除 2000 年数据 ----
result = result[result["统计年份"] >= 2001].copy()

# ---- 10. 保存结果 ----
result_en = result.rename(columns={
    "股票代码": "stkcd", "统计年份": "year", "股票简称": "stkname",
    "行业代码C": "indC", "办公地址_经度": "off_lon", "办公地址_纬度": "off_lat",
    "办公地址_省": "off_prov", "办公地址_市": "off_city",
    "供应商有效数": "n_supplier", "客户有效数": "n_customer",
})
result_en.to_stata(outfile, write_index=False, version=118)

csv_file = outfile.replace(".dta", ".csv")
result.to_csv(csv_file, index=False, encoding="utf-8-sig")
print(f"\n结果已保存至: {outfile}")
print(f"CSV 版本: {csv_file}")
print("完成！")

参考文献

邹颖、石福安、祁亚，2026：《以数促联：公共数据开放与企业供应链地理布局》，《世界经济》第 1 期。

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜使用 Python 计算上市公司供应链地理加权距离

相关文章

评论