Python实现省际绿色全要素生产率测算：基于SBM-DDF模型的完整实操指南

做省际绿色全要素生产率（GTFP）测算，最近几年我在学术项目里几乎每周都要碰一次。单纯的TFP只考虑资本、劳动和产出，而GTFP还要把能源消耗和污染排放这类非期望产出拉进生产框架，这就得靠SBM-DDF模型来做前沿效率测算，再用Python批量求解。这篇文章把我用真实省级面板数据跑完整套流程的经验整理成一份能直接复用的实操记录，适合正在做绿色经济、区域发展或者论文实证的同行。内容尽量不堆公式，但关键数学形式和代码逻辑我都会讲透，保证你看完能自己算出一张分省分年的GTFP表。

1. GTFP测算的完整思路：为什么是SBM-DDF模型

1.1 一个省的经济表现，只看GDP远远不够

绿色全要素生产率的内核很简单：在资本、劳动、能源这些要素投入下，既要让期望产出（比如实际GDP）尽可能高，又要让非期望产出（比如二氧化碳排放）尽可能低。传统全要素生产率测的是“投入产出的转换效率”，GTFP则在这个基础上叠加了资源和环境的约束，所以它反映的是一个地区绿色发展的综合效率，而不是单纯的经济规模。

省际GTFP最常用的场景有三类。第一类是横向比较，看哪些省份的绿色发展效率全国领先，哪些省份属于高投入高排放的粗放模式。第二类是纵向动态分析，观察某个地区GTFP在时间维度的变化趋势，判断产业转型和污染治理是否真正见效。第三类是作为被解释变量，在后面接一系列影响因素回归，比如环境规制、产业结构、技术创新对绿色效率的作用方向。

做这类测算遇到的第一道坎就是模型选择。GTFP本质上是在“多投入、多产出、含坏产出”的框架下测算效率，普通的最小二乘法根本没法处理，直接用传统的DEA模型又没法把污染排放放进去，所以绝大多数研究都会选择非径向、非角度的SBM-DDF模型。

1.2 SBM-DDF模型为什么是省际GTFP的主流选择

传统DEA模型（CCR、BCC）有一个明显的缺陷：它是径向模型，假设所有投入或产出按同一比例缩放。可现实里哪个省份会把资本和劳动同时等比例减少？更麻烦的是，传统DEA没法把污染排放作为“坏产出”放进模型，顶多是把污染当作投入处理，这在经济学逻辑上很别扭——污染是生产过程的副产物，不是要素投入。

后来的方向性距离函数（DDF）解决了坏产出的问题，它通过设定方向向量，让模型沿着“期望产出增加、非期望产出减少”的方向搜索改进空间。但DDF本质上仍然是径向的，它忽略了投入产出变量中的松弛变量，导致效率值偏高，甚至会出现“看起来很有效率但其实还有大量浪费”的矛盾结果。

SBM-DDF把SBM的松弛变量思想嵌进DDF框架里。它不再要求所有变量同比例改进，而是允许每个投入、产出、非期望产出各自有不同的松弛改进空间，同时通过方向向量定义改进方向。这样得到的效率值更精细，对省际数据的区分度也更好。所以目前国内做省级GTFP的文章里，SBM-DDF是一个相当稳健的默认选项。

1.3 这套Python流程能解决什么问题

很多课题组还在用DEAP或MaxDEA这类软件做DEA测算，这些软件在简单模型下够用，但一旦涉及多省份多年份的面板数据，问题就来了：数据要反复复制粘贴，结果要手动合并，换个指标变量全流程重来。用Python做则可以做到从数据读取、模型求解到结果输出的全流程自动化。

我这套流程的核心是一个基于Gurobi求解器的Python函数，它可以对每个省份、每个年份形成一个独立的线性规划模型，批量求解几百个决策单元毫无压力。跑完之后直接输出一张标准的长表，包含省份、年份、GTFP效率值三个核心字段，后续不管做面板回归还是画图都能直接接入。整个过程不依赖任何商业DEA软件，只需要Python环境和一个免费学术许可的Gurobi求解器。

2. SBM-DDF模型的数学原型与指标设计

2.1 从“同比例缩放”到“带方向的松弛改进”

理解SBM-DDF不需要太高深的数学基础，关键是抓住几个递进的概念。

最原始的DEA模型假设一个地区的生产活动可以被其他地区的线性组合“参照”。如果当前地区的投入产出组合落在生产前沿面内部，说明它还有改进空间。CCR模型和BCC模型的改进思路是“等比例缩放”，比如所有投入同时减少10%，这是径向思想。

