各位大侠，能不能帮忙看看我的程序，为什么总是出不了结果呢，谢谢各位。
new;closeall; cls;  external proc indices;
clearg der, der2, der22, c, c1, se1, se2, ss,H,R,U,vvv,j2,j12,
       c2,se,Eq1,Eq2,Eq3,Eq4,prob1,kk,EP,xx,xx2;
let pricevar = income work pphhsize ownhouse black other hispanic colledge ppmsacat pphhhead;

let exogvar = ppage pphhhead pphhsize income midwest west south married colledge male black other hispanic work ownhouse ppmsacat spring summer autumn;
datafile="c:\\gauss 9.0\\data7";
save path=D:\gauss;
load path=D:\gauss;
j1=.51        /*step length*/;
j2=0;         /*Hessian controller*/;
j3=0;j4=0;    /*other weights*/;
j6=01;        /*results display controller*/;
j8=01         /*step length controller: 1=adjust, 0=fixed */;
j20=0;
critic=0.01;  /*convergent criteria*/;
@***************************@;
/***Starting values ***/;

numvar_x=rows(exogvar)+2;
value_x=zeros(numvar_x,1);
numvar_p=rows(pricevar)+2;
value_p=zeros(numvar_p,1);
alpha=0;
s11=2;
s12=0;
s22=1;
b=value_x|value_p|alpha|s11|s22|s12;

/******************read in data************************/;
n=10675;
load data[n,31]=data7.txt;
      

numobs=rows(data);
t=numobs;
st=1;
caseid=data[st:n,1];
ppage=data[st:n,2];
pphhhead=data[st:n,3];
pphhsize=data[st:n,4];
ppdualin=data[st:n,5];
ppmsacat=data[st:n,6];
q5expend=data[st:n,7];
texpend=data[st:n,8];
tpounds=data[st:n,9];
income=data[st:n,10];
work=data[st:n,11];
ownhouse=data[st:n,12];
northeast=data[st:n,13];
midwest=data[st:n,14];
south=data[st:n,15];
west=data[st:n,16];
married=data[st:n,17];
colledge=data[st:n,18];
male=data[st:n,19];
spring=data[st:n,20];
summer=data[st:n,21];
autumn=data[st:n,22];
winter=data[st:n,23];
white=data[st:n,24];
black=data[st:n,25];
other=data[st:n,26];
hispanic=data[st:n,27];
price=data[st:n,28];
ap6gr2_1=data[st:n,29];
ap6gr2_2=data[st:n,30];
ap6gr2_3=data[st:n,31];

        
         
        lnprice=ln(price);
        purdum=tpounds .gt 0;
        nopurchase=1-purdum;
        y=tpounds;
        
        demand_rhs=ones(numobs,1)~data[.,2]~data[.,3]~data[.,4]~data[.,10]~data[.,28]~data[.,14]~data[.,16]~data[.,15]~data[.,17]~data[.,18]~data[.,19]~data[.,25]~data[.,26]~data[.,27]~data[.,11]~data[.,12]~data[.,6]~data[.,20]~data[.,21]~data[.,22];
        price_rhs=ones(numobs,1)~data[.,10]~data[.,3]~data[.,4]~data[.,11]~data[.,12]~data[.,25]~data[.,26]~data[.,27]~data[.,18]~data[.,6];
        
        
        demand_rhs[.,2 4 5 6]=ln(demand_rhs[.,2 4 5 6]);
        price_rhs[.,2 3 4]=ln(price_rhs[.,2 3 4]);
        
numbetas=cols(demand_rhs);
numgammas=cols(price_rhs);
numpar=rows(b);
k=numpar;
cc=1;
numdis1=trunc(rows(b)./7);
remdis1=rows(b)-(numdis1*7);


PROC LLf(b);
clearg var_x,var_p,cov_xp,betab,gammab,nn,alpha,za,xb,
       f0,error_x,error_p,error_xp,sigma,LL,rho;
        
  nn=numbetas+numgammas;
  betab=b[1:numbetas,.];           
  gammab=b[numbetas+1:nn];
  alpha=b[nn+1];
  
if j20 == 1;
/*
  sigma=upmat(xpnd(b[nn+2:nn+4]));
  sigma=sigma'sigma;
  var_x=sigma[1,1];
  var_p=sigma[2,2];
  cov_xp=sigma[1,2];
*/
  
  var_x=b[nn+2]^2;
  var_p=b[nn+3]^2;
  cov_xp=b[nn+4];
  sigma=(var_x~cov_xp)|(cov_xp~var_p);

endif;
  
if j20 ==0;  
  var_x=b[nn+2]^2;
  var_p=b[nn+3]^2;
  cov_xp=0;
  sigma=(var_x~cov_xp)|(cov_xp~var_p);   
endif;
  za=demand_rhs*betab;
  xb=price_rhs*gammab;
  
  error_x=y-za-alpha*lnprice;
  error_p=lnprice-xb;
  error_xp=error_x~error_p;
  
  f0=cdfn((-za-alpha*xb)./sqrt(var_x+alpha*alpha*var_p+2*alpha*cov_xp));
     
  LL=purdum.*(-ln(2*pi)-0.5*ln(det(sigma))
              -0.5*sumc(((error_xp*invpd(sigma)).*error_xp)'))+
     (1-purdum).*(ln(f0));   
                                 
RETP(LL);
ENDP;