简介本资源面向经济学、管理学及区域发展研究领域的学习者与科研人员提供一套完整的全要素生产率测算工具重点解决SBM-GML指数、ML指数及超效率SBM的Matlab实现问题适合具备一定计量基础、希望快速上手效率评价与绿色生产率分析的用户。压缩包共16个文件包含11个m脚本、4份pdf说明文档和1份xlsx示例数据整体约6.17MB脚本覆盖VRS与CRS下含非期望产出的SBM及SBM-GML计算文档则从Matlab安装、代码调用到结果解读逐项配图讲解。目前已有5472人学习下载。资源不仅提供可直接运行的示例数据与完整代码还梳理了SBM-GML、GML-DDF与SBM-DDF三种主流方法的差异并附指标解读图形展示便于读者对照文献理解测算逻辑、复现结果并迁移至自身研究场景。1. SBM-GML指数到底在算什么从Malmquist到绿色全要素生产率的落地逻辑做效率测算的人迟早会撞上SBM-GML指数这个组合。它不是一个新指标而是把两套成熟方法焊在一起超效率SBM负责在每一期截面里给决策单元DMU算出一个能区分有效单元的效率值GMLGlobal Malmquist-Luenberger负责把多期效率值串成一条可分解的指数链。标题里提到的“本包计算的是SBM-GML指数”说的就是用超效率SBM做底层距离函数再用全局前沿构造ML指数最终得到绿色全要素生产率GTFP的变动及其分解项。这套东西解决的核心问题是当你的投入产出数据里包含非期望产出比如二氧化碳、废水、固废传统的Malmquist指数要么处理不了坏产出要么在线性规划无解时直接翻车。SBM-GML把非期望产出塞进方向距离函数的框架里同时用全局前沿避免“技术倒退”的假象。适合谁用做区域绿色发展评估、行业环境效率、企业碳绩效的硕博论文和课题报告尤其是需要把效率变动拆成“技术进步”和“效率追赶”两块的人。我见过太多人卡在第一步用MaxDEA或DEAP跑出一堆数却说不清GML指数大于1到底意味着什么。这篇就按“先立住原理、再动手复现、最后排坑”的顺序把SBM-GML从公式到代码到结果解读讲透。你不需要先成为运筹学专家但得愿意把每个参数的含义搞清楚——否则跑出来的结果自己都不敢信。2. 超效率SBM与GML的耦合原理为什么不能直接用传统Malmquist2.1 超效率SBM处理非期望产出的数学骨架传统SBM模型在效率值为1时无法区分有效DMU的优劣超效率SBM通过把被评价DMU从参考集中剔除让有效单元的效率值可以大于1。当引入非期望产出时目标函数变成最小化投入和期望产出的松弛、最大化非期望产出的松弛。用数学语言说假设有n个DMU每个DMU有m种投入、s种期望产出、h种非期望产出向量形式为x∈R^m、y^g∈R^s、y^b∈R^h。超效率SBM的效率值ρ*的规划式可以写成min ρ (1/m) * Σ( x̄_i / x_ik ) / [ (1/(sh)) * ( Σ( ȳ_r^g / y_rk^g ) Σ( ȳ_t^b / y_tk^b ) ) ] s.t. x̄_i ≥ Σ_{j≠k} λ_j x_ij ȳ_r^g ≤ Σ_{j≠k} λ_j y_rj^g ȳ_t^b ≥ Σ_{j≠k} λ_j y_tj^b λ_j ≥ 0, x̄_i ≥ x_ik, ȳ_r^g ≤ y_rk^g, ȳ_t^b ≥ y_tk^b这里的关键在于非期望产出的方向它被当作“越少越好”的产出所以在约束里是ȳ_t^b ≥ Σλ_j y_tj^b即实际非期望产出不能低于前沿面上的非期望产出。这个符号方向搞反结果会完全颠倒——这是血泪经验里最常见的一类错误。2.2 GML指数如何把多期效率串起来GML的核心是用全局生产可能性集PG替代当期生产可能性集Pt。全局前沿由所有时期的DMU共同构造这样技术退步就不会因为前沿面移动而出现虚假的指数下降。GML指数的定义是GML (1 D^G(x^t, y^t, b^t; g)) / (1 D^G(x^{t1}, y^{t1}, b^{t1}; g))其中D^G是全局方向距离函数。它可以分解为EC效率变化 当期技术效率的比值反映追赶效应TC技术进步 全局前沿与当期前沿的比值反映创新效应分解公式为 GML EC × TC。如果GML1说明绿色全要素生产率提升EC1说明落后单元在追赶前沿TC1说明整体技术前沿在向外推。2.3 为什么选SBM-GML而不是EBM-GML或方向距离函数常见做法里方向距离函数DDF也能处理非期望产出但它允许投入和产出同比例扩张容易高估效率。EBM是SBM和径向模型的混合参数调起来更玄学。SBM-GML的优势在于非径向、非角度直接处理松弛变量对数据量纲不敏感且超效率版本能给出全排序。代价是计算量随DMU数量增加而上升且线性规划求解器需要能处理非线性分式规划——通常通过Charnes-Cooper变换转成线性规划。我一般会建议如果DMU数量少于50、时期少于10年SBM-GML完全够用如果DMU超过200考虑用EBM或两阶段方法降维。这不是理论上的硬边界而是求解时间和数值稳定性的经验阈值。3. 用Python复现SBM-GML从数据准备到指数分解的完整链路3.1 数据格式与预处理投入、期望产出、非期望产出的三表结构数据准备是最容易埋雷的地方。你需要三组变量投入指标如劳动、资本、能源、期望产出如GDP、工业增加值、非期望产出如CO2、SO2、废水。每个DMU在每个时期都要有完整记录缺失值不能简单填零——填零会让该DMU在当期变成“超级有效”结果完全失真。我一般用长格式的DataFrame列名统一为dmu_id, year, input1, input2, ..., output_g1, output_g2, ..., output_b1, output_b2, ...。下面是一个模拟数据生成的代码import numpy as np import pandas as pd np.random.seed(42) dmus [fDMU{i} for i in range(1, 21)] years [2018, 2019, 2020, 2021, 2022] records [] for dmu in dmus: base_input np.random.uniform(50, 150, size2) base_good np.random.uniform(80, 200, size2) base_bad np.random.uniform(10, 60, size2) for t, year in enumerate(years): growth 1 0.03 * t np.random.normal(0, 0.02) records.append({ dmu_id: dmu, year: year, input1: base_input[0] * growth * np.random.normal(1, 0.05), input2: base_input[1] * growth * np.random.normal(1, 0.05), output_g1: base_good[0] * growth * np.random.normal(1, 0.08), output_g2: base_good[1] * growth * np.random.normal(1, 0.08), output_b1: base_bad[0] * (1 0.01 * t) * np.random.normal(1, 0.06), output_b2: base_bad[1] * (1 0.01 * t) * np.random.normal(1, 0.06), }) df pd.DataFrame(records) print(df.head()) print(df.shape)这段代码生成了20个DMU、5个时期、2投入2期望产出2非期望产出的面板数据。关键参数说明growth控制期望产出的增长趋势非期望产出用(10.01*t)模拟缓慢上升随机扰动项的标准差控制在0.05-0.08之间避免数据过于平滑导致所有DMU都在前沿面上。实际数据中你需要检查每个变量的最小值是否为正——SBM模型要求投入和产出均为正值出现零或负数需要先做平移处理。3.2 超效率SBM的线性规划求解用scipy搭建单期效率计算超效率SBM的分式规划可以通过Charnes-Cooper变换转为线性规划。下面用scipy.optimize.linprog实现单期效率计算。为了可读性我把投入产出矩阵整理成numpy数组并写一个函数返回每个DMU的效率值。from scipy.optimize import linprog def super_sbm_undesirable(inputs, good_outputs, bad_outputs): inputs: (n, m) 投入矩阵 good_outputs: (n, s) 期望产出矩阵 bad_outputs: (n, h) 非期望产出矩阵 返回: (n,) 效率值数组 n, m inputs.shape _, s good_outputs.shape _, h bad_outputs.shape efficiencies np.zeros(n) for k in range(n): # 决策变量: [lambda_1..lambda_n, x_bar_1..x_bar_m, y_g_bar_1..y_g_bar_s, y_b_bar_1..y_b_bar_h, t] # 总变量数 n m s h 1 num_vars n m s h 1 c np.zeros(num_vars) # 目标函数: 最小化 (1/m)Σ(x_bar/x_k) / t 的线性化形式 # 经Charnes-Cooper变换后目标为最小化 (1/m)Σ(x_bar_i / x_ik) for i in range(m): c[n i] 1.0 / (m * inputs[k, i]) # 注意: 这里省略了分母的归一化约束完整实现需要加入Σ(y_g_bar/y_gk)Σ(y_b_bar/y_bk) sh 的约束 # 为保持代码简洁以下用简化版演示约束结构 A_eq [] b_eq [] # 约束1: Σ_{j≠k} lambda_j * x_ij - x_bar_i 0 for i in range(m): row np.zeros(num_vars) for j in range(n): if j ! k: row[j] inputs[j, i] row[n i] -1 A_eq.append(row) b_eq.append(0) # 约束2: Σ_{j≠k} lambda_j * y_g_rj - y_g_bar_r 0 -Σ y_g_bar 0 for r in range(s): row np.zeros(num_vars) for j in range(n): if j ! k: row[j] -good_outputs[j, r] row[n m r] 1 A_eq.append(row) b_eq.append(0) # 约束3: Σ_{j≠k} lambda_j * y_b_tj - y_b_bar_t 0 for t in range(h): row np.zeros(num_vars) for j in range(n): if j ! k: row[j] bad_outputs[j, t] row[n m s t] -1 A_eq.append(row) b_eq.append(0) # 约束4: x_bar_i x_ik -x_bar_i -x_ik for i in range(m): row np.zeros(num_vars) row[n i] -1 A_eq.append(row) b_eq.append(-inputs[k, i]) # 约束5: y_g_bar_r y_g_rk for r in range(s): row np.zeros(num_vars) row[n m r] 1 A_eq.append(row) b_eq.append(good_outputs[k, r]) # 约束6: y_b_bar_t y_b_tk -y_b_bar_t -y_b_tk for t in range(h): row np.zeros(num_vars) row[n m s t] -1 A_eq.append(row) b_eq.append(-bad_outputs[k, t]) # 变量下界: lambda0, x_bar0, y_g_bar0, y_b_bar0, t0 bounds [(0, None)] * num_vars res linprog(c, A_ubA_eq, b_ubb_eq, boundsbounds, methodhighs) if res.success: efficiencies[k] res.fun else: efficiencies[k] np.nan return efficiencies这段代码展示了超效率SBM的核心约束结构但为了篇幅做了简化——完整的Charnes-Cooper变换需要加入分母归一化约束并且目标函数要包含非期望产出的松弛项。实际使用时我建议直接用pyDEA或pulp重写或者参考Tone(2001)的原始公式逐项实现。参数说明methodhighs是scipy目前最稳定的线性规划求解器bounds全部设为(0, None)保证非负如果res.success为False说明该DMU在当期无可行解通常是因为数据中存在极端值或零值。3.3 构造全局前沿并计算GML指数分解EC与TC的代码实现有了单期效率值还不够GML需要全局前沿下的距离函数。全局前沿的做法是把所有时期的投入产出数据堆叠在一起形成一个“超级参考集”然后对每个DMU在每个时期计算其相对于全局前沿的距离。下面代码演示如何组织全局数据并计算GML指数。def compute_gml(df, input_cols, good_cols, bad_cols): df: 长格式面板数据 返回: 包含GML、EC、TC的DataFrame years sorted(df[year].unique()) dmus sorted(df[dmu_id].unique()) results [] # 构造全局参考集 global_inputs df[input_cols].values global_good df[good_cols].values global_bad df[bad_cols].values for dmu in dmus: for t_idx in range(len(years) - 1): year_t years[t_idx] year_next years[t_idx 1] # 当期数据 row_t df[(df[dmu_id] dmu) (df[year] year_t)] row_next df[(df[dmu_id] dmu) (df[year] year_next)] if row_t.empty or row_next.empty: continue x_t row_t[input_cols].values yg_t row_t[good_cols].values yb_t row_t[bad_cols].values x_next row_next[input_cols].values yg_next row_next[good_cols].values yb_next row_next[bad_cols].values # 全局前沿下的距离函数简化用全局参考集计算效率值 # 实际应用中需要调用方向距离函数此处用超效率SBM值近似 eff_t_global super_sbm_undesirable( np.vstack([global_inputs, x_t]), np.vstack([global_good, yg_t]), np.vstack([global_bad, yb_t]) )[-1] eff_next_global super_sbm_undesirable( np.vstack([global_inputs, x_next]), np.vstack([global_good, yg_next]), np.vstack([global_bad, yb_next]) )[-1] # 当期前沿下的效率值 df_t df[df[year] year_t] df_next df[df[year] year_next] eff_t_current super_sbm_undesirable( np.vstack([df_t[input_cols].values, x_t]), np.vstack([df_t[good_cols].values, yg_t]), np.vstack([df_t[bad_cols].values, yb_t]) )[-1] eff_next_current super_sbm_undesirable( np.vstack([df_next[input_cols].values, x_next]), np.vstack([df_next[good_cols].values, yg_next]), np.vstack([df_next[bad_cols].values, yb_next]) )[-1] # GML 全局效率比值 gml eff_next_global / eff_t_global if eff_t_global ! 0 else np.nan # EC 当期效率比值 ec eff_next_current / eff_t_current if eff_t_current ! 0 else np.nan # TC GML / EC tc gml / ec if ec ! 0 else np.nan results.append({ dmu_id: dmu, year: f{year_t}-{year_next}, GML: gml, EC: ec, TC: tc }) return pd.DataFrame(results) # 调用示例 input_cols [input1, input2] good_cols [output_g1, output_g2] bad_cols [output_b1, output_b2] gml_results compute_gml(df, input_cols, good_cols, bad_cols) print(gml_results.head(10)) print(gml_results[[GML, EC, TC]].describe())这段代码的逻辑是对每个DMU的相邻两期分别计算其在全局前沿和当期前沿下的效率值然后做比值。参数说明eff_t_global和eff_next_global是全局前沿下的效率值eff_t_current和eff_next_current是当期前沿下的效率值。GML大于1表示绿色全要素生产率提升EC大于1表示追赶效应TC大于1表示技术进步。注意这里的super_sbm_undesirable函数需要传入完整的参考集实际计算时全局参考集包含所有DMU所有时期的数据计算量会显著增加。如果DMU数量超过100建议用并行计算或改用pulp的稀疏矩阵求解。4. 避坑与排查SBM-GML计算中最容易翻车的五个地方4.1 非期望产出的方向搞反导致效率值全部大于1现象跑出来的效率值大量超过1甚至出现3.0以上的极端值。原因非期望产出的约束方向写反了。在SBM模型中非期望产出是“越少越好”所以约束应该是实际非期望产出不能低于前沿面上的非期望产出即ȳ_b ≥ Σλ_j y_bj。如果写成≤模型会认为非期望产出越多越好导致所有DMU都“有效”。解决检查线性规划中非期望产出对应的约束行确保符号方向与理论一致。一个快速验证方法是把某个DMU的非期望产出翻倍如果效率值反而上升说明方向反了。4.2 全局前沿构造时忘记堆叠所有时期数据现象GML指数和EC指数完全相等TC恒等于1。原因全局前沿只用了当期数据没有把所有时期的数据堆叠进去。GML的核心就是全局参考集如果全局集等于当期集TC自然为1。解决在构造global_inputs、global_good、global_bad时确保用的是df的全量数据而不是df[df[year]year_t]。我一般会在代码里加一行断言assert len(global_inputs) len(df)防止自己手滑。4.3 数据中的零值和负值导致线性规划无解现象res.success返回False效率值为nan。原因SBM模型要求投入和产出均为正值数据中出现零或负数会导致约束不可行。常见于非期望产出中的某些年份排放为零或者投入指标做了标准化后出现负值。解决对零值做平移处理比如加一个很小的正数1e-6对负值先检查数据来源如果是标准化导致的换用min-max归一化到[0.1, 1]区间。注意平移会改变效率值的绝对大小但一般不影响排序和指数分解的方向。4.4 DMU数量少于投入产出指标总数导致大量DMU有效现象超过一半的DMU效率值为1超效率SBM也无法区分。原因DMU数量n小于投入数m加产出数s加非期望产出数h自由度不足几乎所有DMU都在前沿面上。这是DEA方法的固有局限。解决要么增加DMU数量比如把地级市换成区县要么减少指标数量用主成分分析降维要么改用随机前沿分析SFA。我一般会确保n ≥ 3×(msh)这是经验法则不是理论硬约束。4.5 指数分解结果与预期相反时先检查效率值的计算顺序现象GML大于1但TC小于1或者EC和TC的乘积不等于GML。原因效率值的计算顺序搞错了。GML eff_next_global / eff_t_globalEC eff_next_current / eff_t_currentTC GML / EC。如果先算TC再算EC或者用错了分子分母结果会完全乱套。解决在代码里加一行验证assert abs(GML - EC * TC) 1e-6如果不成立说明分解逻辑有误。这个断言帮我省了至少三天调试时间。5. 让结果经得起追问SBM-GML的稳健性检验与结果呈现技巧跑出GML指数只是第一步真正让审稿人或评审专家信服的是稳健性检验。我一般会做三件事第一换用EBM-GML重新计算看GML的符号和显著性是否一致第二把样本期拆成两段分别计算GML看技术进步项是否稳定第三对非期望产出做敏感性分析比如把CO2排放换成SO2排放看效率排序的Spearman相关系数是否高于0.8。如果这三项都通过结果基本可以写进论文。结果呈现上不要只放一张GML均值的折线图。我习惯用表格展示每个DMU的GML、EC、TC的均值和标准差再用核密度图展示GML的分布演变。下面是一个结果汇总的代码片段# 按DMU汇总GML、EC、TC的均值和标准差 summary gml_results.groupby(dmu_id).agg( GML_mean(GML, mean), GML_std(GML, std), EC_mean(EC, mean), TC_mean(TC, mean) ).reset_index() # 按年份汇总均值 yearly gml_results.groupby(year).agg( GML_mean(GML, mean), EC_mean(EC, mean), TC_mean(TC, mean) ).reset_index() print(summary.round(4)) print(yearly.round(4)) # 验证分解恒等式 gml_results[check] gml_results[EC] * gml_results[TC] - gml_results[GML] print(分解误差最大值:, gml_results[check].abs().max())这段代码做了两件事按DMU和按年份汇总指数并验证GML EC × TC的恒等式。参数说明groupby后的agg可以一次性算多个统计量check列用于验证分解误差正常情况下应该在1e-6以内。如果误差超过1e-4说明效率值计算中有数值不稳定需要检查线性规划的求解器设置。最后一个技巧如果审稿人问“为什么不用ML指数而用GML”标准回答是——ML指数在跨期比较时可能出现线性规划无解且无法处理非期望产出GML通过全局前沿解决了这两个问题同时保持了可分解性。这个回答我用了不下十次每次都能堵住追问。做效率测算这些年最大的教训就是不要迷信软件输出的结果每一个指数值背后都要能手工推一遍公式。希望帮到你。本文还有配套的精品资源点击获取