基于DEA的GTFP测算：ML与GML指数解析及Python实现

简介面向经济管理、资源环境等领域的科研人员与高校学生资源内提供一套基于MATLAB的GMLGlobal Malmquist-Luenberger指数与MLMalmquist-Luenberger指数测算代码适用于DEA框架下引入非期望产出如污染物排放的绿色全要素生产率GTFP分解可用于生产效率动态变化与可持续性评估。压缩包为RAR格式内仅包含1个m文件体积约1KB代码短小精悍聚焦DEA的GML分解核心算法可直接计算GML/ML指数及其分解项技术创新指数、效率变化指数帮助理解技术进步与效率变动的来源。目前已有2618人学习或下载得到较多使用者验证。通过运行该代码读者可快速掌握基于DEA的GML分解流程免去从零编写程序的成本结合代码逻辑还能进一步拓展到企业效率评价、行业生产率比较、环境政策效果分析等实际应用场景为实证研究与课程设计提供便利。1. GTFP测算为什么绕不开DEA、GML指数与ML指数这条主线做绿色全要素生产率GTFP测算DEA框架下的指数选型基本就两条路ML指数和GML指数。前者把跨期生产率放到当期技术前沿上做后者用全局参考集做DEA框架下的GML分解绕开前者跨期线性规划无解的硬伤。很多新手把GTFP测算当成「找数据、跑软件、出表」三步走真正上手才发现选错指数、跨期无解、分解结果对不上、技术进步方向解释反了这些问题远比跑一次DEA本身更难定位。这篇文章就按「指数在什么前提下成立、GML和ML怎么从公式走到代码、算完怎么验证」的顺序展开面向处理面板数据、想把非期望产出放进模型、又要交出可复现测算代码的从业者。2. ML指数的测算逻辑非期望产出如何进入DEA模型2.1 为什么传统Malmquist指数不能直接用于GTFP测算经典Malmquist指数的距离函数基于产出可自由处置的假设期望产出和非期望产出被同等对待。污染物作为坏产出进入模型时如果用普通DF运算一个通过末端治理减少二氧化硫排放的企业会被系统判定为“产出不足”因为减排行为让总产出径向变小却没有在模型中获得任何奖励。GTFP测算的起点就是修正这一点把非期望产出当作联合生产的副产物减排必须占用真实资源技术在“增产”和“减排”两个方向上同时发挥作用。Färe、Grosskopf、Linde和Pasurka在1989年提出的环境生产可能集是处理这个问题的最小公理框架。它包括三个关键约束一是产出弱可处置减少非期望产出需要放弃部分期望产出二是投入强可处置多余投入可以无偿清除三是联合生产不存在“只有好产出、没有坏产出”的可行生产点。熟悉DEA的人会注意到第三个假设意味着期望产出与非期望产出必须同时出现这直接决定了后续线性规划中约束条件的写法。2.2 定向技术距离函数把“增产”和“减排”放进同一个方向向量在环境生产可能集上测算GTFP方向性距离函数DDF是比谢泼德距离函数更自然的工具。给定方向向量 g(0, -b)一个决策单元 (x, y, b) 沿该方向扩张β倍投影点为 (x, (1β)y, (1-β)b)。这意味着决策单元可以用一个统一的β同时完成两个动作期望产出增加β非期望产出减少β。以t期参考集评价t期第k个决策单元线性规划写为max β s.t. Σ λj xj ≤ xk // 投入不增加 Σ λj yj ≥ (1β) yk // 期望产出扩张 Σ λj bj (1-β) bk // 非期望产出压缩 Σ λj 1 // 若采用VRS λj ≥ 0, β ≥ 0这里非期望产出约束用等式而非不等式是弱可处置性的直接体现。β的最优值衡量决策单元距离前沿的“绿色无效程度”β0说明该DMU已经在由样本构造的技术前沿上β越大则改进空间越大。理解这个规划是后续所有代码的基础ML和GML指数都只是把若干个这样的β值按不同方式组成可跨期比较的指数。2.2.1 投入处理上的两种常见方向向量方向向量如果写成 g(-x, y, -b)模型允许投入同步压缩β的含义就变成“投入产出同时改进的比例”这种写法更接近成本面如果写成 g(0, y, -b)则只评价产出侧效率。做GTFP的文献尤其Chung等1997年提出的原始ML指数默认产出侧方向向量。测算前先确认你采纳的方向否则后续分解中的技术变化项符号可能和文献对不上。2.3 ML指数的四点平均形式与EC/TC分解2.3.1 四种组合的DDF分别测什么ML指数是相邻两期、两个参考集、两个被评单元组合出的四个DDF值的几何平均。记 D_t(x_{t1}) 为用 t 期参考集评价 t1 期生产点的DDF则ML_t^{t1} sqrt{ [(1 D_t(x_t,y_t,b_t)) / (1 D_t(x_{t1},y_{t1},b_{t1}))] × [(1 D_{t1}(x_t,y_t,b_t)) / (1 D_{t1}(x_{t1},y_{t1},b_{t1}))] }四个分量中对角线组合 D_t(x_t) 和 D_{t1}(x_{t1}) 是当期效率交叉组合 D_t(x_{t1}) 和 D_{t1}(x_t) 是跨期效率。跨期组合回答的问题是如果t1期的生产方式被放在t期末的约束条件下离前沿多远这正是“技术进步”得以识别的原因——t1期的生产点若比t期前沿更远说明前沿没有跟上。ML可以分解为效率变化EC和技术变化TCEC (1 D_t(x_t)) / (1 D_{t1}(x_{t1})) TC sqrt{ [(1 D_{t1}(x_t)) / (1 D_t(x_t))] × [(1 D_{t1}(x_{t1})) / (1 D_t(x_{t1}))] }EC大于1说明决策单元在追赶当期前沿TC大于1说明前沿向外扩张。需要留意EC与TC的乘积在数值上严格等于ML这个恒等关系是后面代码自检的关键。2.3.2 跨期可行性两两搭配时的无解风险ML最被人诟病的问题在于交叉组合不保证可行。t期参考集里的投入产出组合可能根本无法容纳t1期的生产点尤其当样本包含技术跳跃、新工艺或产能突变时线性规划的标准形直接报告infeasible。这种现象在真实数据里非常频繁环保约束收紧后某些省市的污染物大幅下降产出结构也随产业结构调整变化两期生产点就可能落在彼此前沿覆盖范围之外。无解不是数据错误而是指数定义本身的缺陷。更麻烦的是即便有可行解ML指数不满足循环性连续两期的ML连乘不等于多年累计生产率变化因为每期参照前沿都在换。这两个缺陷叠加让ML只适合期数少、技术前沿相对平稳的短面板。如果你的面板横跨五年以上或者样本里存在明显的技术断代就需要考虑GML指数。3. GML指数及其GML分解全局参考集如何修复跨期错位3.1 全局生产可能集把1到T期观测并进一个参考集GML指数的核心改动只有一处把参考集从当期换成全局。Oh在2010年给出的定义里全局生产可能集是所有当期生产可能集的并集再取凸包即 P^G conv(P^1 ∪ P^2 ∪ … ∪ P^T)。落到代码层面就是把面板数据里全部时期的观测同时放进参考矩阵配合Σλ1的凸组合约束就自动完成了凸包操作。这个构造带来的第一个好处是数学层面的任何一个被评DMU无论它属于哪个时期本身都包含在全局参考集里。因此至少存在一个平凡可行解即λ取自身、β0线性规划恒有解。GML指数从根上消除了ML的跨期不可行问题这一点在面板期数多、决策单元技术路线分化大时尤其关键。3.2 GML公式与指数性质GML指数的表达式比ML简洁得多不需要交叉四项只计算同一个决策单元在前后两期、全局前沿下的两个β值GML_t^{t1} (1 D_G(x_t,y_t,b_t)) / (1 D_G(x_{t1},y_{t1},b_{t1}))分子是t期生产点距全局前沿的绿色无效程度加一分母是t1期对应值。比值大于1说明该DMU在全局前沿的意义上变得更有效率。这里有一个容易被忽略的点由于参考集里包含未来信息t期决策单元的D_G值通常大于它在当期前沿下的D值GML指数测的是“相对于整个样本期统一基准”的变化而不是“相对于当时技术环境”的变化。3.3 GML分解EC与BPC的拆解方式GML的分解没有用ML里的TC概念而是拆成效率变化EC和最佳实践差距变化BPC。BPCbest practice gap change度量的是目标时期的前沿与全局前沿之间的差距变化公式为EC (1 D_t(x_t)) / (1 D_{t1}(x_{t1})) BPC [(1 D_G(x_t)) / (1 D_t(x_t))] / [(1 D_G(x_{t1})) / (1 D_{t1}(x_{t1}))]EC与ML中的EC含义一致衡量追赶前沿的程度。BPC则比较两期技术水平各自与全局前沿的“距离”变化如果t1期前沿比t期前沿更接近全局前沿意味着发生了技术进步BPC大于1。乘积关系 GML EC × BPC 在构造上严格成立这也是测算代码里最直接的数值校验。3.3.1 BPC不等于前沿移动解释BPC时要谨慎。BPC衡量的是当期前沿与全局前沿之间差距的相对变化不是前沿绝对位置的移动量。如果两期前沿都在移动但幅度相同BPC可能接近1并不等于0。这意味着报告“技术进步率”时应当表述为“该期前沿逼近全局前沿的程度”而不是“技术前沿外移速度”。这个区别在评审较严格的论文里会被挑出来。3.4 ML与GML怎么选一张对照表维度ML指数GML指数参考集每期独立的当期前沿全部时期观测并集跨期不可行解经常出现恒有可行解循环性/累乘性不满足满足分解项EC TCEC BPC技术进步含义前沿外移当期前沿向全局前沿靠拢面板兼容性期数少、技术稳定长面板、存在技术跳跃计算量每个DMU需4次LP每个DMU需4次LP但无重算如果论文或项目对“技术变化”的定义要求严格对应传统Malmquist语义ML更容易解释如果样本期数超过五年或者中间有明显的政策冲击与结构调整GML是稳定得多的选择。实际项目里我一般两种情况都算把两个指数的结果并排报告差异越大越说明技术前沿发生了结构性变化这本身就是一个值得写的发现。4. GTFP测算代码Python实现GML与ML指数4.1 数据布局与列名约定测算代码的输入是一张长表每行代表一个决策单元在某时期的投入产出观测不要求面板完全平衡但每个DMU至少要有连续两期数据才能计算指数。列名约定对后面矩阵拼装影响很大我通常固定为dmu决策单元编号、period时期编号、x1、x2投入、y期望产出、b非期望产出。下面用一组模拟数据演示三期、八个决策单元x1和x2是投入y与x呈正相关b与y正相关并额外叠加部分低效扰动。随机种子固定为42保证任何人都能复现同一张结果表。import numpy as np import pandas as pd rng np.random.default_rng(42) def make_panel(T3, N8): rows [] for t in range(1, T 1): for i in range(1, N 1): x1 rng.uniform(1, 3) x2 rng.uniform(0.5, 2) # 每3个DMU中有一个完全有效其余存在绿色效率损失 slack 0.0 if i % 3 0 else rng.uniform(0.1, 0.5) y 2 * x1 x2 - slack rng.normal(0, 0.03) b 0.6 * y slack * 1.2 rng.uniform(0, 0.1) rows.append([i, t, x1, x2, y, b]) return pd.DataFrame(rows, columns[dmu, period, x1, x2, y, b]) df make_panel() df.head()模拟数据的逻辑说明slack0的DMU刻意落在前沿附近slack0的DMU同时表现为期望产出偏低和非期望产出偏高方向性距离函数会把这种低效识别为β0。如果你要换成真实数据只需把数据文件读成同结构的DataFrame不改变后续函数。4.2 SciPy求解定向距离函数的线性规划DDF的求解用SciPy的linprog。目标函数是最小化负β等价于最大化β变量向量是 [λ_1, …, λ_J, β]其中J是参考集中DMU的数量。投入约束用不等式期望产出约束转成小于等于号非期望产出约束用等式VRS时再追加一个等式Σλ1。from scipy.optimize import linprog def ddf(reference, target, vrsTrue): # reference: 参考集DataFrame # target: 被评DMU的一行含 x1,x2,y,b J len(reference) rb reference[[x1, x2, y, b]].values tg target[[x1, x2, y, b]].values.astype(float) c np.zeros(J 1) c[-1] -1.0 # 目标最大化 β A_ub, b_ub [], [] # 投入约束Σ λj xj xk逐个投入写入 for col_idx in [0, 1]: A_ub.append(list(rb[:, col_idx]) [0.0]) b_ub.append(tg[col_idx]) # 期望产出约束Σ λj yj - β yk yk # 转为标准形-Σ λj yj β yk -yk A_ub.append(list(-rb[:, 2]) [tg[2]]) b_ub.append(-tg[2]) A_eq, b_eq [], [] # 非期望产出弱可处置Σ λj bj β bk bk A_eq.append(list(rb[:, 3]) [tg[3]]) b_eq.append(tg[3]) if vrs: A_eq.append([1.0] * J [0.0]) # 凸组合约束 b_eq.append(1.0) bounds [(0, None)] * J [(0, None)] # λ≥0, β≥0 res linprog(c, A_ubnp.array(A_ub), b_ubnp.array(b_ub), A_eqnp.array(A_eq), b_eqnp.array(b_eq), boundsbounds, methodhighs) if not res.success: return np.nan return float(res.x[-1])代码说明A_ub按行拼装变量顺序固定为“参考集所有λ在前、β在最后”。目标是c[-1]-1因为linprog默认做最小化最大化β要转成最小化负β。β下界设为0禁止出现负的“效率改进”避免前沿外样本被错误解释为超高效。methodhighs是SciPy 1.6之后的推荐求解器对中等规模面板数据足够稳定。4.3 当期参考集与全局参考集的封装同一个ddf函数可以同时服务ML和GML区别只在传入的reference。当期参考集是从面板里筛出指定时期的所有行全局参考集是全部时期的整张表。封装成一个工厂函数避免在循环里反复写筛选逻辑。def get_reference(df, periodNone, use_globalFalse): if use_global: return df.copy() # 全局参考集 return df[df[period] period].copy() # 当期参考集这里不用use_global就等价于ML的参考集构造。需要注意全局参考集不要去掉重复DMU即使某个DMU在每期都出现它在不同时期的投入产出组合也被视为不同的生产观测全部保留才能体现“跨期技术并集”。4.4 循环计算所有DMU的GML、EC与BPC计算循环按相邻期对展开。对每个DMU需要四组DDF值当期前沿下的当期效率、全局前沿下的当期效率、下期前沿下的下期效率、全局前沿下的下期效率只有计算ML时才额外求两个跨期交叉项。def calc_indexes(df, vrsTrue): periods sorted(df[period].unique()) out [] for t, s in zip(periods[:-1], periods[1:]): ref_t get_reference(df, periodt) ref_s get_reference(df, periods) ref_g get_reference(df, use_globalTrue) for dmu in df[dmu].unique(): d_t df[(df[dmu] dmu) (df[period] t)].iloc[0] d_s df[(df[dmu] dmu) (df[period] s)].iloc[0] Dtt ddf(ref_t, d_t, vrs) Dts ddf(ref_t, d_s, vrs) # 跨期ML专用 Dst ddf(ref_s, d_t, vrs) # 跨期 Dss ddf(ref_s, d_s, vrs) Dgt ddf(ref_g, d_t, vrs) Dgs ddf(ref_g, d_s, vrs) row {dmu: dmu, period: f{t}-{s}} # GML 分解 if not any(np.isnan(x) for x in [Dtt, Dss, Dgt, Dgs]): ec (1 Dtt) / (1 Dss) bpc ((1 Dgt) / (1 Dtt)) / ((1 Dgs) / (1 Dss)) row[GML] ec * bpc row[EC] ec row[BPC] bpc # ML 指数 if not any(np.isnan(x) for x in [Dtt, Dss, Dts, Dst]): row[ML] np.sqrt(((1 Dtt) / (1 Dts)) * ((1 Dst) / (1 Dss))) out.append(row) return pd.DataFrame(out) res calc_indexes(df) print(res.round(4).to_string(indexFalse))这段代码里GML没有直接用(1Dgt)/(1Dgs)而是通过ec*bpc算出来目的就是让输出同时带三个列且天然满足GMLEC×BPC。ML列为NaN的位置正是第2章说的跨期不可行解——在当前这个模拟数据里你会看到少数DMU的ML缺失而GML列始终有值这是两个指数差异最直观的体现。4.4.1 非平衡面板的处理如果面板个别DMU某期缺失loc会抛KeyError。常见做法是先把数据按dmu分组只保留连续出现过目标期对的DMU再进循环。还有一种更稳妥的方式是把循环改成按分组迭代每组内部先对period排序再滑窗取相邻期这样天然跳过缺口。5. 结果解读与三个验证技巧让GML/ML测算经得起复核5.1 从环比指数到累计GTFP连乘的口径GML的循环性允许直接累积。以基期GTFP水平为1第t期的累计GTFP是路径上各期GML的连乘。Python里先按dmu分组对GML列做cumprod再把首期水平重标定为1即可。注意面板回归中通常使用累计水平值做被解释变量而做生产率增长的年度分解时使用环比GML更合适。两种口径对应的经济含义不同不能混用。5.2 结果准确性的一条硬校验测算脚本的出口处应该加一条断言abs(ec * bpc - gml) 1e-8。这不是形式主义。出现不等于1e-8的偏差几乎都是参考集切片错误或目标行选取错位比如把t期的Dss误传给了Dtt。另外检查当期有效DMU的β是否等于0若某DMU在t期达到当期前沿有效Dtt应当非常接近0此时它的EC项中分子为1。若β显著大于0却报告该DMU“效率为1”说明方向向量符号或弱可处置约束写反了。5.3 量纲归一化与小数值陷阱非期望产出的量纲如果远大于期望产出比如SO2排放量动辄数十万吨而GDP用亿元计量DDF的β会被小量纲变量主导导致结果几乎只反映一个产出的改进空间。处理办法是对所有投入产出列做(x - min) / (max - min) 1的归一化把数据平移到[1,2]区间。归一化不影响指数排序和分解方向但会让β值和效率分数更容易解释。测算报告里注明归一化方法也是评审常问的细节。最后说一个实操习惯把ddf函数单独存成模块测算脚本只负责数据预处理和循环逻辑。这样换数据时不会动到底层约束矩阵验证过一遍的求解器逻辑可以长期复用。指数算完后先用第5章的两条校验把脚本锁死再去做累计GTFP和后续计量分析会踏实很多。本文还有配套的精品资源点击获取