企业官网建设流程全解析

简介：面向经济管理、资源环境等领域的科研人员与高校学生，资源内提供一套基于MATLAB的GML（Global Malmquist-Luenberger）指数与ML（Malmquist-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_{t+1}) 为用 t 期参考集评价 t+1 期生产点的DDF，则：

ML_t^{t+1} = sqrt{ [(1 + D_t(x_t,y_t,b_t)) / (1 + D_t(x_{t+1},y_{t+1},b_{t+1}))] × [(1 + D_{t+1}(x_t,y_t,b_t)) / (1 + D_{t+1}(x_{t+1},y_{t+1},b_{t+1}))] }

四个分量中，对角线组合 D_t(x_t) 和 D_{t+1}(x_{t+1}) 是当期效率，交叉组合 D_t(x_{t+1}) 和 D_{t+1}(x_t) 是跨期效率。跨期组合回答的问题是：如果t+1期的生产方式被放在t期末的约束条件下，离前沿多远？这正是“技术进步”得以识别的原因——t+1期的生产点若比t期前沿更远，说明前沿没有跟上。

ML可以分解为效率变化EC和技术变化TC：

EC = (1 + D_t(x_t)) / (1 + D_{t+1}(x_{t+1})) TC = sqrt{ [(1 + D_{t+1}(x_t)) / (1 + D_t(x_t))] × [(1 + D_{t+1}(x_{t+1})) / (1 + D_t(x_{t+1}))] }

EC大于1说明决策单元在追赶当期前沿；TC大于1说明前沿向外扩张。需要留意，EC与TC的乘积在数值上严格等于ML，这个恒等关系是后面代码自检的关键。

2.3.2 跨期可行性：两两搭配时的无解风险

ML最被人诟病的问题在于交叉组合不保证可行。t期参考集里的投入产出组合，可能根本无法容纳t+1期的生产点，尤其当样本包含技术跳跃、新工艺或产能突变时，线性规划的标准形直接报告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^{t+1} = (1 + D_G(x_t,y_t,b_t)) / (1 + D_G(x_{t+1},y_{t+1},b_{t+1}))

分子是t期生产点距全局前沿的绿色无效程度加一，分母是t+1期对应值。比值大于1，说明该DMU在全局前沿的意义上变得更有效率。这里有一个容易被忽略的点：由于参考集里包含未来信息，t期决策单元的D_G值通常大于它在当期前沿下的D值，GML指数测的是“相对于整个样本期统一基准”的变化，而不是“相对于当时技术环境”的变化。

3.3 GML分解：EC与BPC的拆解方式

GML的分解没有用ML里的TC概念，而是拆成效率变化EC和最佳实践差距变化BPC。BPC（best practice gap change）度量的是目标时期的前沿与全局前沿之间的差距变化，公式为：

EC = (1 + D_t(x_t)) / (1 + D_{t+1}(x_{t+1})) BPC = [(1 + D_G(x_t)) / (1 + D_t(x_t))] / [(1 + D_G(x_{t+1})) / (1 + D_{t+1}(x_{t+1}))]

EC与ML中的EC含义一致，衡量追赶前沿的程度。BPC则比较两期技术水平各自与全局前沿的“距离”变化：如果t+1期前沿比t期前沿更接近全局前沿，意味着发生了技术进步，BPC大于1。乘积关系 GML = EC × BPC 在构造上严格成立，这也是测算代码里最直接的数值校验。

3.3.1 BPC不等于前沿移动

解释BPC时要谨慎。BPC衡量的是当期前沿与全局前沿之间差距的相对变化，不是前沿绝对位置的移动量。如果两期前沿都在移动但幅度相同，BPC可能接近1，并不等于0。这意味着报告“技术进步率”时，应当表述为“该期前沿逼近全局前沿的程度”，而不是“技术前沿外移速度”。这个区别在评审较严格的论文里会被挑出来。

3.4 ML与GML怎么选：一张对照表

如果论文或项目对“技术变化”的定义要求严格对应传统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(T=3, N=8): 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()

模拟数据的逻辑说明：slack=0的DMU刻意落在前沿附近，slack>0的DMU同时表现为期望产出偏低和非期望产出偏高，方向性距离函数会把这种低效识别为β>0。如果你要换成真实数据，只需把数据文件读成同结构的DataFrame，不改变后续函数。

4.2 SciPy求解定向距离函数的线性规划

DDF的求解用SciPy的linprog。目标函数是最小化负β，等价于最大化β；变量向量是 [λ_1, …, λ_J, β]，其中J是参考集中DMU的数量。投入约束用不等式，期望产出约束转成小于等于号，非期望产出约束用等式，VRS时再追加一个等式Σλ=1。

from scipy.optimize import linprog def ddf(reference, target, vrs=True): # 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_ub=np.array(A_ub), b_ub=np.array(b_ub), A_eq=np.array(A_eq), b_eq=np.array(b_eq), bounds=bounds, method='highs') if not res.success: return np.nan return float(res.x[-1])

代码说明：A_ub按行拼装，变量顺序固定为“参考集所有λ在前、β在最后”。目标是c[-1]=-1，因为linprog默认做最小化，最大化β要转成最小化负β。β下界设为0，禁止出现负的“效率改进”，避免前沿外样本被错误解释为超高效。method='highs'是SciPy 1.6之后的推荐求解器，对中等规模面板数据足够稳定。

4.3 当期参考集与全局参考集的封装

同一个ddf函数可以同时服务ML和GML，区别只在传入的reference。当期参考集是从面板里筛出指定时期的所有行；全局参考集是全部时期的整张表。封装成一个工厂函数，避免在循环里反复写筛选逻辑。

def get_reference(df, period=None, use_global=False): 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, vrs=True): periods = sorted(df['period'].unique()) out = [] for t, s in zip(periods[:-1], periods[1:]): ref_t = get_reference(df, period=t) ref_s = get_reference(df, period=s) ref_g = get_reference(df, use_global=True) 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(index=False))

这段代码里GML没有直接用(1+Dgt)/(1+Dgs)，而是通过ec*bpc算出来，目的就是让输出同时带三个列，且天然满足GML=EC×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和后续计量分析，会踏实很多。

本文还有配套的精品资源，点击获取