杜老师好，规模报酬可变参数下的结果是不是有点问题？

首先感谢杜老师的答复。您是指你的ado命令的代码？你的sbmeff中vrs部分的5个等式约束的代码如下：
beq=1 \ J(nx+ny+nb,1,0) \ 1
Aeq1=1,J(1,nx,0),1:/((ny+nb)*Y),1:/((ny+nb)*B),J(1,rows(Xref),0)
mat l Aeq1
Aeq2=-X',I(nx),J(nx,ny+nb,0),Xref'
Aeq3=-Y',J(ny,nx,0),-I(ny),J(ny,nb,0),Yref'		 
Aeq4=-B',J(nb,nx+ny,0),I(nb),Bref'
Aeq5=J(1,1+nx+ny+nb,0),J(1,rows(Xref),1)

我买的红兰数据的matlab的这几个约束矩阵的代码如下：
f1=[zeros(1,NT) -(1/m)*DX(i+(t-1)*n,:) zeros(1,s1) zeros(1,s2) 1];
Aeq1=[X(:,1:NT) eye(m) zeros(m,s1) zeros(m,s2) -X(:,i+(t-1)*n);
     Y(:,1:NT) zeros(s1,m) -eye(s1) zeros(s1,s2) -Y(:,i+(t-1)*n);
     Z(:,1:NT) zeros(s2,m) zeros(s2,s1) eye(s2) -Z(:,i+(t-1)*n);
     zeros(1,NT) zeros(1,m) (1/(s1+s2))*DY(i+(t-1)*n,:) (1/(s1+s2))*DZ(i+(t-1)*n,:) 1;
     ones(1,NT) zeros(1,m+s1+s2) -1];
beq1=[zeros(m,1);
      zeros(s1,1);
      zeros(s2,1);
      1;0];
主要区别就在于那个lambda求和等于1的vrs约束矩阵是1在前，还是0在前。
您的是：Aeq5=J(1,1+nx+ny+nb,0),J(1,rows(Xref),1)，0在前
红兰的是：ones(1,NT) zeros(1,m+s1+s2) -1]，1在前。

我认为，您的aeq2-4分别对应他的aeq1的1-3行，您的aeq1对应他的第4行，您的aeq5对应他的aeq的最后一行。
最大的差别就是1在前还是0在前。
还望杜老师指点。

杜老师好，在您解释的基础上，我又认真地对比了一下代码。然后，我迷茫了。我想把对比的情况给您说说，然后也有一个小的请求，希望您能帮我一下。
由于过程比较长，我先说结论和我的请求。
结论就是，对比后，我觉得好像如您所说，好像两边的代码都没有错。于是，我怀疑是我当初测试的时候，stata的参数弄错了什么，所以我又测试了一篇，还是发现，只有vrs结果对不上。
因此，我有一个想法，其实，早以前就想，但因为stata的ado命令是不是还需要编译？我本来想改一下您的sbmeff，然后看看结果，但我做不到。所以，不知道杜老师能不能给我编译一个，把Aeq5=J(1,1+nx+ny+nb,0),J(1,rows(Xref),1)中，这两项调整位置后的sbmefftest命令给我，我再测试一下呢？（如果行的话，我的邮箱LazingLuo@foxmail.com）
非常感谢。

下面，就是我对比的一些过程。因为，我不懂stata语法，也不懂matlab语法，还看不太明白很复杂的数学推导，所以，您看看我下面对比的，想的有没有什么问题。只有vrs结果不一样的原因到底在哪里。

