绿色全要素生产率，python，采用全局生产技术的弱可处置性非径向方向距离函数(NDDF)，可调方向权重，DDF，DEA

---

一、NDDF是什么？

1. 采用方法

一文详细说明SBM、SBM-DDF、DDF、NDDF、ML指数是什么
利用python的pulp库进行CCR、BCC、超效率模型的数学建模
在本文使用的公式是采用全局生产技术的弱可处置性非径向方向距离函数
李江龙（2018）

由上图可知，存在 $\omega 、\lambda 、\beta 、g$ 这三个变量
其中 $\lambda$ 计算前沿面，可以理解为测算投入产出权重以获得一个前沿面
$\beta$ 为各个投入产出缩放比例
$\omega$ 为权重向量(denotes the normalized weight vector that is relevant to the numbers of inputs and outputs)
$g$为方向向量，当$g^x=x^t,g^y=y^t$则为(DDF)。特例如上文 ，$g=(0,0,-E,Y,-D,-S,-W)$
$\overrightarrow{ND}$得出的是无效率程度，值越大越无效率

2.具体参数解释

问题：上文公式问题，投入应该是减少而到达前沿面，${\sum }_{i=1}^N{\sum }_{t=1}^T{\lambda }_{i,t}{E}_{i,t}$为前沿面，$E-{\beta }_E{g}_E$当期减少多少投入能达到前沿面，于是带入方向向量后解释不通。
$$\begin{gathered} 文中方向向量定义为：g=(0,0,-E,Y,-D,-S,-W) \\ 但公式中出现：\sum _{i=1}^N\sum _{t=1}^T{\lambda }_{i,t}{E}_{i,t}\le E-{\beta }_E{g}_E \\ 将方向向量带入得：\sum _{i=1}^N\sum _{t=1}^T{\lambda }_{i,t}{E}_{i,t}\le E+{\beta }_{E}E \end{gathered}$$

Zhou（2012）提到得公式更好理解，此公式中所有约束条件都为“+”，因此带入方向向量时候是正确的。

文中提到 “If the direction vector is set equal to (-F, 0, 0), Eq. (7) is essentially equivalent to an input-oriented DEA model with undesirable outputs.”
$$\begin{gathered} g_1=(-F,0,0) \\ \overrightarrow{D_1}(F,E,C;g_1)=max{w}_{1F}{\beta }_{1F} \\ s.t.\sum _{n=1}^N{z}_{1n}{F}_{1n}\le (1-{\beta }_{1F})F \\ s.t.\sum _{n=1}^N{z}_{1n}{E}_{1n}\le E \\ s.t.\sum _{n=1}^N{z}_{1n}{C}_{1n}=C \\ {z}_{1n}\ge 0,n=1,2,...,N,{\beta }_{1F}\ge 0 \end{gathered}$$

由此可见当有方向向量时，公式中应都设为“+”号。若公式设好，则方向向量应取值为$0或1$而没有$-1$，再乘以原本的投入产出值。

Chung（1997）

当 $k=1,m=1,t=1;y_{k^{\prime }m}^t$ 意味着第一个时期的第一个DMU的第一个期望产出。

3.强弱可处置性

参考
Y. H. Chung, R. F ̈ are and S. Grosskopf.Productivity and Undesirable Outputs: A Directional Distance Function Approach.*Journal of Environmental Management * (1997) 51, 229–240
成刚. 数据包络分析方法与MaxDEA软件[M].北京：知识产权出版社，2014.5

二、代码

1.参考与改进

参考代码：
python DEA:强/弱处置性假设下的考虑非期望产出的非径向距离函数
改变部分

        else:			# 每一个方向应该要乘以原本变量
            self.gx = directional_factor[:self.I] * self.inputs
            self.gy = directional_factor[self.I:(self.I+self.R):1] * self.outputs
            self.gb = directional_factor[(self.I+self.R):(self.I+self.R+self.S):1] * self.bad_outs

按照上文参考文献，在约束条件中已确定了方向
那么方向向量部分得到的应该是 $[0\times L,0\times K,1\times E,1\times D,1\times S,1\times W]$

---

import numpy as np
import pandas as pd
import pulp

