报错如下：
ERROR: Invocation of unresolved module VAR_DUYC(B).
ERROR: NLPQN call: Error in argument module HES.
ERROR: (execution) Unknown or error operation executed.

希望哪位好心的大神能帮忙分析解决，很着急用，谢谢！因为程序很长，所以打字将程序分段落的目的描述一下，这是一个非线性优化估计参数的问题。
首先，constraints总讲明了含待估参数的矩条件
其次，在目标函数中，是要实现含待估参数的目标函数最小化，其中目标函数中用到了在

程序代码如下：
**estimatioin;

proc iml;

ma         =1;        /*Choose if you want MA(1) transitory income shock*/
        
load moments;
m        =moments[,1];                                        /*Vector of covariances*/
var =moments[,2:ncol(moments)-1];        /*Weighting matrix*/

t                =moments[1,ncol(moments)];        /*Lenth of first difference*/
intyr        =moments[3,ncol(moments)];
endyr        =moments[4,ncol(moments)];
k             =ma+1;

**********************************************;

**define the constraints;

start var_duyc(b);        
        
        theta =b[k-1];
        vxi          =b[k];
        vzt   =b[1+k:k+t1];
        veps  =b[1+k+t1:k+t1+t2];
        phi          =b[1+k+t1+t2];
        psi   =b[1+k+t1+t2+1];

var_duy =j(t,t,0);
var_duc =j(t,t,0);
cov_duyc=j(t,t,0);
var_duyc=j(2*t,2*t,0);

/*
var_duy[i,j]: cov(var_duy,var_duy[j]);
vzt:          variance of zeta (permanent shocks);
veps:          variance of epsilon (transitory shocks);
vxi:           variance of xi (heterogeneity of consumption)
phi:           parameter of permanent shocks;
psi:          parameter of transitory shocks;
theta:    serial correlation of transitory shocks.
*/

**matrix of Var(duy): var_duy;
**j<=2;
var_duy[1,1]=vzt[1]+veps[1]+(1-theta)**2*veps[1]+theta**2*veps[1];
var_duy[2,2]=vzt[2]+veps[2]+(1-theta)**2*veps[1]+theta**2*veps[1];
var_duy[1,2]=-(1-theta)**2*veps[1];

j=3;
do while (j<=t);
        var_duy[j,j]=vzt[j]+veps[j]+(1-theta)**2*veps[j-1]+theta**2*veps[j-2];
        var_duy[j-1,j]=-(1-theta)*veps[j-1]+theta*(1-theta)*veps[j-2];
        var_duy[j-2,j]=-theta*veps[j-2];
        j=j+1;
end;

i=2;
do while (i<=t);
        j=i;
        do while (j<=t);
                var_duy[j,i-1]=var_duy[i-1,j];
                j=j+1;
        end;
        i=i+1;
end;

**matrix of Var(duc): var_duc;
j=1;
do while (j<=t);
        var_duc[j,j]=phi**2*vzt[j]+psi**2*veps[j]+vxi;
        j=j+1;
end;

**matrix of Cov(duc, duy): cov_duyc;
j=1;
do while (j<=t);
        cov_duyc[j,j]        =phi*vzt[j]+psi*veps[j];
        j=j+1;
end;

j=2;
do while (j<=t);
        cov_duyc[j-1,j]=phi*(1-theta)*veps[j-1];
        j=j+1;
end;

j=3;
do while (j<=t);
        cov_duyc[j-2,j]=-phi*theta*veps[j-2];
        j=j+1;
end;

var_duyc[1:t,1:t]                 =var_duc;
var_duyc[t+1:2*t,1:t]         =cov_duyc;
var_duyc[1:t,t+1:2*t]         =t(cov_duyc);
var_duyc[t+1:2*t,t+1:2*t]=var_duy;
var_duyc_vct                         =vech(var_duyc);                /*transform the matrix to a vector*/
return(var_duyc_vct);
finish var_duyc;

*************************************************************************;

**define the objective function;

start f(b);
run var_duyc(b);
fmm        =var_duyc_vct;

*j(nrow(var),col(var),0)||;
peso=diag(var);
f         =t(m-fmm)*ginv(peso)*(m-fmm);
return(f);
finish f;

***********************************************************************;

*****************starting values for iteration;

    b0                      =j(1+k+t1+t2+1,1,0);
    b0[k-1]                 =-0.1;               /*theta*/
    b0[k]                   =0.01;               /*variance of xi*/
    b0[1+k:k+t1]             =0.03*j(t1,1,1);     /*variance of the perm shocks*/
    b0[1+k+t1:k+t1+t2]        =0.03*j(t2,1,1);     /*variance of the trans shocks*/
    b0[1+k+t1+t2]           =1;                                /*coeff of partial insurance perm shock*/
    b0[1+k+t1+t2+1]                =0;                                /*coeff of partial insurance of trans shock*/

optn= j(1,10,.);
CALL NLPQN(rc,b,"f",b0,optn) nlc="var_duyc(b)";