具体如下。
第一，红兰数据的matlab的求解过程，应该是和知乎上一篇文章一样，把目标函数做了一个Charnes-Cooper变换。我是参考了知乎上https://zhuanlan.zhihu.com/p/598312233，这篇文章的推导。但我不确定您的代码是否使用了这个变换。
第二，在规模报酬可变部分，我在GitHub上看到的您的几个关键变量的代码如下
c=(1,-1:/(nx*X),J(1,ny+nb+rows(Xref),0))
beq=1 \ J(nx+ny+nb,1,0) \ 1
Aeq1=1,J(1,nx,0),1:/((ny+nb)*Y),1:/((ny+nb)*B),J(1,rows(Xref),0)
Aeq2=-X',I(nx),J(nx,ny+nb,0),Xref'
Aeq3=-Y',J(ny,nx,0),-I(ny),J(ny,nb,0),Yref'		 
Aeq4=-B',J(nb,nx+ny,0),I(nb),Bref'
Aeq5=J(1,1+nx+ny+nb,0),J(1,rows(Xref),1)

红兰数据的matlab的代码如下：
f1=[zeros(1,n) -(1/m)*DX(i+(t-1)*n,:) zeros(1,s1) zeros(1,s2) 1];
beq1=[zeros(m,1) zeros(s1,1) zeros(s2,1) 1 0];
Aeq1=[X(:,n*(t-1)+1:n*t) eye(m) zeros(m,s1) zeros(m,s2) -X(:,i+(t-1)*n);
                Y(:,n*(t-1)+1:n*t) zeros(s1,m) -eye(s1) zeros(s1,s2) -Y(:,i+(t-1)*n);
                Z(:,n*(t-1)+1:n*t) zeros(s2,m) zeros(s2,s1) eye(s2) -Z(:,i+(t-1)*n);
                zeros(1,n) zeros(1,m) (1/(s1+s2))*DY(i+(t-1)*n,:) (1/(s1+s2))*DZ(i+(t-1)*n,:) 1;
                ones(1,n) zeros(1,m+s1+s2) -1];

简单，说明一下matlab代码中几个变量的含义。m就是投入变量的数量，等同于您的nx，s1就是产出变量的数量，等同于您的ny，s2就是非期望产出变量的数量，等同于您的nb，XYZ分别为投入、产出和非期望产出。

下面，我把代码一一对应一下。
先是目标函数
stata：c=(1,-1:/(nx*X),J(1,ny+nb+rows(Xref),0))
matlab：f1=[zeros(1,n) -(1/m)*DX(i+(t-1)*n,:) zeros(1,s1) zeros(1,s2) 1];
不知道是不是做了变换的问题，目标函数的区别有点大，但也有一些相似的地方。比如，J(1,ny+nb+rows(Xref),0)和zeros(1,n) zeros(1,s1) zeros(1,s2)三个部分，都是那么多个0；-1:/(nx*X)和-(1/m)*DX(i+(t-1)*n,:)也类似。只是这些矩阵值的顺序不太一样。

然后是等式约束，您的代码是分开写的，matlab的代码是合并写的，但完全是一一对应的。
第一，Aeq1=1,J(1,nx,0),1:/((ny+nb)*Y),1:/((ny+nb)*B),J(1,rows(Xref),0)
对应倒数第二行，zeros(1,n) zeros(1,m) (1/(s1+s2))*DY(i+(t-1)*n,:) (1/(s1+s2))*DZ(i+(t-1)*n,:) 1。
您的第一项“1”对应下面的最后一项“1”；
您的最后一项，J(1,rows(Xref),0)，对应下面的第一项，zeros(1,n)；
中间部分正好一一对应：J(1,nx,0)对应zeros(1,m)；
                  1:/((ny+nb)*Y)对应(1/(s1+s2))*DY(i+(t-1)*n,:)；
                  1:/((ny+nb)*B)对应(1/(s1+s2))*DZ(i+(t-1)*n,:)；

似乎好像是stata中把等式约束的左右两边调换了一个位置。下面的几个等于约束也是类似的。具体如下。

第二，Aeq2=-X',I(nx),J(nx,ny+nb,0),Xref'
对应第一行，X(:,n*(t-1)+1:n*t) eye(m) zeros(m,s1) zeros(m,s2) -X(:,i+(t-1)*n)
您的第一项，-X'，对应下面的最后一项，-X(:,i+(t-1)*n)；
您的最后一项，Xref'，对应下面的第一项，X(:,n*(t-1)+1:n*t)
中间部分正好一一对应：I(nx)对应eye(m)；
                  J(nx,ny+nb,0)对应zeros(m,s1) zeros(m,s2)

