Undefined symbol G0025
没有定义的符合。

最好把程序贴上来

不好意思，程序较长，分几次发吧，错误用红色表示出来，在5楼和6楼，帮我看看吧，非常感谢啊

@ This Gauss program simulates draws from Bayesian posterior distribution
          using  importance sampling @
    output file = junk reset;
    library pgraph,optmum;
    graphset; optset;  

format 16,8;

rndseed 9137841;   @ resets random number generator so that identical sample is used
                                    each time @

kex = 4;  @ use kex = 1 for example 1
                           kex = 2 for example 2
                           kex = 3 for example 3
                           kex = 4 for oil price data @
if kex < 4;
     #include readdat;
else;
     #include ..\data2001\reoil01.g;
endif;

@ set parameters that describe the data  @
    k = cols(x);                @ k is the number of nonlinear variables @
    klin = 1+cols(xlin);         @ klin is the total number of variables, nonlinear, linear
                                              and constant term @
    xwhole = ones(n,1)~xlin;   @ xwhole is (n x klin) @
  

/* ===========================================================
             CALCULATE AND REPORT SOME BASIC STATISTICS
=============================================================*/
xbar = meanc(x);
sigx = sqrt(meanc((x - xbar')^2));
"number of observations";;n;
"means of nonlinear explanatory variables";xbar';
"standard deviations";sigx';
xlinbar = meanc(xlin);
siglinx = sqrt(meanc((xlin - xlinbar')^2));
"means of all explanatory variables";xlinbar';
"standard deviations";siglinx';
ybar = meanc(y);
sigy = sqrt(meanc((y-ybar)^2));
"mean and standard deviation of dependent variable is";
ybar;;sigy;

xx0 = xwhole'*xwhole;
xx0inv = invpd(xx0);

/* ===========================================================
             SET INITIAL VALUES FOR GLOBAL VARIABLES
==============================================================*/
 kqopt = 2;
          @ kqopt = 1 means evaluate covariance by iterating as in
                      Theorem 2.2
            kqopt = 2 means evaluate covariance directly from Table 1
             option 2 is faster bust only works for k <= 5 @
   kmle = 1;    @ kmle = 1 will call numerical optimization routines to find MLE
                           kmle = 2 will skip this step @
   kc = 2;
          @ kc = 2 echos parameter values when proc is evaluated
            kc =1 produces no echo @
    gamx = ones(k,1);
          @ gamx can be used to restrict some of the weighting coefficients
            to be zero, if desired @
 
   ks = 0;
          @ ks = # of smoothed inferences desired
               ks = 0 means only evaluate likelihood function @
   xs = 0;
          @ xs will be a (ks x k) matrix of values for the vector x at which
            mean is to be evaluated for smoothed inference @
  xswhole = 0;
          @ xswhole will be a (ks x klin) matrix for all explanatory variables including
               constant term at which mean is to be evaluated for smoothed inference @
 _pdate = "";

/* ===========================================================
             INCLUDE NEEDED PROCEDURES
=============================================================*/
ll = 0;
proc(2) = lsq(yq,x1);
   @ This proc performs an OLS regression of yq on x1
        and puts value of log likelihood in the global scalar ll @
local b,eps,sigols,sigmle,ss,nq;
  b = invpd(x1'*x1)*x1'*yq;
   eps = yq - x1*b;
  sigols = sumc(eps^2)/(rows(x1)-cols(x1));
   sigmle = sumc(eps^2)/rows(x1);
   ss = sigols*invpd(x1'*x1);
   ll = -(n/2)*(1 + ln(2*pi)) - (n/2)*ln(sigmle);
 if kc > 1;
    ""; "ols coefficients";b';
     "ols see";;sqrt(sigols);
    "mle see";;sqrt(sigmle);
    "standard errors";sqrt(diag(ss))';
    "ols log likelihood";; ll;"";
endif;
retp(b,sigmle^.5);
endp;

 {bhat,sigmle} = lsq(y,xwhole);
ll0 = ll;        @ this saves ll0 for use in adjusting f(y,theta) by scale of f(y) @

#include proccov;
#include procs;
 #include fullproc;

proc (0) = graph1dim(nofx,sigmult,xfix);
     @ this proc fills the (ks x nk) matrix xs and the (ks x klin) matrix xswhole
          with values needed to calculate the
          effect of changing variable nofx with all others fixed at xfix.
          The variable nofx is varied from its value in xfix plus or minus
           sigmult times its standard deviation @
      ks = 61;   @ ks is the number of function evaluations to be performed @
      xs = ones(ks, klin - 1) .* xfix';
      xs[.,nofx] = seqa(xfix[nofx,1] - sigmult*sigx[nofx,1],
                                   (2*sigmult*sigx[nofx,1])/(ks-1), ks);
      xswhole = ones(ks,1) ~ xs;
      xs = xs[.,1:k];
endp;

proc (3) = graph2dim(nofx1,nofx2,sigmult,thx,xfix);
     @ this proc calculates the effects of changing variables nofx1 and nofx2 over
          plus or minus sigmult times their standard deviations while holding other
          variables constant at xfix. The proc returns x3 (the different values for
          variable nofx1 considered), y3 (the different values for variable nofx2
          considered), and z3 (values of the function mu(x).
          Graph output with command contour(x3',y3,z) @
    local iks, x3, y3, z3, cjunk;
      ks = 61;
      z3 = zeros(ks,ks);
      xs = ones(ks, klin - 1).*xfix';
      iks = 1;
      do until iks > ks;
         xs[.,nofx1] =  xfix[nofx1,1]-sigmult*sigx[nofx1,1] +
              ((iks-1)*(2*sigmult*sigx[nofx1,1])/(ks-1))*ones(ks,1) ;
         xs[.,nofx2] = seqa(xfix[nofx2,1]-sigmult*sigx[nofx2,1],
                                        (2*sigmult*sigx[nofx2,1])/(ks-1),ks);
         xswhole = ones(ks,1) ~ xs;
         xs = xs[.,1:k];
         if kmle == 1;
                    z3[.,iks] = conc(thx);
          else;
                   cjunk = full(thx);
                   z3[.,iks] = cjunk[.,1];
          endif;
       iks = iks + 1;
      endo;
      x3 =seqa(xfix[nofx1,1]-sigmult*sigx[nofx1,1],
                      (2*sigmult*sigx[nofx1,1])/(ks-1),ks);
       y3 = seqa(xfix[nofx2,1]-sigmult*sigx[nofx2,1],
                      (2*sigmult*sigx[nofx2,1])/(ks-1),ks);
retp(x3,y3,z3);
endp;

/* ===========================================================
             PERFORM OLS REGRESSION AND TEST FOR NONLINEARITY
=============================================================*/

 {bhat,sigmle} = lsq(y,xwhole);
"";"Testing the null hypothesis of linearity with respect to x";
    "Note-- it is maintained under both null and alternative that model ";
   "is linear with respect to elements of xlin other than x";
ehat = y - xwhole*bhat;
zeta = lm2(ehat,xlin,x); 

if kex ==4;
   "";"Testing the null hypothesis that model with quantity shifts is linear";
   xzwhole = ones(n,1)~qzquant~xlin[.,xp+1:cols(xlin)];
   {vhat,sigv} = lsq(y,xzwhole);
    vhat = y - xzwhole*vhat;
    zeta = lm2(vhat,xzwhole[.,2:cols(xzwhole)],x);

    "";"Testing the null hypothesis that Mork's specification is correct";
    xzwhole = ones(n,1)~xzmork~xlin[.,xp+1:cols(xlin)];
    {vhat,sigv} = lsq(y,xzwhole);
    vhat = y - xzwhole*vhat;
    zeta = lm2(vhat,xzwhole[.,2:cols(xzwhole)],x);

    "";"Testing the null hypothesis that Lee, Ni, Ratti specification is correct";
    xzwhole = ones(n,1)~xlee1~xlin[.,xp+1:cols(xlin)];
    {vhat,sigv} = lsq(y,xzwhole);
    vhat = y - xzwhole*vhat;
    zeta = lm2(vhat,xzwhole[.,2:cols(xzwhole)],x);

    "";"Testing the null hypothesis that Hamilton's specification is correct";
    xzwhole = ones(n,1)~xzham~xlin[.,xp+1:cols(xlin)];
    {vhat,sigv} = lsq(y,xzwhole);
    vhat = y - xzwhole*vhat;
    zeta = lm2(vhat,xzwhole[.,2:cols(xzwhole)],x);

    "";"Testing the null hypothesis that 3-year-high specification is correct";
    xzwhole = ones(n,1)~x3ham~xlin[.,xp+1:cols(xlin)];
    {vhat,sigv} = lsq(y,xzwhole);
    vhat = y - xzwhole*vhat;
    zeta = lm2(vhat,xzwhole[.,2:cols(xzwhole)],x);

     /*i = 1;
 do until i > n;
 "";"leave out ";;qz[i+8,1];
 qblip = zeros(n,1);
 qblip[i,1] = 1;
 {vhat,sigv} = lsq(y,xzwhole~qblip);
 vhat = y - (xzwhole~qblip)*vhat;
 zeta = lm2(vhat,xzwhole[.,2:cols(xzwhole)]~qblip,x);
 i = i + 1;
 endo; */

endif;

[此贴子已经被作者于2009-5-15 9:49:29编辑过]

/* ===========================================================
             FIND MAXIMUM LIKELIHOOD ESTIMATES, IF DESIRED
=============================================================*/

kc = 1;
if kmle == 1;    @ this starts the MLE section @

/* ------------------------------------------------------------------------------------------------------------------------
   === SET STARTING VALUES FOR NUMERICAL OPTIMIZATION === */

th0 = ( 1./(sigx*sqrt(k)) ) | 0.5;
th0';

  
 @ thx = 0.084 | 0.001 | 1.85 ; h = eye(k+1); i = 1; psi = zeros(1,k+2); goto jump4;  @

proc startval;     @ This defines starting value for iteration to be th @
  retp(th0); endp;

/* -----------------------------------------------------------------------------------------------------------------------------
   === SET GLOBAL CONTROL PARAMETERS FOR NUMERICAL
          OPTIMIZATION === */

     #include optmum.ext;
      _opalgr = 2;    @ This chooses BFGS optimization  @
      _opmiter = 150;  @ This controls the maximum number of iterations @
      _opshess = (1.e-4)*eye(k+1);

/* -----------------------------------------------------------------------------------------------------------------------------
   === CALL GAUSS PROCEDURES FOR NUMERICAL OPTIMIZATION === */

        kc=1;
        {thx,f,g,h} =optmum(&conc,startval);

"";"";"Final results";"MLE as parameterized for numerical optimization ";
"Coefficients:";thx';
"";"Value of log likelihood:";;-f;
"";"Gradient vector:";g';

"==================================";
kc = 2;
call conc(thx);

/* ------------------------------------------------------------------------------------------------------------------------
   === CALCULATE STANDARD ERRORS === */
kc = 3;
ks = 0;
thmore = conc(thx);

kc = 1;
g = gradp(&full,thx|thmore);
"";"Numerically calculated gradient:";g;
"";"";
h = (hessp(&full,thx|thmore));
   va = eigrs(h);
   if minc(eigrs(h)) <= 0;
        "Negative of Hessian is not positive definite";
        "Either you have not found local maximum, or else estimates are up "
        "against boundary condition.  In latter case, impose the restricted "
        "params rather than estimate them to calculate standard errors";
    else;
       h = invpd(h);
       std = diag(h)^.5;
       "For vector of coefficients parameterized as follows,";(thx | thmore)';
       "the standard errors are";std';
    endif;
"variance covariance matrix is";
format /m3;
h;
format /m1;
h;

/* -------------------------------------------------------------------------------------------------------------------------
   === GRAPH FUNCTIONS OF INTEREST === */

_pgrid = 3 | 0;

if kex == 1;
   @ plot figure 4 for MLE @
      _ptek = "graph4ml.tkf";
      title("Figure 4 (MLE)");
      call graph1dim(1,2,xbar);
      chis = conc(thx);
      y1hat = 0.6*xs[.,1].*(xs[.,1].>zeros(ks,1)) + 0.2*xbar[2,1]*ones(ks,1);
      xy(xs[.,1],chis~y1hat);
   @ plot figure 5 for MLE @
      _ptek = "graph5ml.tkf";
      title("Figure 5 (MLE)");
      call graph1dim(2,2,xbar);
      chis = conc(thx);
      xy(xs[.,2],chis);

elseif kex == 2;
  /*    _pgrid = 0 | 0;
      _ptek = "graph6.tkf";
      title("Figure 6");
      {x3,y3,z3} = graph2dim(1,2,2,thx,xbar);
      contour(x3',y3,z3);  */
     
elseif kex == 3;
      _ptek = "graph7ml.tkf";
      title("Figure 7 (MLE)");
      call graph1dim(3,2,xbar);
      chis = conc(thx);
      xy(xs[.,3],chis);
endif;

if kex == 3;
    @ This section calculates values used in Figure 8 @
    kc = 1;
    ks = 61;
     "Conditional mean at different dates";
    xs = seqa(xbar[1,1]-2*sigx[1,1],(4*sigx[1,1])/(ks-1),ks)~
    (ones(ks,1)*(-0.7))~(ones(ks,1)*1955);
    chis = conc(thx);
    xs~chis;

  "";
  xs = seqa(xbar[1,1]-2*sigx[1,1],(4*sigx[1,1])/(ks-1),ks)~
    (ones(ks,1)*(12.3))~(ones(ks,1)*1975);
  chis = conc(thx);
  xs~chis;

  "";
  xs = seqa(xbar[1,1]-2*sigx[1,1],(4*sigx[1,1])/(ks-1),ks)~
    (ones(ks,1)*(3.9)~ones(ks,1)*1985);
  chis = conc(thx);
  xs~chis;
endif;

endif;    @ this ends the MLE section @

/* ============================================================
               PERFORM IMPORTANCE SAMPLING, IF DESIRED
=============================================================*/

/* ============================================================
              SET INITIAL VALUES FOR PARAMETERS
   ============================================================ */

jump4:

@ ---------------------- parameters of prior for theta ------------------------------------------ @

meta = 0;                         @ meta is mean of log of eta @
seta = 1.0;                       @ seta is std deviation of log of eta @
mg = ln(1./(sigx*sqrt(k)));  @ mg is k x 1 vector of prior means for ln(g) @
taug = 1.0*ones(k,1);          @ taug is prior std. deviation for ln(g) @
mth = mg | meta;
sigth = taug | seta;

@ ---------------------- parameters of prior for psi ------------------------------------------------ @
nsig = 0.25;                         @ prior is (1/sigma^2) ~ gamma(nsig, lamsig) @
lamsig = nsig*0.5*(sigy^2);
hbet = n;    
mbet = zeros(klin,1);
mbet[1,1] = meanc(y);

      @ prior is beta ~ N(mbet,hbet*sigma^2*invpd(xwhole'*xwhole))  @

@ ----------parameters that determine importance sampling density -------------- @
montdf = 2;                  @ montdf is degrees of freedom for student t
                                         generated importance  density @
montfact = sqrt(2);                  @ montfact is factor by which standard deviations are
                                                    increased to obtain std. deviation of importance density @
pprob = 0.5;          @ pprob = probability of drawing from the Student t versus
                                  the spread-out prior  @

"";"prior distributions:";
"    ln(eta) ~ Normal(";;meta;;",";;seta;"^2)";
ii = 1;
do until ii > k;
"    ln(g(";;ii;;")) ~ Normal(";;mg[ii,1];;",";;taug[ii,1];;"^2)";
ii = ii+1;
endo;
"    (1/sigma^2) | eta ~ Gamma(";;nsig;;",";; lamsig;;")";
"    beta | sigma,theta ~ N(mbet,hbet*sigma*invpd(xwhole'*xwhole)) ";
"    mbet:";;mbet';
"    hbet:";;hbet;
"";

"Importance density:";
"    degreees of freedom of t distribution";;montdf;
"    fraction of observations from t distribution";;pprob;
"    factor by which std. deviation of importance density exceeds prior";;montfact;

@ ----------------------- parameters to control monte carlo runs -----------------------@

nmonte = 20000;         @ nmonte is number of monte carlo draws generated @
"Number of monte carlo draws";;nmonte;
"";

if kmle == 1;
   thmle = abs(thx);
   pmix = h[1:k+1,1:k+1];
else;
    "you must input a mean and variance for the student t mixing distribution";
endif;

kmle = 2;
pmix = pmix*montfact^2;

uxmix = zeros(nmonte,1);    @uxmix can be used to keep track of
                                                 which component of importance density was used @

@ --------  calculation of constant terms for importance density ------------------- @

tprec = invpd(pmix);
tdet = detl;
tc = gamma( (montdf+k+1)/2 )/gamma(montdf/2);
tc = tc / (   sqrt(tdet) * (montdf * 3.14159)^( (k+1) / 2)  );
nc =1./( sqrt(2*3.14159).*sigth);   @ note that nc is (k +1) x 1 vector @
lnnc =ln(nc);

pmix = chol(pmix);

/* =============================================================
             INCLUDE NECESSARY PROCEDURES
   ============================================================= */

psi = 0; i = 0;                        @ needed to avoid tripping GAUSS error compiling code @
#include bayproc2;

proc normgam(sigtry,bettry,ngam,lamgam,betsy,varsy);
    @ this proc evaluates the log of the product of a gamma(ngam,lamgam) density
     for (1/sigtry) and a N(betsy,sigtry*varsy) density for bettry at the
     points sigtry and bettry @
 local f1,f2,m1,m2;
        if ngam > 10;
          f1 = ngam*ln(lamgam) - (ngam - 1)*ln(sigtry) - (lamgam/sigtry) - lnfact(ngam-1);
        else;
           f1 = ngam*ln(lamgam) - (ngam - 1)*ln(sigtry) - (lamgam/sigtry) - ln(gamma(ngam));
        endif;
 m1 = invpd(varsy);
 m2 = detl;
 f2 = (-klin/2)*ln(2*3.1415927*sigtry) -(1/2)*ln(m2);
        f2 = f2 - (1/(2*sigtry))*(bettry-betsy)'*m1*(bettry-betsy);
retp(f1 + f2);
endp;

proc priadj(th);
     @ this proc adds log of prior p(theta) to log of f(y|X,theta) to
          arrive at log of f(y,theta|X) @
local q;
if ndpchk(16);
       "underflow before calling priadj";
endif;
   q = -((ln(th) - mth)^2)./(2*(sigth^2));
   q = q + lnnc - ln(th);
   q = sumc(q);
if kc > 2;
      "log of prior is";;q;
endif;
q = q + baby(th); 
if ndpchk(16);
   "underflow after calling priadj";
endif;
retp(q);
endp; 

proc stut(nstu,ndif);
      @ this proc generates a (1 x nstu) vector student t with ndf degrees of freedom @
       local stuf,stux;
       stuf = rndn(1,nstu);
       stux = rndn(ndif,1);
       stux  = sumc(stux^2)/ndif;
        stux = sqrt(stux);
       stuf = stuf/stux;
retp(stuf);
endp;

proc mix(i);
      @ this proc generates mixture of prior and MLE-weighted student t variable
    global:
        thmle = (k+1) x 1 vector of ML estimates
 pmix = (k+1) x (k+1) matrix of cholesky decomposition of montfact^2 times
                     asymptotic  variance-covariance matrix
        pprob = probability of drawing from the Student t versus the spread-out prior
     output:
         prmix = (k+1) x 1 vector of variables generated from mixture of spread-out
                       Student t and prior @
local prmix,umix;
        
        umix = rndu(1,1);
       
         if umix < pprob;
              uxmix[i,1] = 1;
              prmix = pmix'*stut(k+1,montdf)' + thmle;
         else;
              uxmix[i,1] = 0;
              prmix = mth + montfact*sigth.*rndn(k+1,1);
              prmix = exp(prmix);
        endif;
retp(prmix);
endp;

proc impor(thx);
      @ this proc calculates the value of the importance sampling density
           at the point thx @
local val1, val2,val;
    if ndpchk(16);
         "underflow before calling impor";
    endif;
    val1 = - ((ln(thx) - mth)^2) ./ (2 * ((montfact*sigth)^2) );
    val1 = val1 - ln(thx) + lnnc - ln(montfact);
    val1 = sumc(val1);
    val1 = exp(val1);
    if ndpchk(16);
          "underflow after calculating val1";
    endif;
    val2 = (thx - thmle)'*tprec*(thx - thmle);
    val2 = 1 + val2/montdf;
    val2 = tc*val2^( -(k+1+montdf)/2 );
    if ndpchk(16);
           "underflow after calculating val2";
    endif;
    val = pprob*val2 + (1 - pprob)*val1; 
    if ndpchk(16);
          "underflow after calling impor";
          "theta is";thx';
    endif;
    if kc > 2;
         "contribution of lognormal component to log importance density";;val1;
         "contribution of student t component to log importance density";;val2;
   endif;
retp(val);
endp;

thet = zeros(nmonte,k+1);
psi = zeros(nmonte,klin+1);

wtchk = zeros(nmonte,2);
thwt = ones(nmonte,1);

i = 1;
do until i > nmonte;
    tryagain:
           if ndpchk(16);
                    "underflow occurred prior to calling mix";
                    ndpclex;
              endif;
    thet[i,.] = mix(i)';
    if ndpchk(16);
                    "underflow occurred just after calling mix";
                    ndpclex;
              endif;

    if minc(thet[i,.]') < 0;
          goto tryagain;
    endif;

    beep1 = priadj(thet[i,.]') - ll0;
     if ndpchk(16);
                    "underflow occurred just after calling priadj";
                    "beep1 is";;beep1;
                    ndpclex;
              endif;
      beep2 = impor(thet[i,.]');
       if ndpchk(16);
                    "underflow occurred just after calling impor";
                    "beep2 is";;beep2;
                    ndpclex;
              endif;
     @ control for underflows @
      if beep1 < -500;
              thwt[i,1] = 0;
      else;
             thwt[i,1] = exp(beep1 - ln(beep2)  );
             wtchk[i,1] = exp(beep1);
             wtchk[i,2] = beep2;
 
        endif;
   
       if ndpchk(16);
                    "underflow occurred right after exp and ln";
                    "beep1 is";;beep1;
                    "beep2 is";;beep2;
                    ndpclex;
              endif;
i = i + 1;
endo;

if ndpchk(16);
    "underflow occurred prior to tabulating monte carlo";
    ndpclex;
endif;

nwt = sumc(thwt);
"sum of weights for Monte Carlo Draws is";;nwt;
"If this is small, there may be a numerical problem with the simulation";

" estimated mean for theta parameters";
thm = sumc (thet .* thwt);
thm = thm/nwt;
thm';
"estimated standard errors";
thv = (thet - thm')'*( (thet - thm') .* thwt) ;
thv = diag(thv/nwt);
sqrt(thv)';
"standard error of estimated mean";
thu = ( (thet - thm') .* thwt )'*( (thet - thm') .* thwt );
thu = diag(thu);
thu = thu ./ (  (sumc(thwt))^2 );
sqrt(thu)';
"Geweke's measure of relative efficiency";
thv'./(nmonte * thu)';

"";" estimated mean for psi parameters";
psim = sumc (psi .* thwt);
psim = psim/nwt;
psim';
"estimated standard errors";
psiv = (psi - psim')'*( (psi - psim') .* thwt) ;
psiv = diag(psiv/nwt);
sqrt(psiv)';
"standard error of estimated mean";
psiu = ( (psi - psim') .* thwt )'*( (psi - psim') .* thwt );
psiu = diag(psiu);
psiu = psiu ./ (  (sumc(thwt))^2 )';
sqrt(psiu)';
"Geweke's measure of relative efficiency";
psiv'./(nmonte * psiu)';

hchk = thwt~wtchk~thet~seqa(1,1,nmonte);
hchk = sortc(hchk,1);
hchk[.,1] = hchk[.,1]/nwt;
hchk[.,2] = hchk[.,2]/sumc(hchk[.,2]);
hchk[.,3] = hchk[.,3]/sumc(hchk[.,3]);
"";"Fraction of distribution accounted for by:";
"          Most influential observation";;hchk[nmonte,1];
"          50 most influential observations";;sumc(hchk[nmonte-49:nmonte,1]);
"";"Fifty most influential observations are as follows";
"weight   f(y,thet)   I(thet)   thet  # in simulation";
format 12,4;
hchk[nmonte-49:nmonte,.];
format 16,8;

if ndpchk(16);
    "underflow occurred as a result of tabulating monte carlo";
    ndpclex;
endif;

@ free up some memory space @
clear hchk, qz*,uxmix,wtchk,xham,xlee1,xmork,xzham,xzmork,xzwhole,_op*;  （line686）出错指示

[此贴子已经被作者于2009-5-15 9:37:21编辑过]

/* ===========================================================
                GRAPH FUNCTIONS OF INTEREST
  ============================================================*/
sigstep = 2;     @ sigstep is the multiple of standard deviations used in graphing
                             inferred function @
critsig = 0.05;   @ call to proc seesmooth will calculate (1 - critsig)*100 percent
                              confidence intervals @

proc(3) = seesmooth(critsig);
@ this proc generates smoothed inferences and 95% confidence intervals for the
function evaluated at the points contained in the global ks x k matrix xs @
local i, j, itest, isum, cjunk, chis, chierr, chisavg, xlow, xhigh;
chis = zeros(ks,nmonte);   @ col i contains smoothed estimates for generated
                                               parameter from monte carlo run i @
chierr = zeros(ks,nmonte); @ col i contains smoothed estimates for generated
                                                parameter from monte carlo run i plus measurement
                                                inference noise @

i = 1;
do until i > nmonte;
     cjunk = full( (thet[i,.] ~ psi[i,.])');
      chis[.,i] = cjunk[.,1];
      chierr[.,i] = chis[.,i] + sqrt(cjunk[.,2]).*rndn(ks,1);
      if ndpchk(16);
               "underflow occurred analyzing realization # ";;i;
               "theta and psi vectors";
               thet[i,.];
               psi[i,.];
              ndpclex;
     endif;
i = i+1;
endo;
chisavg = (chis*thwt)/nwt;

xlow = zeros(ks,1);
xhigh = zeros(ks,1);
j = 1;
do until j > ks;
   cjunk = sortc( (chierr[j,.]' ~ thwt) ,1);
    itest = 1;
    isum = 0;
    start1:
     if (isum/nwt) < (critsig/2);
           isum = isum + cjunk[itest,2];
           itest = itest + 1;
           goto start1;
     endif; 
   xlow[j,1] = cjunk[itest,1];
   itest = nmonte;
   isum = 0;
   start2:
   if (isum/nwt) < (critsig/2);
           isum = isum + cjunk[itest,2];
           itest = itest - 1;
           goto start2;
   endif;
   xhigh[j,1] = cjunk[itest,1];
j = j+1;
endo;
retp(chisavg,xlow,xhigh);
endp;

/*--------------------------- GRAPH EFFECTS OF VARIABLES #1 AND #2 ------------------------*/
if kex == 1;
ks = 61;
call graph1dim(1,sigstep,xbar);  @ this sets xs so that variable 1 changes
                                                        with variable 2 fixed @
{chisavg,xlow,xhigh} = seesmooth(critsig);
let _pltype[4,1] = 4 6 3  3;
ytrue = 0.6*(xs[.,1].*(xs[.,1].>zeros(ks,1))) + 0.2*xs[.,2];   （line762）出错信息
title("Figure 4");
_ptek = "fig4.tkf";
xy(xs[.,1],ytrue~chisavg~xlow~xhigh);

call graph1dim(2,sigstep,xbar);  @ this sets xs so that variable 2 changes
                                                         with variables 1 and 2 fixed @
{chisavg,xlow,xhigh} = seesmooth(critsig);
let _pltype[3,1] =  6 3  3;
title("Figure 5");
_ptek = "fig5.tkf";
xy(xs[.,2],chisavg~xlow~xhigh);

/*---------------------------  GRAPH CONTOUR PLOT ----------------------------------- */
elseif kex == 2;
ks = 61;
z3ex = zeros(ks,ks);  （line778）出错信息
izz = 1;
do until izz > nmonte;
"izz=";;izz;
     {x3z,y3z,z3z} = graph2dim(1,2,2,(thet[izz,.]~psi[izz,.])',xbar);
     z3ex = z3ex + thwt[izz,1]*z3z;
izz = izz+1;
endo;
z3ex = z3ex/nwt;
 _pgrid = 0 | 0;
 _ptek = "graph6.tkf";
 title("Figure 6");
contour(x3z',y3z,z3ex);

/*---------------------------- GRAPH EFFECTS OF TIME ---------------------------------- */
elseif kex == 3;
ks = 61;
call graph1dim(3,sigstep,xbar);  @ this sets xs so that variable 3 varies with variable 2 fixed @
{chisavg,xlow,xhigh} = seesmooth(critsig);

"";"Effect of time";
"    Bayes  lower  upper";
xjunk = chisavg ~xlow ~ xhigh;
xjunk;"";
let _pltype[3,1] =  6 3  3;
_ptek="fig7.tkf";
title("Figure 7");
xy(xs[.,3],chisavg~xlow~xhigh);

/*---------------------------- GRAPH EFFECTS OF LAGS 3 AND 4 OF OIL -------------------------------- */
elseif kex == 4;
scale(-15 | 16, -1 | 2);
@ graph effects of lag 3 for 3 different values of lag 4 @
_ptek = "lags3_4.tkf";
title("Effect of lag 3 for alternative lags 4");
ks = 61;
   xgoose = xlinbar;
   xgoose[4,1] = -5;
   call graph1dim(3,1.5,xgoose);
   {dn5,xlow,xhigh} = seesmooth(critsig);

   xgoose = xlinbar;
   xgoose[4,1] = 0;
   call graph1dim(3,1.5,xgoose);
   {dn0,xlow,xhigh} = seesmooth(critsig);

   xgoose = xlinbar;
   xgoose[4,1] = 5;
   call graph1dim(3,1.5,xgoose);
   {up5,xlow,xhigh} = seesmooth(critsig);

"";"Effect of o(t-3) when     o(t-4)=-5     o(t-4)=0    o(t-4)=5";
xs[.,3]~dn5~dn0~up5;
"";
let _pltype[3,1] = 1 6 3;
xy(xs[.,3],dn5~dn0~up5);

ifig = 1;
do until ifig > xp;
      _ptek = chrs(48+ifig) $+ "bafig.tkf";
       title("Bayes Lag " $+chrs(48+ifig));

ks = 61;
call graph1dim(ifig,1.5,xlinbar);  @ this sets xs so that variable ifig varies with
                                                         others fixed @

{chisavg,xlow,xhigh} = seesmooth(critsig);

"";"Effect of variable";; ifig;
"    Bayes  lower  upper";
xjunk = chisavg ~xlow ~ xhigh;
xjunk;"";
let _pltype[3,1] =  6 3  3;
xy(xs[.,ifig],chisavg~xlow~xhigh);

ifig = ifig+1;
endo;

endif;