格林第五版16.6例子，中自助法估计。谁能分析我的程序一下对不对？
%lad函数来源于Jplv7:James P. LeSage

clc;clear
samplecode_path ='D:\MatlabClass\格林5\格林7网站数据\TableF7-2.csv';%读入csv数据。
[ndata,txt1,MixedData1]=xlsread(samplecode_path);
Y=ndata(:,1);
K=ndata(:,2);
L=ndata(:,3);
N=ndata(:,4);
X=[ones(length(ndata),1) log(K./N) log(L./N)];
y=log(Y./N);
[b,bint,r,rint,stats] =regress(y,X);

ee=r'*r;
sig_2=ee/(rows(X)-cols(X(:,1:end)));%式子4-7s_2=ee/（n-k）;s_2的平方根就是回归标准误。本例为428-5
x=X(:,1:end);
b_se=sqrt(diag(sig_2*inv(x'*x)));%见第五版P53页。
t=b./b_se;%第五版4-14P55页。
c_mat=sig_2*inv(x'*x);%第五版49-50页。
b0=b;
b1=b0;
%b1=[1;2;3];
XX=[y X];
X=XX%XX(1:4,:)%;
%%%%LAD估计
b= lad(y,XX(:,2:4))
%%%%以下或不对。
NR=500;
ind=unidrnd(25,25,NR);
X0=X(:,2:4);
q=y;
BL = lad(q,X0) ;
bl=BL.beta
VB = zeros(3,3);%Init(3,3,0.)
for i=1:NR
   NR;
    aaaa=ind(:,i);
    x=X0(ind(:,i),:);
    BLr = lad(q,x);
    blr=BLr.beta;
    d = blr - bl ; %原文为D。
    VB = VB + 1/NR * d*d';
end
b_se=sqrt(diag(VB))