用Python的PuLP库搞定NDDF模型：一个环境经济学研究生的避坑实战笔记
环境经济学研究生的NDDF模型实战：用Python和PuLP破解绿色全要素生产率测算难题
第一次在顶刊论文里看到"基于全局生产技术的弱可处置性非径向方向距离函数"这个术语时，我的大脑瞬间宕机。作为环境经济学研二学生，导师要求我复现一篇关于中国区域绿色全要素生产率的论文，而文中的NDDF模型就像天书——公式里满是β、λ、g这些神秘符号，参考文献动辄指向1997年的经典文献。更绝望的是，作者根本没提供代码。经过两个月挣扎，我终于从理论到代码完整走通了NDDF建模全流程，这份笔记记录了我掉过的坑和爬出来的经验，特别适合被复杂模型困扰的硕博生参考。
1. 理解NDDF：从公式恐惧到概念拆解
1.1 方向距离函数家族图谱
第一次接触NDDF前，需要理清几个关键概念的关系：
- DEA：数据包络分析，效率评价的经典方法
- DDF：方向距离函数，考虑期望产出增加与非期望产出减少
- NDDF：非径向方向距离函数，允许投入产出项按不同比例调整
关键区别：传统DDF要求所有变量同比例变化，而NDDF通过权重向量ω实现差异化调整，更符合现实经济系统。
1.2 弱可处置性：环境约束的核心体现
这个概念困扰我最久。简单说：
- 强可处置性：可以无成本处理非期望产出（如随意排放污染物）
- 弱可处置性：减少污染需要牺牲部分资源（如安装净化设备会占用生产资金）
在代码中体现为对非期望产出的约束条件写法不同：
if self.disp == "weak disposability": prob += pulp.lpSum([self.weights[j0] * self.bad_outs.values[j0][s] for j0 in self._j]) == self.bad_outs.values[j0][s] - self.betab[s] * self.gb.values[j0][s] elif self.disp == "strong disposability": prob += pulp.lpSum([self.weights[j0] * self.bad_outs.values[j0][s] for j0 in self._j]) >= self.bad_outs.values[j0][s] - self.betab[s] * self.gb.values[j0][s]1.3 全局生产技术 vs 当期生产技术
在构建生产前沿面时：
| 技术类型 | 参考集构成 | 适用场景 |
|---|---|---|
| 全局生产技术 | 所有时期DMU的并集 | 跨期效率比较 |
| 当期生产技术 | 仅同期DMU | 截面分析 |
选择全局技术时，代码中约束条件的求和范围需要覆盖所有时期：
pulp.lpSum([(self.weights[j0] * self.inputs.values[j0][i]) for j0 in self._j]) # 遍历所有DMU不论时期2. PuLP建模关键：方向向量与约束条件
2.1 方向向量g的陷阱
原始论文给出g=(0,0,-E,Y,-D,-S,-W)，但在代码实现时发现矛盾：
# 错误理解会导致约束条件符号混乱 sum(λ*E) ≤ E - β*(-E) # 实际变为 sum(λ*E) ≤ E + βE # 正确写法应保持约束条件统一符号 sum(λ*E) ≤ E + βE # 此时g中E分量应取1而非-1提示：建议参考Zhou(2012)的公式形式，所有约束条件统一使用"+"号，方向向量取值0或1乘以原始变量
2.2 PuLP变量设置技巧
在构建线性规划问题时：
# 变量定义 self.betax = pulp.LpVariable.dicts( "scalingFactor_x", (self._i), lowBound=0, upBound=1 # 投入缩放比例β∈[0,1] ) self.betay = pulp.LpVariable.dicts( "scalingFactor_y", (self._r), lowBound=0 # 期望产出可无限扩张 ) self.betab = pulp.LpVariable.dicts( "scalingFactor_b", (self._s), lowBound=0, upBound=1 # 非期望产出最多削减至0 )2.3 权重向量的经济含义
ω向量决定各变量的重要性分配：
# 典型权重设置方案 weight = [ 1/9, 1/9, 1/9, # 三种投入要素权重相同 1/3, # 期望产出权重最大 1/9, 1/9, 1/9 # 三种非期望产出权重相同 ]3. 完整代码实现与调试
3.1 数据准备规范
建议使用pandas.DataFrame组织数据：
import pandas as pd inputs = pd.DataFrame({ '资本': [20,30,40,20,10,11,12,14], '劳动': [300,200,100,200,400,222,321,231], '能源': [20,30,40,20,10,11,12,14] }) outputs = pd.DataFrame({'GDP': [20,30,40,30,50,21,32,42]}) bad_outputs = pd.DataFrame({ 'SO2': [10,20,10,10,10,12,-2,-1], 'COD': [5,15,5,5,5,8,-1,0], 'CO2': [50,80,50,50,50,55,10,5] })3.2 模型初始化参数
关键参数组合示例：
directional_factor = [1,1,1,1,1,1,1] # 方向向量全取1 weight_vector = [1/9]*3 + [1/3] + [1/9]*3 # 权重分配 dea_problem = DEAProblem( inputs=inputs, outputs=outputs, bad_outs=bad_outputs, weight_vector=weight_vector, directional_factor=directional_factor, returns="VRS", # 可变规模报酬 disp="weak disposability" )3.3 结果解读与可视化
运行模型后：
results = dea_problem.solve() efficiency = pd.DataFrame.from_dict(results[1], orient='index', columns=['无效率值'])无效率值越大表示DMU效率越低，可通过matplotlib绘制效率分布：
import matplotlib.pyplot as plt plt.figure(figsize=(10,6)) efficiency.plot(kind='bar', legend=False) plt.title('各DMU无效率值分布') plt.ylabel('无效率值β') plt.xlabel('DMU编号') plt.grid(axis='y')4. 进阶应用：绿色全要素生产率测算
4.1 跨期效率比较
基于全局技术计算各期效率后，可构建Luenberger生产率指数：
# 计算t到t+1期的生产率变化 GML = 0.5 * ((β_t - β_t1) + (β_t_global - β_t1_global))4.2 常见问题排查
调试中遇到的典型错误：
- 无可行解：检查方向向量符号是否与约束条件一致
- 效率值全为0：确认权重向量ω是否合理分配
- 结果波动异常：检查数据单位是否统一
4.3 效率分解框架
绿色全要素生产率可分解为：
- 技术效率变化(EC)
- 技术进步(TC)
- 规模效率变化(SEC)
实现代码需要保存各期前沿面信息进行对比分析。
在完成第一个完整测算后，我终于理解导师说的"模型是工具，经济含义才是核心"。NDDF给出的无效率值β需要结合具体投入产出变量解释，比如β=0.3意味着在现有技术下，该地区可以再减少30%的能源消耗和污染排放，同时保持GDP不变。这种定量结果对政策制定才有参考价值。