class DEAProblem:
    def __init__(
        self,
        inputs,
        outputs,
        bad_outs,
        weight_vector,
        directional_factor=None,
        returns="CRS",
        disp="weak disposability",
        in_weights=[0, None],
        out_weights=[0, None],
        badout_weights=[0, None],
    ):
        self.inputs = inputs
        self.outputs = outputs
        self.bad_outs = bad_outs
        self.returns = returns
        self.weight_vector = (
            weight_vector  # weight vector in directional distance function
        )
        self.disp = disp
        self.directional_factor = directional_factor
        
        # directional_factor 方向向量需要输入列表形式，并且若不等于0，则用1，公式方向已给定。
        # 投入产出传入数据皆为DataFrame，权重和方向向量为列表形式
        
        self.J, self.I = self.inputs.shape  # no of DMUs, inputs
        _, self.R = self.outputs.shape  # no of outputs
        _, self.S = self.bad_outs.shape  # no of bad outputs
        self._i = range(self.I)  # inputs
        self._r = range(self.R)  # outputs
        self._s = range(self.S)  # bad_output
        self._j = range(self.J)  # DMUs
        if directional_factor == None:  #（方向向量未定义则与原本投入产出相同）
            self.gx = self.inputs
            self.gy = self.outputs
            self.gb = self.bad_outs
        else:			# 每一个方向应该要乘以原本变量
            self.gx = directional_factor[:self.I] * self.inputs
            self.gy = directional_factor[self.I:(self.I+self.R):1] * self.outputs
            self.gb = directional_factor[(self.I+self.R):(self.I+self.R+self.S):1] * self.bad_outs


        self._in_weights = in_weights  # input weight restrictions
        self._out_weights = out_weights  # output weight restrictions
        self._badout_weights = badout_weights  # bad output weight restrictions

        # creates dictionary of pulp.LpProblem objects for the DMUs
        self.dmus = self._create_problems()

    def _create_problems(self):
        """
        Iterate over the DMU and create a dictionary of LP problems, one
        for each DMU.
        """

        dmu_dict = {}
        for j0 in self._j:
            dmu_dict[j0] = self._make_problem(j0)
        return dmu_dict

    def _make_problem(self, j0):
        """
        Create a pulp.LpProblem for a DMU.
        """

        # Set up pulp 设置变量
        prob = pulp.LpProblem("".join(["DMU_", str(j0)]), pulp.LpMaximize)
        self.weights = pulp.LpVariable.dicts(
            "Weight", (self._j), lowBound=self._in_weights[0]
        )
        self.betax = pulp.LpVariable.dicts(
            "scalingFactor_x", (self._i), lowBound=0, upBound=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
        )

        # Set returns to scale
        if self.returns == "VRS":
            prob += pulp.lpSum([self.weights[j] for j in self.weights]) == 1

        # Set up objective function
        prob += pulp.lpSum(
            [(self.weight_vector[i] * self.betax[i]) for i in self._i]
            + [(self.weight_vector[self.I + r] * self.betay[r]) for r in self._r]
            + [
                (self.weight_vector[self.I + self.R + s] * self.betab[s])
                for s in self._s
            ]
        )
# Set up constraints  # 括号外j0为该函数传入的j0
        # Set up constraints
        for i in self._i:
            prob += (
                pulp.lpSum(
                    [(self.weights[j0] * self.inputs.values[j0][i]) for j0 in self._j]
                )
                <= self.inputs.values[j0][i] - self.betax[i] * self.gx.values[j0][i]
            )
        for r in self._r:
            prob += (
                pulp.lpSum(
                    [(self.weights[j0] * self.outputs.values[j0][r]) for j0 in self._j]
                )
                >= self.outputs.values[j0][r] + self.betay[r] * self.gy.values[j0][r]
            )

        if self.disp == "weak disposability":
            for s in self._s:  # 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":
            for s in self._s:  # 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]
                )

        return prob

    def solve(self):
        """
        Iterate over the dictionary of DMUs' problems, solve them, and collate
        the results into a pandas dataframe.
        """

        sol_status = {}
        sol_weights = {}
        sol_efficiency = {}
        sol_objective = {}
        sol_directional ={}
        if self.directional_factor == None:
            for ind, problem in list(self.dmus.items()):
                problem.solve()
                sol_status[ind] = pulp.LpStatus[problem.status]
                sol_weights[ind] = {}
                sol_objective[ind] = problem.objective
                for v in problem.variables():
                    sol_weights[ind][v.name] = v.varValue
                sol_efficiency[ind] = pulp.value(problem.objective)
            return sol_status, sol_efficiency, sol_weights, sol_objective
        else:
            for ind, problem in list(self.dmus.items()):
                problem.solve()
                sol_status[ind] = pulp.LpStatus[problem.status]
                sol_weights[ind] = {}
                sol_objective[ind] = problem.objective
                sol_directional[ind] = self.directional_factor
                for v in problem.variables():
                    sol_weights[ind][v.name] = v.varValue
                sol_efficiency[ind] = pulp.value(problem.objective)
            return sol_status, sol_efficiency, sol_weights, sol_objective, sol_directional

sol_status 每个DMU是否计算成功
sol_efficiency 得出方程的解，测量的是无效率值
sol_weights Weight求得是$\lambda$ ; scalingFactor_x、 scalingFactor_y、 scalingFactor_b求的是目标函数中的$\beta$
sol_objective 目标函数的设定情况
sol_directional方向向量的设定情况