DDF在这个基础上加了一个方向向量：投入往什么方向减少、期望产出往什么方向增加、非期望产出往什么方向减少。你可以理解为给每个决策单元指定了一个“改进路径”，改进时并不强制等比例，而是沿着方向向量走。

SBM-DDF更进一步：它把每个变量的松弛量直接放进目标函数，松弛量越大说明改进空间越大，无效率程度越高。因为每个变量的松弛量可以自行变化，所以模型不再受“等比例”约束，这就是它比传统径向模型更精确的原因。

2.2 SBM-DDF的数学形式与方向向量选择

我采用的SBM-DDF是一个相对简洁、便于编程的形式。假设有K个决策单元（省份），每个决策单元有N种投入、M种期望产出、I种非期望产出。对第k个决策单元，模型求解以下最优化问题：

[ \max ; \frac{1}{N+M+I}\left(\sum_{n=1}^{N}\frac{s_n^x}{g_n^x} + \sum_{m=1}^{M}\frac{s_m^y}{g_m^y} + \sum_{i=1}^{I}\frac{s_i^b}{g_i^b}\right) ]

约束条件为：

[ \sum_{k'=1}^{K}z_{k'}x_{k'n} + s_n^x = x_{kn},\quad n=1,\dots,N ]

[ \sum_{k'=1}^{K}z_{k'}y_{k'm} - s_m^y = y_{km},\quad m=1,\dots,M ]

[ \sum_{k'=1}^{K}z_{k'}b_{k'i} + s_i^b = b_{ki},\quad i=1,\dots,I ]

[ z_{k'} \ge 0,\quad s_n^x \ge 0,\quad s_m^y \ge 0,\quad s_i^b \ge 0 ]

这里(s_n^x)表示第n种投入的松弛量，也就是投入冗余；(s_m^y)表示第m种期望产出的松弛量，也就是产出不足；(s_i^b)表示第i种非期望产出的松弛量，也就是污染冗余。方向向量(g)在代码里通常直接取当前决策单元的投入产出值，即(g^x=x_k)、(g^y=y_k)、(g^b=b_k)，这样松弛比例可以被解释为相对自身的改进空间。

模型的目标函数越大，说明这个省份相对于生产前沿面的改进空间越大，无效率程度越高。因此GTFP效率值定义为：

[ GTFP = 1 - \theta^* ]

其中(\theta^*)是目标函数的最优值。GTFP的取值在0到1之间，越接近1说明绿色发展效率越高。不同文献对权重系数的处理略有差异，有的把三类松弛项分别取平均再加总，有的像我在这里一样统一用(1/(N+M+I))加权。这个差异不影响省份之间的排序关系，但论文里要和前文公式保持一致。

2.3 省际GTFP的投入产出指标怎么搭

指标选取是整个测算中最容易出问题的地方。一个省的数据指标选得不好，模型算出来再漂亮也没有意义。以下是我建议的基础指标体系，也是目前文献中使用频率最高的一套。

资本存量最常用的做法是永续盘存法：(K_t = I_t + (1-\delta)K_{t-1})，折旧率取10.96%，基期资本存量可以从张军等文献的测算结果里找到。实际GDP需要按GDP平减指数换算成不变价，这样才能保证省际和跨年份可比。

CO2排放量我建议直接用IPCC推荐的方法估算：根据各省煤炭、原油、天然气的消费量和对应的碳排放系数计算。如果用CEADs数据库，可以直接拿到省际化石燃料燃烧的CO2排放数据，会省去很多重复计算的工作。要注意能源消费总量和CO2排放不能同时用“全国总量拆分”这种高度推算的数据，否则在做省际比较的时候会出现虚假差异。

3. Python实现全流程：环境、数据与代码

3.1 环境准备：Python安装与VSCode配置一步到位

很多刚接触的人第一步就卡在环境上。我的建议是不要手搓裸Python环境，直接用Anaconda创建虚拟环境，这样pandas、numpy这些常用库一次性装好。

具体步骤很简单。先从官网下载Anaconda安装包，装好之后打开终端，创建一个专门用于GTFP计算的虚拟环境：

conda create -n gtfp python=3.10 conda activate gtfp

接着安装几个核心库：

pip install pandas numpy openpyxl pip install gurobipy

Gurobi的安装比一般库稍微麻烦一点，因为它需要许可证。学术用户可以直接在官网申请免费学术许可，然后在终端里用grbgetkey命令激活。激活完成后，Python里执行import gurobipy不报错就说明环境通了。

编辑器方面我用的是VSCode配合Python扩展。我的习惯是打开项目根目录后立刻用命令面板选择解释器：Ctrl+Shift+P，输入Python: Select Interpreter，选择刚才建好的gtfp环境。这一步非常关键，选错了解释器会导致代码里安装的库全部找不到，也会出现后面要讲到的Pylance解析路径报错。

3.2 数据整理：把省市年份指标做成标准的“长面板”

环境配好之后，真正花时间的其实是数据整理。我强烈建议把Excel数据整理成“长面板”格式，也就是一列省份、一列年份、后面是各个指标，每一行对应一个省份某一年份的观测。

结构化之后大概长这样：

整理数据时有几个容易忽略的细节。第一，省份名称和编码一定要统一，同一个省份不能有时叫“北京”、有时叫“北京市”。第二，单位必须统一，GDP有的数据库是亿元、有的是万元，能源有的是万吨标准煤、有的是吨标准煤，不统一后面计算结果会离谱。第三，方向向量会出现除零风险，比如某些早期年份CO2排放量为0，或者某个省份某项指标缺失严重，这种数据要么插补，要么在计算时对该指标加一个极小正值。

对于缺失值，我一般先用线性插值处理量和漏掉的比例，如果某省份连续多年数据缺失，直接删除该省份的这几个年份，而不是硬插值。否则构造出的虚假数据会让生产前沿面失真。

3.3 核心代码：SBM-DDF求解函数与面板循环

数据准备好后，核心代码就是一个求解函数加一个面板循环。下面这个函数是整套流程的心脏，它接收一个决策单元集合的投入矩阵、期望产出矩阵和非期望产出矩阵，对每个省份构建一个线性规划模型并求解。

import numpy as np import pandas as pd import gurobipy as gp from gurobipy import GRB def solve_sbm_ddf(X, Y, B): K = X.shape[0] N = X.shape[1] M = Y.shape[1] I = B.shape[1] # 方向向量默认取当前决策单元的投入产出值 gx = X.values.copy() gy = Y.values.copy() gb = B.values.copy() scores = [] for k in range(K): model = gp.Model("dmu_%d" % k) model.setParam("OutputFlag", 0) # 不输出求解日志 # 定义松弛变量和强度变量 sx = model.addVars(N, lb=0, name="sx") sy = model.addVars(M, lb=0, name="sy") sb = model.addVars(I, lb=0, name="sb") z = model.addVars(K, lb=0, name="z") # 目标函数：三类松弛比例的平均值 obj = gp.LinExpr() for n in range(N): obj += sx[n] / gx[k, n] for m in range(M): obj += sy[m] / gy[k, m] for i in range(I): obj += sb[i] / gb[k, i] model.setObjective(obj / (N + M + I), GRB.MAXIMIZE) # 投入约束：参考组合加投入冗余等于当前投入 for n in range(N): lhs = gp.quicksum(z[k2] * X.iloc[k2, n] for k2 in range(K)) model.addConstr(lhs + sx[n] == X.iloc[k, n]) # 期望产出约束：参考组合减去产出不足等于当前产出 for m in range(M): lhs = gp.quicksum(z[k2] * Y.iloc[k2, m] for k2 in range(K)) model.addConstr(lhs - sy[m] == Y.iloc[k, m]) # 非期望产出约束：参考组合加污染冗余等于当前污染 for i in range(I): lhs = gp.quicksum(z[k2] * B.iloc[k2, i] for k2 in range(K)) model.addConstr(lhs + sb[i] == B.iloc[k, i]) model.optimize() if model.status == GRB.OPTIMAL: scores.append(1 - model.objVal) else: scores.append(np.nan) model.dispose() return np.array(scores)

这个函数里最容易搞错的就是约束条件的符号。投入约束的(s_x)表示冗余，所以参考组合加上冗余等于当前值；期望产出约束的(s_y)表示不足，参考组合减去不足等于当前值；非期望产出约束的(s_b)表示污染冗余，参考组合加上冗余等于当前值。符号反了求解出来的目标函数就会变成负数，效率值自然也不对。

下面是面板数据调用的完整主流程：

df = pd.read_excel("province_panel.xlsx") results = [] # 如果要做当期前沿，则按年份分组求解 for year in sorted(df["year"].unique()): sub = df[df["year"] == year].reset_index(drop=True) X = sub[["capital", "labor", "energy"]] Y = sub[["gdp"]] B = sub[["co2"]] eff = solve_sbm_ddf(X, Y, B) temp = pd.DataFrame({ "province": sub["province"], "year": year, "GTFP": eff }) results.append(temp) result_df = pd.concat(results, ignore_index=True) result_df.to_csv("gtfp_results.csv", index=False, encoding="utf-8-sig")

跑完之后会得到一个类似下面的结果：

如果是用全局前沿测算，也就是把所有年份的所有省份放在一起构建参照集，只需要不按年份分组、直接用全量数据调用solve_sbm_ddf即可。两种方式得到的结果含义不同，论文里写清楚用的是哪一种就好。我个人的习惯是：如果研究重点是GTFP增长率和GML指数，就用全局前沿；如果只是评价每年各省的效率排名，用当期前沿更稳妥。

4. 避坑指南：GTFP测算中的典型问题与解决实录

4.1 VSCode环境报错“cannot be resolved against python helper roots”的根治

这个报错我在第一次配置VSCode时踩过，排查了很久。它在代码里表现为Python插件的智能提示和跳转全部失效，某些内置函数也被标红，但代码实际运行时又能正常跑。核心原因是Pylance语言服务器解析Python辅助路径时出了问题，最常见的触发场景是：项目目录里存在多个虚拟环境，或者.vscode/settings.json里配置了不存在的python.analysis.extraPaths。

修复办法按顺序来。第一，在命令面板里执行Python: Select Interpreter，选择当前项目对应的虚拟环境解释器。第二，打开项目根目录下的.vscode/settings.json，如果看到类似这样的配置，检查路径是否存在：

{ "python.analysis.extraPaths": [ "C:/Users/xxx/anaconda3/envs/gtfp/Lib/site-packages" ] }

如果路径不存在或者和当前环境不匹配，直接删除整个settings.json，然后执行Developer: Reload Window。第三，如果报错还没消失，说明Pylance缓存损坏，需要在VSCode命令面板里执行Python: Clear Cache and Reload Window。做完这三步，这个报错基本不会再出现。

4.2 效率值全为1：先怀疑参照集，再怀疑方向向量

如果测出来的所有省份GTFP都是1，那不是皆大欢喜，而是模型出问题了。最常见的原因是当期参照集里决策单元数量太少。假设某年只有6个省份进入测算，而模型本身有3种投入、1种期望产出、1种非期望产出，加起来的维度让每个省份都容易构造出自己的前沿组合，结果大家都变成有效率单元。

我一般把决策单元数量控制在投入产出总维度的3倍以上。如果样本确实少，可以考虑改用全局前沿测算，把多年份的数据合并成一个参照集，这样决策单元数量成倍增加。另一个容易被忽略的原因是方向向量取值问题：如果把方向向量取得过大，松弛比例会整体变小，效率值趋近于1；取得过小，效率值又普遍偏低。如果采用的是一般文献里的做法，即方向向量取当前决策单元自身的投入和产出值，多数情况下效率值会落在0.6到1之间，有明显区分度。

4.3 Gurobi许可证缺失与开源求解器替代方案

如果没有申请到Gurobi学术许可证，运行代码时会直接报错：

GurobiError: License expired or invalid

两个解决办法。第一个是去官网申请学术许可，用学校邮箱一般几小时就能批下来，然后执行grbgetkey激活。第二个是用开源求解器替代，比如把线性规划部分改用cvxpy配合HiGHS求解器，代码思路几乎一致，只是变量和约束的写法从model.addVar变成cp.Variable，对熟悉Python的人反而更直观。

我个人的习惯是在正式测算时用Gurobi，因为它处理几百个线性规划模型时速度更快更稳定；但如果只是教学演示，或者别人要我快速看一下结果，我会直接写一个cvxpy版本，完全避开许可证问题。

4.4 0值、缺失值和异常值对模型的影响

数据里的0值是个大麻烦。方向向量直接取当前DMU的投入产出值时，如果某个省份某一年的CO2排放量为0，那么目标函数里s_b / g_b就变成除以0，整个模型直接报错。我处理时会在数据预处理阶段检查每个指标的最小值，如果发现0值，就给该指标加一个极小正值，比如1e-6，避免分母为零。

异常值也一样要警惕。某省某年GDP录入时多打了一个零，这个省份会变成异常前沿点，把其他省份的效率值整体压低。我会在跑模型之前先做一个简单的描述性统计，画出每个变量的箱线图，找出明显偏离均值的记录逐一核对原始数据。

4.5 GTFP测算常见问题速查表

我个人在实际操作中的体会是，GTFP测算真正考验人的不是模型推导，而是对数据的态度。SBM-DDF本身是一个非常成熟的模型，Python实现也不复杂，但面对真实省际数据时，你会遇到单位不统一、缺失值处理、异常年份、方向向量的选取这些从论文里根本学不到的细节问题。把这套流程固化下来，后续再做任何类型的DEA分析，都会顺畅很多。最后再分享一个实用的小技巧：输出效率值之后，先做一张按年份顺序排列的各省热力图，可以非常快速地发现数据异常和模型结果异常，这个习惯帮我避免过不止一次“跑完才发现某年数据带错了单位”的惨痛教训。