SBM-GML指数实战：从模型构造到Python代码实现与避坑指南

简介本资源面向需要测算全要素生产率的科研人员与研究生提供基于Matlab的SBM-GML指数、ML指数及超效率SBM完整代码包可计算VRS与CRS下含非期望产出的效率值并依据投入产出数据生成GML指数。压缩包共16个文件约6.17MB包含11个m脚本文件、4份pdf说明文档和1份xlsx示例数据脚本按模型模块拆分文档覆盖理论介绍、Matlab安装、整体操作步骤与指标解读图形展示。已有5472人学习下载。资源对代码与结果均做了图文注释配套示例数据可跟随说明完整跑通流程帮助不熟悉Matlab的用户快速上手同时梳理了SBM-GML、GML-DDF与SBM-DDF三种GML计算方法的差异便于读者理解模型选择依据并复现结果。1. SBM-GML指数到底在算什么从ML指数失效到绿色全要素生产率的测度逻辑如果你跑过绿色全要素生产率GTFP的测算大概率遇到过这样的场景用传统的Malmquist-LuenbergerML指数算出来的结果要么线性规划无解要么在跨期方向距离函数下出现技术倒退的假象要么效率值大于1让人怀疑人生。SBM-GML指数就是在这个背景下被广泛采用的替代方案——它把非径向、非角度的SBM方向距离函数和Global Malmquist-Luenberger指数结合起来既处理了投入产出的松弛问题又规避了传统ML指数在跨期比较中的不可行解和不可传递性。这套方法在环境经济学、区域绿色发展评估、碳排放效率测算里已经是主流工具但真正动手算的时候数据包络分析DEA的建模细节、方向向量的设定、全局前沿的构造方式每一步都会直接影响最终结果的可靠性。这篇文章面向的是需要自己动手算出SBM-GML指数、并且希望结果经得起审稿人推敲的研究生和青年学者我会把从模型设定到代码实现再到结果验证的完整路径拆开讲清楚。2. SBM-GML指数的模型构造方向距离函数、全局前沿与指数分解2.1 为什么传统ML指数在环境约束下容易翻车传统Malmquist指数基于Shephard距离函数要求投入产出同比例变化这在处理非期望产出如CO₂、SO₂时非常别扭。Chung等人在1997年提出Malmquist-Lenberger指数引入方向距离函数DDF允许在增加期望产出的同时减少非期望产出。但ML指数有两个硬伤第一它使用当期前沿contemporaneous frontier跨期比较时可能出现线性规划无可行解第二ML指数不满足传递性也就是说从t到t1的指数乘以t1到t2的指数不等于t到t2的指数。这两个问题在面板数据较长时尤其明显很多论文里出现的效率值异常波动根源就在这里。GML指数Global Malmquist-Lenberger的核心改进是构造一个全局前沿global frontier即所有时期的决策单元DMU共同构成一个生产可能集。这样一来所有DMU都在同一个前沿面上比较既避免了不可行解又天然满足传递性。Oh2010证明了GML指数可以分解为效率变化EC和技术变化TC且EC和TC的乘积严格等于GML指数。这个性质在实证分析中非常重要因为你可以清楚地看到效率提升和技术进步各自的贡献。2.2 SBM方向距离函数的设定与方向向量选择SBMSlack-Based Measure的优势在于把投入和产出的松弛量直接纳入目标函数而不是像径向模型那样只考虑比例改进。结合方向距离函数后SBM-DDF的形式如下假设有n个DMU每个DMU有m种投入x、s种期望产出y、k种非期望产出b。方向向量g(g_x, g_y, g_b)通常取g(-x, y, -b)表示投入减少、期望产出增加、非期望产出减少的方向。SBM-DDF的效率值通过求解以下线性规划得到import numpy as np from scipy.optimize import linprog def sbm_ddf(x, y, b, X, Y, B, gx, gy, gb): 求解单个DMU的SBM方向距离函数值 x, y, b: 当前DMU的投入、期望产出、非期望产出向量 X, Y, B: 所有DMU的投入、期望产出、非期望产出矩阵 gx, gy, gb: 方向向量 返回: beta值效率损失程度beta越小越有效 n X.shape[0] m X.shape[1] s Y.shape[1] k B.shape[1] # 决策变量: [lambda_1...lambda_n, beta, sx_1...sx_m, sy_1...sy_s, sb_1...sb_k] # 目标: min beta - (1/(msk)) * (sum(sx/x) sum(sy/y) sum(sb/b)) # 简化版: 只最小化beta松弛量通过约束体现 c np.zeros(n 1 m s k) c[n] 1 # beta系数 # 等式约束: 投入、期望产出、非期望产出的前沿构造 A_eq [] b_eq [] # 投入约束: X^T * lambda sx x - beta * gx for i in range(m): row np.zeros(n 1 m s k) row[:n] X[:, i] row[n] gx[i] row[n 1 i] 1 A_eq.append(row) b_eq.append(x[i]) # 期望产出约束: Y^T * lambda - sy y beta * gy for j in range(s): row np.zeros(n 1 m s k) row[:n] -Y[:, j] row[n] gy[j] row[n 1 m j] 1 A_eq.append(row) b_eq.append(-y[j]) # 非期望产出约束: B^T * lambda sb b - beta * gb for l in range(k): row np.zeros(n 1 m s k) row[:n] B[:, l] row[n] gb[l] row[n 1 m s l] 1 A_eq.append(row) b_eq.append(b[l]) # 变量边界: lambda 0, beta 0, 松弛量 0 bounds [(0, None)] * n [(0, None)] [(0, None)] * (m s k) res linprog(c, A_eqA_eq, b_eqb_eq, boundsbounds, methodhighs) if res.success: return res.x[n] # beta值 else: raise ValueError(线性规划求解失败检查数据是否存在异常值)这段代码的核心逻辑是通过构造一个凸锥前沿找到当前DMU在方向向量上能达到的最优改进程度。beta值越大说明该DMU离前沿越远效率越低。方向向量gx、gy、gb的选择直接影响结果——常见做法是取gxx、gyy、gbb即按当前DMU的实际投入产出设定方向这样beta就解释为可改进的比例。也有文献取gx1、gy1、gb1此时beta是绝对改进量。两种设定下GML指数的数值会不同但排序通常一致。我一般建议按实际值设定方向向量因为这样更符合“比例改进”的经济含义。2.3 全局前沿的构造与GML指数的分解计算全局前沿的构造很直接把所有时期的所有DMU放在一起形成一个大的生产可能集。假设有T个时期每个时期有n个DMU那么全局前沿就是T×n个DMU共同构成的前沿面。计算GML指数时需要分别计算四个方向距离函数值当期前沿下的t期和t1期值全局前沿下的t期和t1期值。GML指数的公式为GML (1 D^G(t, t1)) / (1 D^G(t, t))其中D^G表示全局前沿下的方向距离函数值。进一步分解为EC (1 D^t(t, t1)) / (1 D^t(t, t)) TC [(1 D^G(t, t1)) / (1 D^t(t, t1))] × [(1 D^t(t, t)) / (1 D^G(t, t))]EC衡量的是效率追赶效应TC衡量的是技术前沿移动效应。GML EC × TC。def gml_index(X_list, Y_list, B_list, gx, gy, gb): 计算GML指数及其分解 X_list, Y_list, B_list: 各时期的投入、期望产出、非期望产出矩阵列表 返回: GML, EC, TC 的时间序列 T len(X_list) n X_list[0].shape[0] # 构造全局前沿数据 X_global np.vstack(X_list) Y_global np.vstack(Y_list) B_global np.vstack(B_list) gml_list [] ec_list [] tc_list [] for t in range(T - 1): X_t, Y_t, B_t X_list[t], Y_list[t], B_list[t] X_t1, Y_t1, B_t1 X_list[t1], Y_list[t1], B_list[t1] gml_t [] ec_t [] tc_t [] for i in range(n): # 当期前沿下的t期和t1期 d_t_t sbm_ddf(X_t[i], Y_t[i], B_t[i], X_t, Y_t, B_t, gx, gy, gb) d_t_t1 sbm_ddf(X_t1[i], Y_t1[i], B_t1[i], X_t, Y_t, B_t, gx, gy, gb) # 全局前沿下的t期和t1期 d_g_t sbm_ddf(X_t[i], Y_t[i], B_t[i], X_global, Y_global, B_global, gx, gy, gb) d_g_t1 sbm_ddf(X_t1[i], Y_t1[i], B_t1[i], X_global, Y_global, B_global, gx, gy, gb) # GML指数 gml (1 d_g_t1) / (1 d_g_t) # 效率变化 ec (1 d_t_t1) / (1 d_t_t) # 技术变化 tc gml / ec gml_t.append(gml) ec_t.append(ec) tc_t.append(tc) gml_list.append(gml_t) ec_list.append(ec_t) tc_list.append(tc_t) return np.array(gml_list), np.array(ec_list), np.array(tc_list)这段代码的关键在于全局前沿下的方向距离函数值d_g_t和d_g_t1必须用同一套全局数据计算否则传递性不成立。另外方向向量gx、gy、gb在整个计算过程中必须保持一致不能在不同时期用不同的方向向量否则指数分解会失去意义。实际跑数据时我习惯先把所有时期的投入产出数据标准化到同一量纲避免因为单位差异导致线性规划数值不稳定。3. 数据准备与代码实现从原始面板到SBM-GML结果的全流程3.1 投入产出指标体系怎么搭才经得起审稿SBM-GML指数的结果可靠性七成取决于指标体系的设计。常见的绿色全要素生产率测算框架里投入指标一般包括劳动力从业人员数、资本存量永续盘存法计算、能源消费万吨标准煤。期望产出是GDP或工业增加值非期望产出是CO₂排放量、SO₂排放量、废水排放量等。这里有几个容易踩的坑资本存量的折旧率取多少很多文献取9.6%但如果你研究的是特定行业这个值可能需要调整。能源消费是取实物量还是标准煤建议统一折算成标准煤否则不同能源品种的加总没有意义。非期望产出的处理上CO₂排放量需要用排放系数法自己算不能直接用能源消费量代替。数据来源方面省级面板一般用《中国统计年鉴》《中国能源统计年鉴》和各省统计年鉴。地级市面板的数据可得性差一些可能需要从《中国城市统计年鉴》和各省市统计公报里手动整理。我一般会先做一个数据完整性检查如果某个城市某年的SO₂排放量缺失要么用插值法补要么直接剔除该样本不要用均值填充因为均值填充会人为降低该DMU的效率波动导致GML指数被低估。3.2 用Python跑通SBM-GML的最小可复现示例下面是一个完整的可运行示例用模拟数据演示从数据构造到GML指数计算的全过程。你可以直接把这段代码复制到Jupyter Notebook里跑替换成自己的数据即可。import numpy as np import pandas as pd # 1. 构造模拟数据 np.random.seed(42) n_dmu 30 # 30个决策单元 T 5 # 5个时期 X_list, Y_list, B_list [], [], [] for t in range(T): # 投入: 劳动力、资本、能源 labor np.random.uniform(100, 500, n_dmu) capital np.random.uniform(200, 800, n_dmu) energy np.random.uniform(50, 300, n_dmu) X np.column_stack([labor, capital, energy]) # 期望产出: GDP Y np.random.uniform(500, 2000, n_dmu).reshape(-1, 1) # 非期望产出: CO2、SO2 co2 energy * np.random.uniform(1.5, 2.5, n_dmu) so2 energy * np.random.uniform(0.01, 0.05, n_dmu) B np.column_stack([co2, so2]) X_list.append(X) Y_list.append(Y) B_list.append(B) # 2. 设定方向向量 # 取所有时期投入产出的均值作为方向向量 X_all np.vstack(X_list) Y_all np.vstack(Y_list) B_all np.vstack(B_list) gx X_all.mean(axis0) gy Y_all.mean(axis0) gb B_all.mean(axis0) # 3. 计算GML指数 gml, ec, tc gml_index(X_list, Y_list, B_list, gx, gy, gb) # 4. 整理结果 # gml形状: (T-1, n_dmu) # 计算每个时期的平均GML指数 for t in range(T-1): print(f时期 {t}-{t1}: GML均值{gml[t].mean():.4f}, fEC均值{ec[t].mean():.4f}, TC均值{tc[t].mean():.4f}) # 计算每个DMU的累积GML指数 cum_gml np.prod(gml, axis0) print(f\n累积GML指数范围: [{cum_gml.min():.4f}, {cum_gml.max():.4f}]) print(f累积GML指数均值: {cum_gml.mean():.4f})这段代码跑出来的结果GML均值应该在1附近波动。如果所有时期的GML均值都远大于1或远小于1说明方向向量设定有问题或者数据本身存在系统性偏差。正常情况下GML指数在1附近波动EC和TC的乘积严格等于GML。你可以用np.allclose(gml, ec * tc)验证一下如果返回True说明分解正确。参数说明gx、gy、gb的方向向量选择会影响beta的绝对值但不影响GML指数的排序。如果你想让结果更直观可以把方向向量设为所有DMU的均值这样beta就解释为“相对于平均水平的改进空间”。另外sbm_ddf函数里的methodhighs是scipy的最新线性规划求解器比旧的simplex方法更稳定建议保留。3.3 结果验证GML指数算出来之后怎么判断靠不靠谱算完GML指数第一件事是检查传递性。取任意三个连续时期t、t1、t2验证GML(t,t2)是否等于GML(t,t1)×GML(t1,t2)。如果不等说明全局前沿的构造有问题或者方向距离函数在跨期比较时出现了数值误差。第二件事是检查EC和TC的乘积是否等于GML这个在代码里已经保证了但如果你自己改了公式一定要重新验证。第三件事是看GML指数的分布如果大部分DMU的GML都大于1说明整体生产率在提升如果EC普遍小于1而TC大于1说明效率在退步但技术在进步这种“技术驱动型”增长在东部沿海地区比较常见。还有一个容易被忽略的验证步骤把GML指数的结果和传统ML指数对比。如果两者排序差异很大说明你的数据里存在明显的不可行解问题这时候GML指数的结果更可信。如果两者排序基本一致说明不可行解问题不严重但GML指数的传递性优势仍然存在。4. 避坑与排查SBM-GML指数计算中最容易翻车的五个地方4.1 线性规划无解现象、原因与解决现象跑linprog的时候返回res.success False或者beta值异常大比如大于10。原因通常有三个一是数据里有负值或零值SBM模型要求所有投入产出为正二是某个DMU在所有时期都是极端值导致前沿面被拉伸三是方向向量设得太小导致约束条件过于严格。解决办法先检查数据把所有零值替换为一个很小的正数比如1e-6把所有负值取绝对值或剔除该样本。如果某个DMU确实是极端值可以考虑用Winsorize缩尾处理把上下1%的值替换为1%分位数和99%分位数。4.2 指数分解不成立EC×TC≠GML的排查路径现象np.allclose(gml, ec * tc)返回False。原因通常是当期前沿和全局前沿的方向距离函数值计算时用了不同的方向向量或不同的数据矩阵。排查路径第一步检查sbm_ddf函数调用时传入的X、Y、B矩阵是否一致第二步检查方向向量gx、gy、gb是否在所有调用中保持不变第三步检查全局前沿数据是否包含了所有时期的DMU。如果这三步都没问题那可能是数值精度问题把np.allclose的容差从1e-8放宽到1e-6再试。4.3 结果全大于1或全小于1方向向量与量纲的隐藏陷阱现象所有DMU的GML指数都大于1.5或都小于0.5。原因方向向量设得太大或太小导致beta值被系统性放大或缩小。比如如果gx取的是所有DMU投入的最大值而某个DMU的投入只有最大值的十分之一那beta就会很小GML指数就会接近1。解决办法方向向量取所有DMU的均值或者取当前DMU的实际值。另外检查投入产出的量纲是否统一——如果劳动力单位是“人”资本单位是“亿元”能源单位是“万吨标准煤”那方向向量的三个分量差异会很大建议先做标准化处理。4.4 非期望产出处理不当CO₂排放算错导致GML虚高现象GML指数普遍偏高TC分量异常大。原因非期望产出的数据质量差或者排放系数用错了。比如用能源消费量直接代替CO₂排放量忽略了不同能源品种的排放因子差异。解决办法CO₂排放量 Σ(能源消费量 × 折标准煤系数 × 碳排放系数 × 氧化率)。原煤、焦炭、汽油、柴油的排放系数都不同不能统一用一个系数。另外非期望产出的方向向量gb建议取实际值不要取均值因为非期望产出的减少空间通常比期望产出的增加空间小。4.5 面板数据跨期比较全局前沿构造的三个常见错误现象GML指数在某个时期突然跳变或者EC和TC的走势完全相反。原因全局前沿构造时把不同时期的DMU混在一起但没有考虑技术异质性。比如2005年的DMU和2020年的DMU放在同一个前沿面上2020年的DMU自然更靠近前沿导致2005年的DMU效率被低估。解决办法如果研究时期跨度较大超过10年建议分阶段构造全局前沿或者用窗口DEA的方法每3-5年一个窗口。另外确保所有时期的DMU数量一致如果有城市在某个时期被合并或拆分需要做数据调整。5. 进阶技巧用超效率SBM-GML处理有效DMU的排序问题标准SBM模型有一个固有缺陷所有有效DMU的效率值都是1无法区分它们之间的优劣。如果你需要对这些有效DMU进行排序比如评选绿色发展的标杆城市标准SBM-GML就不够用了。超效率SBMSuper-SBM的思路是在计算某个DMU的效率时把它从参考集中剔除然后用剩余DMU构造前沿面。如果该DMU仍然有效它的效率值就会大于1从而实现排序。把超效率SBM和GML指数结合就是超效率SBM-GML指数。实现上只需要在sbm_ddf函数里加一个参数exclude_idx在构造X、Y、B矩阵时把当前DMU排除掉。但要注意超效率模型在全局前沿下可能会出现无可行解的情况尤其是当某个DMU在所有时期都是唯一有效的时候。这时候需要退回到标准SBM-GML或者用松弛量来辅助排序。def super_sbm_ddf(x, y, b, X, Y, B, gx, gy, gb, exclude_idx): 超效率SBM方向距离函数 exclude_idx: 需要排除的DMU索引 # 排除当前DMU mask np.ones(X.shape[0], dtypebool) mask[exclude_idx] False X_ex X[mask] Y_ex Y[mask] B_ex B[mask] # 调用标准SBM-DDF return sbm_ddf(x, y, b, X_ex, Y_ex, B_ex, gx, gy, gb)这个函数的逻辑很简单把当前DMU从参考集中剔除然后用剩余DMU计算beta。如果beta仍然为0即该DMU在剔除后仍然有效说明它是“超效率”的可以给它一个大于1的效率值。实际跑的时候我一般会先用标准SBM-GML算一遍找出所有有效DMU再用超效率SBM-GML对这些有效DMU重新排序。这样既保证了整体结果的稳定性又解决了有效DMU的区分问题。还有一个实用技巧如果你发现超效率SBM-GML的结果和标准SBM-GML的排序差异很大不要慌这通常说明你的数据里存在“一枝独秀”的DMU——它在所有时期都远离其他DMU。这时候建议检查一下这个DMU的数据是否真实或者考虑用Malmquist指数的Bootstrap方法做置信区间估计看看排序差异是否在统计上显著。最后说一个我自己的习惯每次跑完SBM-GML我都会把结果和原始数据放在一起做散点图看看GML指数高的DMU是不是真的在投入产出上有优势。有好几次我发现某个城市的GML指数异常高结果一查数据发现它的CO₂排放量少了一个数量级——原来是单位写错了。这种错误光看代码是发现不了的必须回到数据本身。希望帮到你。本文还有配套的精品资源点击获取