对到这里，我大概明白您说的「参照单元我是放在后面，也就是lambda决策变量是放在后面的」。
就是说，您是把X(:,n*(t-1)+1:n*t)或Xref'放在了“+”后面，把-X'或-X(:,i+(t-1)*n)放在了第一项。

第三，Aeq3=-Y',J(ny,nx,0),-I(ny),J(ny,nb,0),Yref'
对应第二行，Y(:,n*(t-1)+1:n*t) zeros(s1,m) -eye(s1) zeros(s1,s2) -Y(:,i+(t-1)*n);
您的第一项，-Y'，对应下面的最后一项，-Y(:,i+(t-1)*n)；
您的最后一项，Yref'，对应下面的第一项，Y(:,n*(t-1)+1:n*t)
中间部分正好一一对应：J(ny,nx,0)对应zeros(s1,m)；
                  -I(ny)对应-eye(s1)
                  J(ny,nb,0)对应zeros(s1,s2)

第四，Aeq4=-B',J(nb,nx+ny,0),I(nb),Bref'
对应第三行，Z(:,n*(t-1)+1:n*t) zeros(s2,m) zeros(s2,s1) eye(s2) -Z(:,i+(t-1)*n);
您的第一项，-B'，对应下面的最后一项，-Z(:,i+(t-1)*n)；
您的最后一项，Bref'，对应下面的第一项，Z(:,n*(t-1)+1:n*t)
中间部分正好一一对应：J(nb,nx+ny,0)对应zeros(s2,m) zeros(s2,s1)；
                  I(nb)对应eye(s2)

第五，即规模报酬可变的约束，Aeq5=J(1,1+nx+ny+nb,0),J(1,rows(Xref),1)
对应最后一行，ones(1,n) zeros(1,m+s1+s2) -1
显然，这里您是的两个部分，matlab的是三个部分。因为matlab的处理方式是将等于右边都移到左边，右边保留一个0。这点体现在beq中。那么，按上面的逻辑，好像对应关系也没有错。
您的第一项，是没有，有话理论上应该是“-1”，对应下面的最后一项“-1”
您的最后一项，J(1,rows(Xref),1)，对应下面的第一项，ones(1,n)
剩下的就是您的J(1,1+nx+ny+nb,0)，对应下面的，zeros(1,m+s1+s2)

最后就剩下beq了。
您的是：beq=1 \ J(nx+ny+nb,1,0) \ 1
matlab的是：beq1=[zeros(m,1) zeros(s1,1) zeros(s2,1) 1 0];
我对照了您在crs中的beq=1 \ J(nx+ny+nb,1,0)，显然，最后的那个1是对应规模报酬可变约束的1。
您的第一个1，正好对应matlab中后面的1。
matlab中最后的0是对应规模报酬可变约束，和上面一样将原来的1移到了等于同一边。

这样看起来，我大概明白您的意思，感觉好像也没有错。
但我还是在想，规模报酬可变的约束为∑λ=1，这里只有λ的事，是不是不存在放前面后面的问题呢？

我刚刚又测试了一次，发现还是crs结果一样，vrs不一样。那问题一定就出现在新增的这条约束条件上。
原来我是想把您的sbmeff代码修改后，再算一次的。但我发现，好像stata的ado命令也是需要编译的，是吗？
我在安装好的命令那看，就找不到对应的代码了。所以，我也修改不了，也测试不了。

不知道杜老师能不能给我编译一个，把Aeq5=J(1,1+nx+ny+nb,0),J(1,rows(Xref),1)中，这两项调整位置后的sbmefftest命令给我，我再测试一下呢？
非常感谢。

kerrydu/gtfpch

内容风险标识

评论 (9)

kerrydu/gtfpch .gitee-modal { width: 500px !important; }

内容风险标识