1.1约束条件关键代码解释

        for i in self._i:
            prob += (
            pulp.lpSum([(self.weights[j0] * self.inputs.values[j0][i]) for j0 in self._j]) #不要求设定时间，遍历所有时期的所有dmu，将所有变量都累加。i是第i个投入变量。
            <= self.inputs.values[j0][i] - self.betax[i] * self.gx.values[j0][i]) #这里 j0 是函数_make_problem参数中的j0，是上面函数第j0个变量。

2.示例

X = pd.DataFrame(
    np.array(
        [
            [20, 300,20],
            [30, 200,30],
            [40, 100,40],
            [20, 200,20],
            [10, 400,10],
            [11, 222,11],
            [12, 321,12],
            [14, 231,14],
        ]
    )
)
y = pd.DataFrame(np.array([[20], [30], [40], [30], [50], [21], [32], [42]]))
b = pd.DataFrame(np.array([[10], [20], [10], [10], [10], [12], [-2], [-1]]))
weight = [0, 0,1 / 3 ,1 / 3, 1 / 3]
directional = [0,0,1,1,1]
solve1 = DEAProblem(X, y, b, weight, directional_factor=directional, disp="weak disposability").solve()

3.实例演示

G = pd.read_excel(r'F:\绿色全要素测算.xlsx',sheet_name='Sheet1')
G.head(50)

其中k、l、e是投入要素(inputs)，y是期望产出(output)，s、f、w是非期望产出(bad output)

input1 = G.iloc[:,0:3]
input1
output = G.iloc[:,3:4]
output
b_output = G.iloc[:,4:7]
b_output
weight = [1/9,1/9,1/9,1/3,1/9,1/9,1/9] #令投入要素、期望产出和非期望产出比例相同
directional_facto = [1,1,1,1,1,1,1] #全设为1，因为在公式中已经定好(-K,-L,-E,Y,-S,-F,-W)投入要素和非期望产出是越少越好。
solve11 = DEAProblem(input1,output,b_output,weight_vector=weight,directional_factor=directional_facto,disp="weak disposability").solve()

status = pd.DataFrame.from_dict(solve11[0], orient="index", columns=["status"])
efficiency = pd.DataFrame.from_dict(solve11[1], orient="index", columns=["efficiency"])
weights = pd.DataFrame.from_dict(solve11[2], orient="index")
results = pd.concat([status, efficiency, weights], axis=1)
results.to_excel(r'F:\绿色全要素测算1.xlsx') 将结果输出

此为得出的效率值

此为目标函数的$\beta$值

需要说明一点是该代码得出的为单个dmu在全局生产技术下所得到的非效率值。在结果中看不出时期和第几个dmu，因此若有分时期和dmu位数可以后期手动加上。因为得出的结果是根据输入投入产出的顺序得到的，所以按输入时的时期和dmu在结果前添加两列即可
后续会根据该模型继续改进

三、绿色指标

邵帅（2022）利用Luenberger 生产率进行定义t+1 期的碳排放绩效（TFCEP）为

若前一期的无效率值大于后一期，则说明全要素碳排放生产率得到了改善
杨翔（2015）利用全域GMalmquist-Luenberger测算碳生产率

李江龙（2018）

四、非全局生产技术的弱可处置性非径向方向距离函数(NDDF)

张宁（2022）采用两期模型，同样利用两期非径向Luenberger生产率指数定义碳要素生产率。

---

参考文献

1. 李江龙，徐斌. 诅咒”还是“福音”:资源丰裕程度如何影响 中国绿色经济增长?[J]. 经济研究，2018(9)
2. 邵帅，范美婷，杨莉莉. 经济结构调整、绿色技术进步 与中国低碳转型发展 ——基于总体技术前沿和空间溢出效应视角的经验考察[J].管理世界, 2022(2)
3. 杨翔，李小平，周大川. 中国制造业碳生产率的差异与收敛性研究[J]. 数量经济技术经济研究，2015(12)
4. 张宁. 碳全要素生产率、低碳技术创新和节能 减排效率追赶* ———来自中国火力发电企业的证据[J]. 经济研究，2022(2)
5. P. Zhou, B.W. Ang, H. Wang., 2012. Energy and CO2 emission performance in electricity generation: A non-radial directional distance function approach. European Journal of Operational Research 221, 625-635
6. Y. H. Chung, R. F ̈ are and S.,1997. Grosskopf.Productivity and Undesirable Outputs: A Directional Distance Function Approach.Journal of Environmental Management 51, 229–240
7. 成刚. 数据包络分析方法与MaxDEA软件[M].北京：知识产权出版社，2014.5