做省际绿色全要素生产率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}{NMI}\left(\sum_{n1}^{N}\frac{s_n^x}{g_n^x} \sum_{m1}^{M}\frac{s_m^y}{g_m^y} \sum_{i1}^{I}\frac{s_i^b}{g_i^b}\right) ]约束条件为[ \sum_{k1}^{K}z_{k}x_{kn} s_n^x x_{kn},\quad n1,\dots,N ][ \sum_{k1}^{K}z_{k}y_{km} - s_m^y y_{km},\quad m1,\dots,M ][ \sum_{k1}^{K}z_{k}b_{ki} s_i^b b_{ki},\quad i1,\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^xx_k)、(g^yy_k)、(g^bb_k)这样松弛比例可以被解释为相对自身的改进空间。模型的目标函数越大说明这个省份相对于生产前沿面的改进空间越大无效率程度越高。因此GTFP效率值定义为[ GTFP 1 - \theta^* ]其中(\theta^*)是目标函数的最优值。GTFP的取值在0到1之间越接近1说明绿色发展效率越高。不同文献对权重系数的处理略有差异有的把三类松弛项分别取平均再加总有的像我在这里一样统一用(1/(NMI))加权。这个差异不影响省份之间的排序关系但论文里要和前文公式保持一致。2.3 省际GTFP的投入产出指标怎么搭指标选取是整个测算中最容易出问题的地方。一个省的数据指标选得不好模型算出来再漂亮也没有意义。以下是我建议的基础指标体系也是目前文献中使用频率最高的一套。指标类型变量名单位常用来源投入资本存量亿元2000年不变价基于永续盘存法自己测算投入劳动力万人《中国统计年鉴》投入能源消费总量万吨标准煤《中国能源统计年鉴》期望产出实际GDP亿元2000年不变价《中国统计年鉴》非期望产出CO2排放量万吨IPCC系数法自行估算或CEADs数据库资本存量最常用的做法是永续盘存法(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 python3.10 conda activate gtfp接着安装几个核心库pip install pandas numpy openpyxl pip install gurobipyGurobi的安装比一般库稍微麻烦一点因为它需要许可证。学术用户可以直接在官网申请免费学术许可然后在终端里用grbgetkey命令激活。激活完成后Python里执行import gurobipy不报错就说明环境通了。编辑器方面我用的是VSCode配合Python扩展。我的习惯是打开项目根目录后立刻用命令面板选择解释器CtrlShiftP输入Python: Select Interpreter选择刚才建好的gtfp环境。这一步非常关键选错了解释器会导致代码里安装的库全部找不到也会出现后面要讲到的Pylance解析路径报错。3.2 数据整理把省市年份指标做成标准的“长面板”环境配好之后真正花时间的其实是数据整理。我强烈建议把Excel数据整理成“长面板”格式也就是一列省份、一列年份、后面是各个指标每一行对应一个省份某一年份的观测。结构化之后大概长这样provinceyearcapitallaborenergygdpco2北京201012345.61200.54000.07890.02000.0北京201113209.81215.34120.58550.22050.0天津20108765.4780.23200.05200.01800.0整理数据时有几个容易忽略的细节。第一省份名称和编码一定要统一同一个省份不能有时叫“北京”、有时叫“北京市”。第二单位必须统一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, lb0, namesx) sy model.addVars(M, lb0, namesy) sb model.addVars(I, lb0, namesb) z model.addVars(K, lb0, namez) # 目标函数三类松弛比例的平均值 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(dropTrue) 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_indexTrue) result_df.to_csv(gtfp_results.csv, indexFalse, encodingutf-8-sig)跑完之后会得到一个类似下面的结果provinceyearGTFP北京20100.8734北京20110.8851天津20100.7912如果是用全局前沿测算也就是把所有年份的所有省份放在一起构建参照集只需要不按年份分组、直接用全量数据调用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测算常见问题速查表现象可能原因排查方向与解决建议所有效率值都是1决策单元数量过少、方向向量选择不当改用全局前沿检查方向向量是否过大在合理范围内调整模型维度某省效率值长期为0该省投入产出组合严重偏离前沿检查该省数据是否录入异常确认是否纳入所有应纳入的指标求解器报License错误Gurobi许可证未激活或已过期重新执行grbgetkey或改用cvxpyHiGHS开源方案VSCode代码标红但能运行Python解析器路径错误选择正确解释器删除无效extraPaths清空Pylance缓存同一批数据两次跑结果不一致缺失值或0值处理方式不统一将数据预处理步骤固定成脚本确保可复现结果表某年缺失该年某些省份数据为空导致模型求解失败检查该年的空值比例对相关样本做插补删除后再跑我个人在实际操作中的体会是GTFP测算真正考验人的不是模型推导而是对数据的态度。SBM-DDF本身是一个非常成熟的模型Python实现也不复杂但面对真实省际数据时你会遇到单位不统一、缺失值处理、异常年份、方向向量的选取这些从论文里根本学不到的细节问题。把这套流程固化下来后续再做任何类型的DEA分析都会顺畅很多。最后再分享一个实用的小技巧输出效率值之后先做一张按年份顺序排列的各省热力图可以非常快速地发现数据异常和模型结果异常这个习惯帮我避免过不止一次“跑完才发现某年数据带错了单位”的惨痛教训。