尧图网页设计 
简介针对GML指数与ML指数测算需求这份MATLAB源码包提供了完整的GTFP测算与DEA分解实现适用于环境经济效率分析、生产效率评价等场景适合经济学研究者、数据分析人员及需要量化技术进步与效率变化的从业者。压缩包共1个文件为.m格式源码文件大小约1KB代码涵盖数据预处理、DEA模型构建、效率计算、Malmquist指数计算及GML分解等核心环节结构紧凑便于直接调用或二次开发。已有2619人学习下载。通过运行该代码可清晰理解GML指数如何引入非期望产出ML指数如何分解为技术创新与效率变化并获得从原始输入输出数据到最终结果输出的完整测算链路结合描述中的CCR/BCC模型与分解逻辑还能帮助读者深入掌握DEA环境绩效评价的基本原理为可持续性评估、企业效率对标及政策模拟提供可复用的计算工具。1. GML指数和ML指数测算GTFP为什么现在论文都在追GML新手拿到GML指数和ML指数测算GTFP的代码时第一反应是翻几篇文献照着跑。但真正上手会发现ML指数虽然老牌却在传递性和线性规划可行两个地方反复“翻车”GML指数用全局前沿把这两个坑一次填平几乎是现阶段做绿色全要素生产率测算最顺滑的选择。这篇文章把GML和ML的差异、DEA的GML分解怎么落地、代码怎么写、数据里哪些地方容易埋雷讲清楚。正在做学术论文或政策评估的读者手里只要有一份面板数据就可以按下面的流程算出GML指数、EC和TC分解结果再累积成GTFP做后续回归。2. 方向性距离函数与生产前沿ML指数为什么输给GML指数2.1 传统DEA装不下污染变量方向性距离函数才是GTFP的地基传统DEA的CCR和BCC模型只处理“多投入多产出”并且所有产出默认是好产出。GTFP测算里GDP要增加、CO2要减少这两种方向相反的产出放进径向模型时会出现一个荒谬的结果GDP和CO2同时被径向压缩完全没有体现“绿色”的含义。方向性距离函数解决了这个问题。方向性距离函数DDF的直观意思给定一个方向向量 g (g_y, -g_b)让被评价单元往期望产出增加、非期望产出减少的方向投影投影能走多远β就是多远。形式是D(x, y, b; g_y, -g_b) sup { β : (y βg_y, b - βg_b) ∈ P(x) }这里P(x)是当前投入水平下所有可行产出组合的集合。方向向量最常用的取法是用被评价单元自身的期望产出和非期望产出来定义比如 g_y y_0g_b b_0。这个取法让β天然代表“还能同时涨多少产出、降多少污染”。在计算上DDF需要解一个线性规划。投入一般取不等号约束允许不把投入用完非期望产出必须取等号约束这对应文献里的弱可处置性——想减少污染总得牺牲一部分资源配置效率。方向性距离函数这一步做对了后面的GML和ML才有意义。2.2 ML指数的几何平均构造跨期混合距离函数带来的两个病根ML指数是Chung、Färe和Grosskopf在上世纪九十年代提出的全称是Malmquist-Luenberger指数。它用一个几何平均值连接相邻两期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^{t1},y^{t1},b^{t1}) 是用第t期的生产前沿去评价第t1期的生产活动。如果第t1期的产出组合和投入组合落在第t期前沿包络的范围之外这个线性规划就没有可行解。在技术快速变化、产业结构剧变的样本里这类无解根本不是小概率事件。第二几何平均破坏了传递性。ML^{1,3} 并不等于 ML^{1,2} 乘以 ML^{2,3}这意味着不同基期算出来的累计生产率指数会互相矛盾。后面要拿GTFP做回归时传递性一旦失守整个面板数据的可比性就崩了。2.3 GML指数的全局前沿填平LP无解与传递性两个坑GML指数Global Malmquist-Luenberger Index由Oh提出核心改动只有一处把所有时期的生产可能集放在一起取它们的凸包构成一个全局前沿 P^G(x)。然后所有距离函数都基于这个全局前沿计算GML^{t,t1} (1 D^G(x^t,y^t,b^t)) / (1 D^G(x^{t1},y^{t1},b^{t1}))这个改动的价值不能用“技巧”来形容。全局前沿包含了历史所有时期的观测被评价单位无论在哪个时期都能在全局前沿的产出集合里找到参照因此线性规划天然可行——ML最让人头疼的“本次LP无解”在GML里直接消失。传递性也随之恢复。因为GML是同一前沿下的比值相邻两期的GML乘起来中间的项正好约分所以 GML^{1,T} ∏ GML^{t,t1}。这意味着你可以放心地把各期GML连乘得到累计GTFP拿它做门槛回归、面板回归或者政策效应评估都不会出现前后不一致。GML还能做和ML形式十分接近的分解分成EC效率变化和TC技术变化见第4章。经济学含义上EC衡量的是被评价单元向其当期生产前沿靠近的速度TC衡量的是全局前沿相对当期前沿的移动程度。这个分解比ML的几何平均交叉项干净得多结果也更可靠。2.4 选型清单什么时候GML什么时候ML还能用很多新手会问能不能偷懒用ML我的经验是分场景只有3年以内的短面板、技术变动又不明显时ML和GML数值差距很小ML算起来简单一点但任何超过3年、跨机构或跨区域的面板或者要拿GTFP做回归的场景我会直接上GML。审稿人对ML的“不可传递性”和“线性规划不可行”是写在条件反射里的你不想在返修时解释一堆。下面这张表可以直接放进论文或笔记里。维度ML指数GML指数生产前沿每年单独构造一次全部时期共同构造全局前沿跨期混合LP可能无解必然可行传递性不成立成立累计GTFP与基期选择有关与基期无关分解结构几何平均交叉项单一比值EC与TC直接相乘近年文献采用度少数老文章主流和高引文章标题里写“DEA的GML分解”GML本身不是一种新模型而是把DEA算出的方向性距离函数放进一个指数框架。图方便时很多人直接用MaxDEA一键输出但要想搞明白结果是不是对的还是得看懂这一章的距离函数。3. GML指数测算代码实现数据结构、LP求解与参数设置3.1 数据准备最少七列的面板结构GML测算的输入不复杂但最常出问题的就是数据表没排对。一份能直接跑的面板数据至少要有这些列区域省份/企业/行业ID、年份、劳动投入L、资本存量K、能源投入E可选但建议加、期望产出Y实际GDP、非期望产出CCO2或SO2等。数据必须是平衡面板吗不是非平衡面板可以跑但累计GTFP在缺口处要断开重新累积。资本存量千万别直接用当年固定资产投资额——DEA要求的是“存量”概念一般用永续盘存法PIM折算基年选择越早越好折旧率建议在全文中写死并说明。GDP和投资需要统一平减到某一年不变价不然跨期比较全是幻觉。非期望产出的计量单位建议与文献一致。比如CO2通常用万吨不同文献的单位可能不一样。量纲不会改变GML的相对排序因为DEA约束是逐维度写入的但极端数量级会让线性规划求解精度变差。遇到“上一期很明显是有效的怎么这期变成无解”这种玄学先查查是否C列有0值。3.2 方向性距离函数的线性规划变量、目标与约束对每个DMU i在第t期我们要算两类距离函数全局前沿距离 D^G 和当期前沿距离 D^t。核心是同一个LP只是参考集不同。以全局前沿为例设参考集中有J个观测全部时期的DMU启用的变量是β和J个λ权重目标函数 max β约束条件投入约束sum_j λ_j X_{jk} ≤ x_{0k}对所有投入k期望产出约束sum_j λ_j Y_j ≥ (1 β) y_0非期望产出约束sum_j λ_j C_j (1 - β) c_0λ_j ≥ 0β ≥ 0投入约束是不等号期望产出是不等号非期望产出必须是等号。原因前面已经说过污染排放是弱可处置的产出与污染之间做不到完全剥离。当参考集换成某年当期DMU集合时就是当期前沿距离。方向向量在这里取的是被评价DMU自身的产出和污染值 (g_y, g_b) (y_0, c_0)。如果想让β的含义变成“绝对增减量”可以把方向向量换成单位向量(1,1)LP的右端项也要相应改。绝大多数GTFP文献都采用自我方向建议保持主流做法因为β天然就是百分比变化报告结果时不用额外解释。3.3 MATLAB核心代码全局方向性距离函数求解把上面的LP搬进MATLAB要转化成linprog的标准形式。决策变量顺序取 [β; λ_1; ...; λ_J]目标函数 f [-1; zeros(J,1)]因为linprog默认求最小。下面这个函数只做一件事给定参考集和当前被评DMU返回方向性距离β。function [beta, lambda] ddf_solve(X, Y, C, x0, y0, c0, crs) % ddf_solve: 计算方向性距离函数 % 输入: % X: 参考集投入矩阵, n行m列 % Y: 参考集期望产出, n行1列 % C: 参考集非期望产出, n行1列 % x0: 被评价DMU投入, 1行m列 % y0: 被评价DMU期望产出, 标量 % c0: 被评价DMU非期望产出, 标量 % crs: 1-CRS, 0-VRS, 建议默认CRS % 输出: % beta: 方向性距离值(目标函数最优解) % lambda: 参考集内各DMU的权重 n size(X, 1); m size(X, 2); % 决策变量: [beta; lambda(1:n)] f [-1; zeros(n, 1)]; Aineq []; bineq []; % 投入可自由处置: sum(lambda_j * X(j,k)) x0(k) for k 1:m Aineq [Aineq; 0, X(:, k)]; bineq [bineq; x0(k)]; end % 期望产出扩张: sum(lambda_j * Y_j) (1beta) * y0 % 标准形式: -y0 * beta - sum(lambda_j * Y_j) -y0 Aineq [Aineq; -y0, -Y]; bineq [bineq; -y0]; % 非期望产出弱可处置: sum(lambda_j * C_j) (1-beta) * c0 % 标准形式: c0 * beta sum(lambda_j * C_j) c0 Aeq [c0, C]; beq c0; % VRS时追加 sum(lambda) 1 if crs 0 Aeq [Aeq; 0, ones(1, n)]; beq [beq; 1]; end lb [0; zeros(n, 1)]; % beta和lambda都非负 ub []; opts optimoptions(linprog, Display, off, Algorithm, dual-simplex); [x, ~, flag] linprog(f, Aineq, bineq, Aeq, beq, lb, ub, opts); if flag 0 beta x(1); lambda x(2:end); else beta NaN; lambda NaN(n, 1); end end逻辑说明函数不区分全局前沿还是当期前沿区别只在调用者传入的参考集。crs参数控制CRS还是VRSCRS下λ没有任何和的约束VRS下加一个 sum(lambda)1。代码里flag是linprog的求解状态遇到无解返回NaN——在GML里如果还有NaN八成是数据问题。参数说明方向向量取被评估DMU自己的y0和c0所以LP右端项里有y0和c0。若想改成g(1,1)把期望产出约束右端改成 y0 beta非期望产出约束右端改成 c0 - beta 即可但β的经济含义会不同不要混用。3.4 组装GML、EC、TC遍历面板的完整主程序有了ddf_solve还差最后一步对每个DMU每一年分别调一次全局前沿和当期前沿。主程序如下。%% 主程序: 遍历N个DMU、T年, 计算GML/EC/TC % id: N*T行1列, 区域或企业编号 % year: N*T行1列, 年份 % K, L, E, Y, C: N*T行1列, 已整理好 N length(unique(id)); T length(unique(year)); D_G zeros(N, T); D_t zeros(N, T); for i 1:N for t 1:T ridx find(id i year t); % 被评价观测的行号 x0 [K(ridx), L(ridx), E(ridx)]; y0 Y(ridx); c0 C(ridx); % 全局前沿: 参考集是全体年份 X_all [K, L, E]; D_G(i, t) ddf_solve(X_all, Y, C, x0, y0, c0, 1); % 当期前沿: 参考集只取第t年的观测 t_idx find(year t); D_t(i, t) ddf_solve(X_all(t_idx, :), Y(t_idx), C(t_idx), ... x0, y0, c0, 1); end end %% 由距离函数组装指数 EC zeros(N, T-1); TC zeros(N, T-1); GML zeros(N, T-1); for i 1:N for t 1:T-1 EC(i, t) (1 D_t(i, t)) / (1 D_t(i, t1)); TC(i, t) ((1 D_G(i, t)) / (1 D_t(i, t))) * ... ((1 D_t(i, t1)) / (1 D_G(i, t1))); GML(i, t) (1 D_G(i, t)) / (1 D_G(i, t1)); end end逻辑说明主程序先算两类D值再按第2章的公式组装三个指数。GML直接等于EC乘TC写GML只是为了输出便利。D_t和D_G的维度都是N×T注意D_t是“当期前沿测当期”不要当成“测跨期”——GML不需要跨期混合距离函数这正是它稳健的地方。参数说明这里crs固定传1CRS。如果你想试VRS把两处ddf_solve的最后一个参数改成0。几乎所有的GTFP文献都默认CRS代际可比的结果大多建立在CRS上。E这一列如果缺失就把x0里去掉EX_all也去掉E列剩下K和L两个投入也完全没问题。4. DEA的GML分解与累计GTFP从环比指数到论文汇报口径4.1 GML分解公式EC和TC各自衡量什么在DEA框架下GML一分为二常见写法是EC (1 D^t_t) / (1 D^{t1}_{t1})TC (1 D^G_t) / (1 D^t_t) × (1 D^{t1}{t1}) / (1 D^G{t1})其中 D^t_t 表示第t期的当期前沿测第t期DMUD^G_t 表示全局前沿测第t期DMU。第一个式子EC诀窍是当DMU向当期前沿靠近时D_t会减小EC就会大于1。第二个式子TC本质上衡量的是“全局前沿和当期前沿之间的距离变化”如果被评DMU处在技术进步很快的行业中它所在的当期前沿快速向全局前沿靠拢TC大于1。注意GML的TC和ML的TC不是同一个东西。ML的TC是两期前沿直接比高低GML的TC是当期前沿与全局前沿的差距变化。所以论文里不要直接拿GML的TC和别人paper里的TC数值做比较量级可能差0.3以上。如果你在复现文献时发现GML的TC整体小于ML原因也在这里不是代码错了。4.2 累计GTFP的定基换算为什么GML能直接连乘GML是环比指数两个相邻年份的GML相乘就等于跨期的累计生产率变化因为全局前沿是同一个。把基期GTFP设为100那么DMU i在第t年的累计GTFP为CumGTFP_{i,t} 100 × ∏_{τ2}^{t} GML_{i,τ-1,τ}写成MATLAB脚本非常简单%% 从GML累计GTFP, 基期100 CumGTFP zeros(N, T); CumGTFP(:, 1) 100; for i 1:N for t 2:T CumGTFP(i, t) CumGTFP(i, t-1) * GML(i, t-1); end end提示非平衡面板的每段连续年份要分开累计缺口处以100重置基期再把各段结果拼接成完整面板。否则累计过程会把缺年份当成GML1算进去低估增长。逻辑说明这个循环等价于 cumprod每一年的累计值等于上一年累计值乘当年的环比GML。基期取100是文献里的通用做法好处是回归系数可以直接解读成“相对基期变化百分之多少”。4.3 分解结果汇报论文常用的四张表有了EC、TC和CumGTFP一封论文结果通常汇报四块内容各年平均GML及其分解分DMU类型的累计GTFP均值对比EC与TC的散点图回归模型里GTFP变量的描述性统计。这里只提格式要求数学上不要漏掉基期说明表格里最好标注“GMLEC×TC”和“基期GTFP100”。在做分组均值对比时注意EC和TC是几何平均后的小数直接求算术平均会高估组内均值最好先求组内几何均值再跨组对比。很多论文的表格里EC和TC的均值与GML均值对不上就是因为把算术平均和几何平均混用了。4.4 GML与ML分解结果对不上的深层原因常见困惑为什么我用同一份数据算出来的ML和GML差异很大原因不在于公式而在于前沿的构造。ML每年重新构造一次前沿且跨期LP存在不可行GML全局前沿包络了所有时期所以它捕捉到的技术进步速度会平和一些。在技术快速进步的行业GML的TC通常小于ML的TC这是一种统计口径差异不是算错了。如果论文需要做稳健性不建议只报ML或只报GML可以把两者都列出来用一句话说明“GML在全局前沿下满足传递性故以GML为主”其余交给附表。这个过程不仅能堵审稿人的嘴也能帮自己发现数据里的异常点——当GML和ML在某一年突然分歧极大时那一年大概率有数据录入错误。5. GML指数测算避坑指南5个高频翻车点与排错方法5.1 现象LP无解beta返回NaN原因当期前沿包络不住极端观测。如果用了ML的跨期混合距离函数技术断层会让第t期前沿无法评价第t1期的产出组合GML的全局前沿包含所有观测原则上不会无解。如果GML下ddf_solve返回NaN通常是参考集里某列全是0Aeq退化成了全零行。解决切换到GML后仍然无解就查数据——某年的Y或C是否有0值某投入列是否存在空值被读成NaN找到问题行后剔除或做缺失填补不要继续带NaN跑循环。5.2 现象β异常大甚至超过10原因方向向量取值有0或极度接近0。如果某DMU某年的CO2排放非常接近0比如新型清洁企业c0≈0会让非期望产出等式约束形同虚设β被拉得异常大。解释上完全不排碳的企业反而被判定“效率极高”这是方向性距离函数的黑匣子效应。解决确保非期望产出为正且量纲合理确实存在0排放观测时把它剔除或者改用单位方向向量 g(1,1) 做非径向DDF并在论文里写清楚β的含义变化。别硬算数据里有一个0值结果就全歪。5.3 现象累计GTFP曲线出现背离常识的下滑原因非平衡面板累计断点处理错了。上面第4章提示过如果在CumGTFP循环里直接乘上缺失年份的GML1相当于把断档期当成了技术进步率为1累计指数会被拉低一格。面板里中间缺两年后面的累计值全部偏低10%以上。解决对每个DMU的连续年份段分段累计断点处重新设基期为100最后用merge拼回完整面板。审稿人最喜欢抓这种细节处理干净能少很多返修。5.4 现象资本存量数据忽高忽低GML结果随折旧率假设大幅波动原因永续盘存法参数不统一。资本K在GTFP里是存量各省/企业口径决定结果。永续盘存法要确定基年资本存量和折旧率δ文献里常取9.6%、10.96%、7%不等。把δ从7%改成11%累计GTFP的年均增速可能差0.5个百分点这个波动足以改变回归显著性。解决主回归用9.6%保持主流口径同时附表给δ7%和δ12%的敏感性分析。别不写参数DEA不帮你处理经济口径问题。折旧率、基期、平减方式这三个参数建议直接写进代码注释。5.5 现象CRS和VRS结果差异巨大甚至符号反转原因样本里规模效率占主导。CRS假设规模报酬不变VRS放开了规模效率。在GTFP测算里大多数文献用CRS因为GML的全局前沿在CRS下数学性质更干净且与方向性距离函数的弱可处置性兼容。如果样本里同时有小微型DMU和大型DMUVRS会让EC里混入规模效率变化GML分解结果和CRS完全不同。解决论文正文以CRS为主稳健性里提一句“改用VRS后结论方向不变”即可。如果你发现VRS下符号都反了先别急着改结论检查是不是样本里存在产出规模极端分化的DMU考虑对样本做缩尾或分组处理。6. 验证GML测算结果的三个技巧复现闭环与检验要点6.1 用成熟软件做交叉验证主程序跑完别急着信。拿自己数据的一个子集比如3个省份、5年丢进MaxDEA或类似工具对比GML和EC/TC的数值。注意MaxDEA的参数区要选对方向性距离函数和全局前沿选项。如果MATLAB和软件结果绝对值一致、相对误差在0.01以内说明LP构造没有问题不一致时优先查非期望产出的等号约束有没有被编译错。6.2 检查方向性距离函数分布区间GML算出的D_G绝大多数情况下应落在0到1之间。把D_G导出后看一眼均值和中位数如果某年某DMU的D_G大于1说明数据里有异常观测回到5.2排查。这个检查比看GML指数本身更早暴露问题——GML是比值结构两边同时出问题时错误会被掩盖D_G是原始层问题藏不住。6.3 从文献里“抄”一组结果做复现测试我每次写新代码版本时都会找一个已发表论文的公开数据表哪怕只有10行按对方的投入产出口径重跑一遍。重点不是抄数据而是强迫自己处理对方表格里那些“看起来不规整”的细节比如折旧率、平减基期、非期望产出的单位。能还原对方的GML均值说明你的代码在别人的数据上也成立还原不出来基本可以断定是自己的口径设定有偏差。这几年我改得最多的不是求解器而是资本存量和价格基期。先把口径写进代码注释再开始跑数据这是最省时间的一条路。希望帮到你。本文还有配套的精品资源点击获取 GML指数测算GTFP全攻略：从方向性距离函数到DEA代码实现与避坑
简介针对GML指数与ML指数测算需求这份MATLAB源码包提供了完整的GTFP测算与DEA分解实现适用于环境经济效率分析、生产效率评价等场景适合经济学研究者、数据分析人员及需要量化技术进步与效率变化的从业者。压缩包共1个文件为.m格式源码文件大小约1KB代码涵盖数据预处理、DEA模型构建、效率计算、Malmquist指数计算及GML分解等核心环节结构紧凑便于直接调用或二次开发。已有2619人学习下载。通过运行该代码可清晰理解GML指数如何引入非期望产出ML指数如何分解为技术创新与效率变化并获得从原始输入输出数据到最终结果输出的完整测算链路结合描述中的CCR/BCC模型与分解逻辑还能帮助读者深入掌握DEA环境绩效评价的基本原理为可持续性评估、企业效率对标及政策模拟提供可复用的计算工具。1. GML指数和ML指数测算GTFP为什么现在论文都在追GML新手拿到GML指数和ML指数测算GTFP的代码时第一反应是翻几篇文献照着跑。但真正上手会发现ML指数虽然老牌却在传递性和线性规划可行两个地方反复“翻车”GML指数用全局前沿把这两个坑一次填平几乎是现阶段做绿色全要素生产率测算最顺滑的选择。这篇文章把GML和ML的差异、DEA的GML分解怎么落地、代码怎么写、数据里哪些地方容易埋雷讲清楚。正在做学术论文或政策评估的读者手里只要有一份面板数据就可以按下面的流程算出GML指数、EC和TC分解结果再累积成GTFP做后续回归。2. 方向性距离函数与生产前沿ML指数为什么输给GML指数2.1 传统DEA装不下污染变量方向性距离函数才是GTFP的地基传统DEA的CCR和BCC模型只处理“多投入多产出”并且所有产出默认是好产出。GTFP测算里GDP要增加、CO2要减少这两种方向相反的产出放进径向模型时会出现一个荒谬的结果GDP和CO2同时被径向压缩完全没有体现“绿色”的含义。方向性距离函数解决了这个问题。方向性距离函数DDF的直观意思给定一个方向向量 g (g_y, -g_b)让被评价单元往期望产出增加、非期望产出减少的方向投影投影能走多远β就是多远。形式是D(x, y, b; g_y, -g_b) sup { β : (y βg_y, b - βg_b) ∈ P(x) }这里P(x)是当前投入水平下所有可行产出组合的集合。方向向量最常用的取法是用被评价单元自身的期望产出和非期望产出来定义比如 g_y y_0g_b b_0。这个取法让β天然代表“还能同时涨多少产出、降多少污染”。在计算上DDF需要解一个线性规划。投入一般取不等号约束允许不把投入用完非期望产出必须取等号约束这对应文献里的弱可处置性——想减少污染总得牺牲一部分资源配置效率。方向性距离函数这一步做对了后面的GML和ML才有意义。2.2 ML指数的几何平均构造跨期混合距离函数带来的两个病根ML指数是Chung、Färe和Grosskopf在上世纪九十年代提出的全称是Malmquist-Luenberger指数。它用一个几何平均值连接相邻两期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^{t1},y^{t1},b^{t1}) 是用第t期的生产前沿去评价第t1期的生产活动。如果第t1期的产出组合和投入组合落在第t期前沿包络的范围之外这个线性规划就没有可行解。在技术快速变化、产业结构剧变的样本里这类无解根本不是小概率事件。第二几何平均破坏了传递性。ML^{1,3} 并不等于 ML^{1,2} 乘以 ML^{2,3}这意味着不同基期算出来的累计生产率指数会互相矛盾。后面要拿GTFP做回归时传递性一旦失守整个面板数据的可比性就崩了。2.3 GML指数的全局前沿填平LP无解与传递性两个坑GML指数Global Malmquist-Luenberger Index由Oh提出核心改动只有一处把所有时期的生产可能集放在一起取它们的凸包构成一个全局前沿 P^G(x)。然后所有距离函数都基于这个全局前沿计算GML^{t,t1} (1 D^G(x^t,y^t,b^t)) / (1 D^G(x^{t1},y^{t1},b^{t1}))这个改动的价值不能用“技巧”来形容。全局前沿包含了历史所有时期的观测被评价单位无论在哪个时期都能在全局前沿的产出集合里找到参照因此线性规划天然可行——ML最让人头疼的“本次LP无解”在GML里直接消失。传递性也随之恢复。因为GML是同一前沿下的比值相邻两期的GML乘起来中间的项正好约分所以 GML^{1,T} ∏ GML^{t,t1}。这意味着你可以放心地把各期GML连乘得到累计GTFP拿它做门槛回归、面板回归或者政策效应评估都不会出现前后不一致。GML还能做和ML形式十分接近的分解分成EC效率变化和TC技术变化见第4章。经济学含义上EC衡量的是被评价单元向其当期生产前沿靠近的速度TC衡量的是全局前沿相对当期前沿的移动程度。这个分解比ML的几何平均交叉项干净得多结果也更可靠。2.4 选型清单什么时候GML什么时候ML还能用很多新手会问能不能偷懒用ML我的经验是分场景只有3年以内的短面板、技术变动又不明显时ML和GML数值差距很小ML算起来简单一点但任何超过3年、跨机构或跨区域的面板或者要拿GTFP做回归的场景我会直接上GML。审稿人对ML的“不可传递性”和“线性规划不可行”是写在条件反射里的你不想在返修时解释一堆。下面这张表可以直接放进论文或笔记里。维度ML指数GML指数生产前沿每年单独构造一次全部时期共同构造全局前沿跨期混合LP可能无解必然可行传递性不成立成立累计GTFP与基期选择有关与基期无关分解结构几何平均交叉项单一比值EC与TC直接相乘近年文献采用度少数老文章主流和高引文章标题里写“DEA的GML分解”GML本身不是一种新模型而是把DEA算出的方向性距离函数放进一个指数框架。图方便时很多人直接用MaxDEA一键输出但要想搞明白结果是不是对的还是得看懂这一章的距离函数。3. GML指数测算代码实现数据结构、LP求解与参数设置3.1 数据准备最少七列的面板结构GML测算的输入不复杂但最常出问题的就是数据表没排对。一份能直接跑的面板数据至少要有这些列区域省份/企业/行业ID、年份、劳动投入L、资本存量K、能源投入E可选但建议加、期望产出Y实际GDP、非期望产出CCO2或SO2等。数据必须是平衡面板吗不是非平衡面板可以跑但累计GTFP在缺口处要断开重新累积。资本存量千万别直接用当年固定资产投资额——DEA要求的是“存量”概念一般用永续盘存法PIM折算基年选择越早越好折旧率建议在全文中写死并说明。GDP和投资需要统一平减到某一年不变价不然跨期比较全是幻觉。非期望产出的计量单位建议与文献一致。比如CO2通常用万吨不同文献的单位可能不一样。量纲不会改变GML的相对排序因为DEA约束是逐维度写入的但极端数量级会让线性规划求解精度变差。遇到“上一期很明显是有效的怎么这期变成无解”这种玄学先查查是否C列有0值。3.2 方向性距离函数的线性规划变量、目标与约束对每个DMU i在第t期我们要算两类距离函数全局前沿距离 D^G 和当期前沿距离 D^t。核心是同一个LP只是参考集不同。以全局前沿为例设参考集中有J个观测全部时期的DMU启用的变量是β和J个λ权重目标函数 max β约束条件投入约束sum_j λ_j X_{jk} ≤ x_{0k}对所有投入k期望产出约束sum_j λ_j Y_j ≥ (1 β) y_0非期望产出约束sum_j λ_j C_j (1 - β) c_0λ_j ≥ 0β ≥ 0投入约束是不等号期望产出是不等号非期望产出必须是等号。原因前面已经说过污染排放是弱可处置的产出与污染之间做不到完全剥离。当参考集换成某年当期DMU集合时就是当期前沿距离。方向向量在这里取的是被评价DMU自身的产出和污染值 (g_y, g_b) (y_0, c_0)。如果想让β的含义变成“绝对增减量”可以把方向向量换成单位向量(1,1)LP的右端项也要相应改。绝大多数GTFP文献都采用自我方向建议保持主流做法因为β天然就是百分比变化报告结果时不用额外解释。3.3 MATLAB核心代码全局方向性距离函数求解把上面的LP搬进MATLAB要转化成linprog的标准形式。决策变量顺序取 [β; λ_1; ...; λ_J]目标函数 f [-1; zeros(J,1)]因为linprog默认求最小。下面这个函数只做一件事给定参考集和当前被评DMU返回方向性距离β。function [beta, lambda] ddf_solve(X, Y, C, x0, y0, c0, crs) % ddf_solve: 计算方向性距离函数 % 输入: % X: 参考集投入矩阵, n行m列 % Y: 参考集期望产出, n行1列 % C: 参考集非期望产出, n行1列 % x0: 被评价DMU投入, 1行m列 % y0: 被评价DMU期望产出, 标量 % c0: 被评价DMU非期望产出, 标量 % crs: 1-CRS, 0-VRS, 建议默认CRS % 输出: % beta: 方向性距离值(目标函数最优解) % lambda: 参考集内各DMU的权重 n size(X, 1); m size(X, 2); % 决策变量: [beta; lambda(1:n)] f [-1; zeros(n, 1)]; Aineq []; bineq []; % 投入可自由处置: sum(lambda_j * X(j,k)) x0(k) for k 1:m Aineq [Aineq; 0, X(:, k)]; bineq [bineq; x0(k)]; end % 期望产出扩张: sum(lambda_j * Y_j) (1beta) * y0 % 标准形式: -y0 * beta - sum(lambda_j * Y_j) -y0 Aineq [Aineq; -y0, -Y]; bineq [bineq; -y0]; % 非期望产出弱可处置: sum(lambda_j * C_j) (1-beta) * c0 % 标准形式: c0 * beta sum(lambda_j * C_j) c0 Aeq [c0, C]; beq c0; % VRS时追加 sum(lambda) 1 if crs 0 Aeq [Aeq; 0, ones(1, n)]; beq [beq; 1]; end lb [0; zeros(n, 1)]; % beta和lambda都非负 ub []; opts optimoptions(linprog, Display, off, Algorithm, dual-simplex); [x, ~, flag] linprog(f, Aineq, bineq, Aeq, beq, lb, ub, opts); if flag 0 beta x(1); lambda x(2:end); else beta NaN; lambda NaN(n, 1); end end逻辑说明函数不区分全局前沿还是当期前沿区别只在调用者传入的参考集。crs参数控制CRS还是VRSCRS下λ没有任何和的约束VRS下加一个 sum(lambda)1。代码里flag是linprog的求解状态遇到无解返回NaN——在GML里如果还有NaN八成是数据问题。参数说明方向向量取被评估DMU自己的y0和c0所以LP右端项里有y0和c0。若想改成g(1,1)把期望产出约束右端改成 y0 beta非期望产出约束右端改成 c0 - beta 即可但β的经济含义会不同不要混用。3.4 组装GML、EC、TC遍历面板的完整主程序有了ddf_solve还差最后一步对每个DMU每一年分别调一次全局前沿和当期前沿。主程序如下。%% 主程序: 遍历N个DMU、T年, 计算GML/EC/TC % id: N*T行1列, 区域或企业编号 % year: N*T行1列, 年份 % K, L, E, Y, C: N*T行1列, 已整理好 N length(unique(id)); T length(unique(year)); D_G zeros(N, T); D_t zeros(N, T); for i 1:N for t 1:T ridx find(id i year t); % 被评价观测的行号 x0 [K(ridx), L(ridx), E(ridx)]; y0 Y(ridx); c0 C(ridx); % 全局前沿: 参考集是全体年份 X_all [K, L, E]; D_G(i, t) ddf_solve(X_all, Y, C, x0, y0, c0, 1); % 当期前沿: 参考集只取第t年的观测 t_idx find(year t); D_t(i, t) ddf_solve(X_all(t_idx, :), Y(t_idx), C(t_idx), ... x0, y0, c0, 1); end end %% 由距离函数组装指数 EC zeros(N, T-1); TC zeros(N, T-1); GML zeros(N, T-1); for i 1:N for t 1:T-1 EC(i, t) (1 D_t(i, t)) / (1 D_t(i, t1)); TC(i, t) ((1 D_G(i, t)) / (1 D_t(i, t))) * ... ((1 D_t(i, t1)) / (1 D_G(i, t1))); GML(i, t) (1 D_G(i, t)) / (1 D_G(i, t1)); end end逻辑说明主程序先算两类D值再按第2章的公式组装三个指数。GML直接等于EC乘TC写GML只是为了输出便利。D_t和D_G的维度都是N×T注意D_t是“当期前沿测当期”不要当成“测跨期”——GML不需要跨期混合距离函数这正是它稳健的地方。参数说明这里crs固定传1CRS。如果你想试VRS把两处ddf_solve的最后一个参数改成0。几乎所有的GTFP文献都默认CRS代际可比的结果大多建立在CRS上。E这一列如果缺失就把x0里去掉EX_all也去掉E列剩下K和L两个投入也完全没问题。4. DEA的GML分解与累计GTFP从环比指数到论文汇报口径4.1 GML分解公式EC和TC各自衡量什么在DEA框架下GML一分为二常见写法是EC (1 D^t_t) / (1 D^{t1}_{t1})TC (1 D^G_t) / (1 D^t_t) × (1 D^{t1}{t1}) / (1 D^G{t1})其中 D^t_t 表示第t期的当期前沿测第t期DMUD^G_t 表示全局前沿测第t期DMU。第一个式子EC诀窍是当DMU向当期前沿靠近时D_t会减小EC就会大于1。第二个式子TC本质上衡量的是“全局前沿和当期前沿之间的距离变化”如果被评DMU处在技术进步很快的行业中它所在的当期前沿快速向全局前沿靠拢TC大于1。注意GML的TC和ML的TC不是同一个东西。ML的TC是两期前沿直接比高低GML的TC是当期前沿与全局前沿的差距变化。所以论文里不要直接拿GML的TC和别人paper里的TC数值做比较量级可能差0.3以上。如果你在复现文献时发现GML的TC整体小于ML原因也在这里不是代码错了。4.2 累计GTFP的定基换算为什么GML能直接连乘GML是环比指数两个相邻年份的GML相乘就等于跨期的累计生产率变化因为全局前沿是同一个。把基期GTFP设为100那么DMU i在第t年的累计GTFP为CumGTFP_{i,t} 100 × ∏_{τ2}^{t} GML_{i,τ-1,τ}写成MATLAB脚本非常简单%% 从GML累计GTFP, 基期100 CumGTFP zeros(N, T); CumGTFP(:, 1) 100; for i 1:N for t 2:T CumGTFP(i, t) CumGTFP(i, t-1) * GML(i, t-1); end end提示非平衡面板的每段连续年份要分开累计缺口处以100重置基期再把各段结果拼接成完整面板。否则累计过程会把缺年份当成GML1算进去低估增长。逻辑说明这个循环等价于 cumprod每一年的累计值等于上一年累计值乘当年的环比GML。基期取100是文献里的通用做法好处是回归系数可以直接解读成“相对基期变化百分之多少”。4.3 分解结果汇报论文常用的四张表有了EC、TC和CumGTFP一封论文结果通常汇报四块内容各年平均GML及其分解分DMU类型的累计GTFP均值对比EC与TC的散点图回归模型里GTFP变量的描述性统计。这里只提格式要求数学上不要漏掉基期说明表格里最好标注“GMLEC×TC”和“基期GTFP100”。在做分组均值对比时注意EC和TC是几何平均后的小数直接求算术平均会高估组内均值最好先求组内几何均值再跨组对比。很多论文的表格里EC和TC的均值与GML均值对不上就是因为把算术平均和几何平均混用了。4.4 GML与ML分解结果对不上的深层原因常见困惑为什么我用同一份数据算出来的ML和GML差异很大原因不在于公式而在于前沿的构造。ML每年重新构造一次前沿且跨期LP存在不可行GML全局前沿包络了所有时期所以它捕捉到的技术进步速度会平和一些。在技术快速进步的行业GML的TC通常小于ML的TC这是一种统计口径差异不是算错了。如果论文需要做稳健性不建议只报ML或只报GML可以把两者都列出来用一句话说明“GML在全局前沿下满足传递性故以GML为主”其余交给附表。这个过程不仅能堵审稿人的嘴也能帮自己发现数据里的异常点——当GML和ML在某一年突然分歧极大时那一年大概率有数据录入错误。5. GML指数测算避坑指南5个高频翻车点与排错方法5.1 现象LP无解beta返回NaN原因当期前沿包络不住极端观测。如果用了ML的跨期混合距离函数技术断层会让第t期前沿无法评价第t1期的产出组合GML的全局前沿包含所有观测原则上不会无解。如果GML下ddf_solve返回NaN通常是参考集里某列全是0Aeq退化成了全零行。解决切换到GML后仍然无解就查数据——某年的Y或C是否有0值某投入列是否存在空值被读成NaN找到问题行后剔除或做缺失填补不要继续带NaN跑循环。5.2 现象β异常大甚至超过10原因方向向量取值有0或极度接近0。如果某DMU某年的CO2排放非常接近0比如新型清洁企业c0≈0会让非期望产出等式约束形同虚设β被拉得异常大。解释上完全不排碳的企业反而被判定“效率极高”这是方向性距离函数的黑匣子效应。解决确保非期望产出为正且量纲合理确实存在0排放观测时把它剔除或者改用单位方向向量 g(1,1) 做非径向DDF并在论文里写清楚β的含义变化。别硬算数据里有一个0值结果就全歪。5.3 现象累计GTFP曲线出现背离常识的下滑原因非平衡面板累计断点处理错了。上面第4章提示过如果在CumGTFP循环里直接乘上缺失年份的GML1相当于把断档期当成了技术进步率为1累计指数会被拉低一格。面板里中间缺两年后面的累计值全部偏低10%以上。解决对每个DMU的连续年份段分段累计断点处重新设基期为100最后用merge拼回完整面板。审稿人最喜欢抓这种细节处理干净能少很多返修。5.4 现象资本存量数据忽高忽低GML结果随折旧率假设大幅波动原因永续盘存法参数不统一。资本K在GTFP里是存量各省/企业口径决定结果。永续盘存法要确定基年资本存量和折旧率δ文献里常取9.6%、10.96%、7%不等。把δ从7%改成11%累计GTFP的年均增速可能差0.5个百分点这个波动足以改变回归显著性。解决主回归用9.6%保持主流口径同时附表给δ7%和δ12%的敏感性分析。别不写参数DEA不帮你处理经济口径问题。折旧率、基期、平减方式这三个参数建议直接写进代码注释。5.5 现象CRS和VRS结果差异巨大甚至符号反转原因样本里规模效率占主导。CRS假设规模报酬不变VRS放开了规模效率。在GTFP测算里大多数文献用CRS因为GML的全局前沿在CRS下数学性质更干净且与方向性距离函数的弱可处置性兼容。如果样本里同时有小微型DMU和大型DMUVRS会让EC里混入规模效率变化GML分解结果和CRS完全不同。解决论文正文以CRS为主稳健性里提一句“改用VRS后结论方向不变”即可。如果你发现VRS下符号都反了先别急着改结论检查是不是样本里存在产出规模极端分化的DMU考虑对样本做缩尾或分组处理。6. 验证GML测算结果的三个技巧复现闭环与检验要点6.1 用成熟软件做交叉验证主程序跑完别急着信。拿自己数据的一个子集比如3个省份、5年丢进MaxDEA或类似工具对比GML和EC/TC的数值。注意MaxDEA的参数区要选对方向性距离函数和全局前沿选项。如果MATLAB和软件结果绝对值一致、相对误差在0.01以内说明LP构造没有问题不一致时优先查非期望产出的等号约束有没有被编译错。6.2 检查方向性距离函数分布区间GML算出的D_G绝大多数情况下应落在0到1之间。把D_G导出后看一眼均值和中位数如果某年某DMU的D_G大于1说明数据里有异常观测回到5.2排查。这个检查比看GML指数本身更早暴露问题——GML是比值结构两边同时出问题时错误会被掩盖D_G是原始层问题藏不住。6.3 从文献里“抄”一组结果做复现测试我每次写新代码版本时都会找一个已发表论文的公开数据表哪怕只有10行按对方的投入产出口径重跑一遍。重点不是抄数据而是强迫自己处理对方表格里那些“看起来不规整”的细节比如折旧率、平减基期、非期望产出的单位。能还原对方的GML均值说明你的代码在别人的数据上也成立还原不出来基本可以断定是自己的口径设定有偏差。这几年我改得最多的不是求解器而是资本存量和价格基期。先把口径写进代码注释再开始跑数据这是最省时间的一条路。希望帮到你。本文还有配套的精品资源点击获取