企业官网建设流程全解析

简介：本资源面向需要测算全要素生产率的科研人员与研究生，提供基于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-Luenberger（ML）指数算出来的结果，要么线性规划无解，要么在跨期方向距离函数下出现技术倒退的假象，要么效率值大于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到t+1的指数乘以t+1到t+2的指数，不等于t到t+2的指数。这两个问题在面板数据较长时尤其明显，很多论文里出现的效率值异常波动，根源就在这里。

GML指数（Global Malmquist-Lenberger）的核心改进是构造一个全局前沿（global frontier），即所有时期的决策单元（DMU）共同构成一个生产可能集。这样一来，所有DMU都在同一个前沿面上比较，既避免了不可行解，又天然满足传递性。Oh（2010）证明了GML指数可以分解为效率变化（EC）和技术变化（TC），且EC和TC的乘积严格等于GML指数。这个性质在实证分析中非常重要，因为你可以清楚地看到效率提升和技术进步各自的贡献。

2.2 SBM方向距离函数的设定与方向向量选择

SBM（Slack-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/(m+s+k)) * (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_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs') if res.success: return res.x[n] # beta值 else: raise ValueError("线性规划求解失败，检查数据是否存在异常值")

这段代码的核心逻辑是：通过构造一个凸锥前沿，找到当前DMU在方向向量上能达到的最优改进程度。beta值越大，说明该DMU离前沿越远，效率越低。方向向量gx、gy、gb的选择直接影响结果——常见做法是取gx=x、gy=y、gb=b，即按当前DMU的实际投入产出设定方向，这样beta就解释为可改进的比例。也有文献取gx=1、gy=1、gb=1，此时beta是绝对改进量。两种设定下GML指数的数值会不同，但排序通常一致。我一般建议按实际值设定方向向量，因为这样更符合“比例改进”的经济含义。

2.3 全局前沿的构造与GML指数的分解计算

全局前沿的构造很直接：把所有时期的所有DMU放在一起，形成一个大的生产可能集。假设有T个时期，每个时期有n个DMU，那么全局前沿就是T×n个DMU共同构成的前沿面。计算GML指数时，需要分别计算四个方向距离函数值：当期前沿下的t期和t+1期值，全局前沿下的t期和t+1期值。

GML指数的公式为：

GML = (1 + D^G(t, t+1)) / (1 + D^G(t, t))

其中D^G表示全局前沿下的方向距离函数值。进一步分解为：

EC = (1 + D^t(t, t+1)) / (1 + D^t(t, t)) TC = [(1 + D^G(t, t+1)) / (1 + D^t(t, t+1))] × [(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[t+1], Y_list[t+1], B_list[t+1] gml_t = [] ec_t = [] tc_t = [] for i in range(n): # 当期前沿下的t期和t+1期 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期和t+1期 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(axis=0) gy = Y_all.mean(axis=0) gb = B_all.mean(axis=0) # ========== 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}->{t+1}: GML均值={gml[t].mean():.4f}, " f"EC均值={ec[t].mean():.4f}, TC均值={tc[t].mean():.4f}") # 计算每个DMU的累积GML指数 cum_gml = np.prod(gml, axis=0) 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函数里的method='highs'是scipy的最新线性规划求解器，比旧的simplex方法更稳定，建议保留。

3.3 结果验证：GML指数算出来之后怎么判断靠不靠谱

算完GML指数，第一件事是检查传递性。取任意三个连续时期t、t+1、t+2，验证GML(t,t+2)是否等于GML(t,t+1)×GML(t+1,t+2)。如果不等，说明全局前沿的构造有问题，或者方向距离函数在跨期比较时出现了数值误差。第二件事是检查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就不够用了。超效率SBM（Super-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], dtype=bool) 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₂排放量少了一个数量级——原来是单位写错了。这种错误，光看代码是发现不了的，必须回到数据本身。希望帮到你。

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