/***************************************************************************/
/*                                                                         */
/*   Title: DOITC.G                                                        */
/*          includes doitc (procedure), callopt (sub-procedure)            */

/* OLD version saved 3/19 before elimination of obsfact.                   */
/*                                                                         */
/* Written: John H. Cochrane                                               */
/*      by: University of Chicago Graduate School of Business              */
/*          1101 E. 58th Street, Rosenwald 2xx                             */
/*          Chicago, Ill 60637                                             */
/*          ph#: (312) 702-3059                                            */
/*          e-mail: john.cochrane@gsb.uchicago.edu                         */
/*                                                                         */
/* Purpose: DOITC sets up the data, factors, and weighting matrix to       */
/*          do the iterated search of the technology &/or consumption      */
/*          parameters, plus estimation of additional factors.  CALLOPTC   */
/*          is the procedure, called by DOITC, to do the actual iteration. */
/*          CALLOPTC uses GAUSS application program OPTMUM to search for   */
/*          optimal parameter values.                                      */
/*                                                                         */
/*                                                                         */
/*                                                                         */
/***************************************************************************/

/* ------------------------------------ */
/* Part I: (1) Call GAUSS library files */
/*         (2) Reset graphic globals    */
/* ------------------------------------ */

library pgraph,optmum,gmm;
graphset;
#include optmum.ext;

/* -----------------------------------*/
/* Part II: Procedures                */
/* -----------------------------------*/

/* ------------------------------------------------------------------------- */
/* CALLOPTC.G                                                                */
/*                                                                           */
/* USAGE:                                                                    */
/*  {alpha,delta,mpk,jt,covs,b,f,allr,ik,reti} =                             */
/*  callopt(alpha,delta,mpk,selpar,invest,rassetx,xtrfact,inst,rf,W,howto,b);*/
/*                                                                           */
/* INPUTS: (see description under doitc proc)                                */
/*OUTPUTS: (see description under doitc proc)                                */
/*                                                                           */
/* -------------------------------------------------------------------------- */

proc(16) = calloptc(alpha,delta,mpk,gam,selpar,invest,rassetx,xtrfact,inst,
                                                         rf,W,howto,b);
    local ik,tpar,jt,gradnt,retcode,reti,m,f,lgt,covs,allr,ccode,neg_m;
    ccode = 0;
    if maxc(howto[2].==6|7|8|9);
        ccode = 1;                             /* is consumption in here */
    endif;

    ik = makeik(delta,invest);                         /* create i/k series */
    tpar = encodec(alpha,delta,mpk,gam,selpar,ccode);
    call passtof(invest,rassetx,xtrfact,inst,
                       rf,ik,W,selpar,howto,alpha,delta,mpk,b);

    /* use optmum to find production function parameters that minimize JT */
    if 1;
        output off;
        {tpar,jt,gradnt,retcode} = optmum(&findjtc,tpar);
        output on;
        if retcode /= 0 ;
            "Optmum terminated abnomally. retcode=" retcode;
        endif;
    endif;


   if 0;
    /* diagnostic: print JT on a grid */

     clearg grider,i,j,tp1,tp2,jtmat;
     grider = (0.80|0.85|0.90|0.95|0.975|1.00|1.025|1.05|1.10|1.15|1.20);
     tp1 = grider*tpar[1];
     tp2 = grider*tpar[2];
     jtmat = zeros(rows(tp1),rows(tp2));
     i = 1;
     do while i <= rows(tp1);
         j = 1;
         do while j <= rows(tp2);
             jtmat[i,j] = findjtc(tp1[i]|tp2[j]);
             j = j+1;
         endo;
         i = i+1;
     endo;
     format /RD 10,4;
     "Results on a grid";
     100*((0~tp2')|(tp1~jtmat));
     graphset;
     _plegstr = ftos(grider[1],"%*.*lf",4,3);
     j = 2;
     do while j <= rows(grider);
         _plegstr = _plegstr $+ "\000" $+ ftos(grider[j],"%*.*lf",4,3);
         j = j+1;
     endo;
     _plegctl = 1;
     xy(100*tp2,100*jtmat);
     graphset;
     surface(tp2',tp1,jtmat);
 endif;


    {reti,ik,m,f,lgt,b,covs,allr,neg_m,w} = passbakc;
    {alpha,delta,mpk,gam} = decodec(tpar,selpar,ccode,alpha,delta,mpk,gam);
    retp(alpha,delta,mpk,gam,jt,covs,b,f,allr,ik,reti,tpar,neg_m,lgt,m,w);
endp;

/* ------------------------------------------------------------------------- */
/* DOITC.G                                                                   */
/*                                                                           */
/*  USAGE:                                                                   */
/*   {alpha1,delta1,mpk1,gam1,b1} =                                          */
/*      doit(alpha,delta,mpk,selpar,howto,invest,rassetx,xtrfact,inst,       */
/*                            rf,nwlag,step2,ergrf,invlbl,xtrlbl);  */
/* INPUTS:                                                                   */
/*  alpha,delta,mpk:  default parameter vectors, rows() = cols(invest).      */
/*  selpar:           3x1 vector, which of {alpha delta mpk} to search over  */
/*  howto:            3x1 vector controls levels vs. excess returns as factors*/
/*                     whether extra factors are included and cond/uc moments*/
/*  invest:           matrix of investment series 47:1 to 90:4               */
/*  rassetx:          matrix of excess asset returns 47:2 to 90:4            */
/*  xtrfact:          matrix of extra factors to try 47:2                    */
/*  inst:             matrix of instruments, 0 if just constant, 47:1 to 90:4*/
/*                           (program lags twice and chops returns to fit)   */
/*  rf:               risk free rate series to form excess inv. ret. 47:2    */
/*  nwlag:            scalar, # newey west lags to use                       */
/*  step2:            1: do second stage gmm, 0: don't                       */
/*  ergrf:            1: make graphs Er vs. pred. 1st stage|2nd|see?         */
/*  invlbl:           label conformable to cols(invest)                      */
/*  xtrlbl:           label conformable to cols(xtrafact)                    */
/*  xtrcode           factor code to indicate whether to use as RHS variable */
/*                        in OLS regression                                  */
/*                                                                           */
/* OUTPUTS:                                                                  */
/*  {alpha1,delta1,mpk1,b1,gam1} = first stage gmm; use for successsive calls*/
/*                                                                           */
/* Notes:                                                                    */
/* investment returns are automatically added to asset returns for testing.if*/
/* an extra factor is also a return, you must manually add it to rassetx     */
/*                                                                           */
/*                                                                           */
/*                                                                           */
/* ------------------------------------------------------------------------- */

proc(8) =
   doitc(alpha,delta,mpk,selpar,howto,invest,rassetx,xtrfact,inst,
                          rf,nwlag,step2,ergrf,invlbl,xtrlbl,xtrcode,
                          preder2,erx2,grphcode);

clearg ik,tpar,jt,gradnt,retcode,reti,m,f,lgt,b,label,W,s,
    D0,covs,varmat,se,trash,rix,T,erx,erix,covrf,sig,sigin,predr,doff,pval,
    alpha1,delta1,mpk1,b1,ninv,nxtr,npar,nb,parvec,j,jtur,oldhow,chi2,blank,
    table,mask,fmt,allr,dev,preder,
    n,k,const,xx,b2,resid,s2,t,R2,y,select,ninst,cfact,dofobs,gam,gam1,
    stderr,i,olspval,j1,neg_m,tneg_m,junk,ej,pj,pval2,lgt,
    T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,lgt2,lgt3,lgt4,lgt5,lgt6,lgt7,
    lgt8,lgtjunk,deltas,mom2x,iteridx,rhow,lastw,oldchi,bT,allrT,tpT,trash,
    constlbl;

clearg term,varg,va,ve,npar,lam,firststg;

clearg _price;

 /* --------------------------------------------------------------- */
 /* 1. Generate output statements, given values passed by switches  */
 /* --------------------------------------------------------------- */

    tneg_m = 0;

    if rows(invlbl) /= cols(invest);
        "Fix investment label";
    endif;
    if ((howto[2] == 3) or (howto[2] == 5)) and (rows(xtrlbl) /= cols(xtrfact));
        "Fix xtra label";
    endif;

    if maxc(howto[2] .== 1|2|3|6|7) ;
        "Using unscaled factors";
    endif;
    if maxc(howto[2] .== 4|5);
        "using scaled factors";
    endif;
    if howto[2] == 8;
        "using scaled i returns, nonscaled consumption growth";
    endif;
    if howto[3] == 0;
        "Using unconditional momemts";
    endif;
    if howto[3] == 1;
        "Using scaled returns";
    endif;

    /* --------------------------------------------------- */
    /* 4.  do 1st-stage GMM with identity weighting matrix */
    /* --------------------------------------------------- */

    firststg = 1;                    /* flag used in std error, jt computation*/
    w = 0;                           /* findjt makes it an identity matrix */
    b = 0;                           /* only needs b when doing grads */
    gam = -1;                        /* hard wire starting value.  */

    {alpha,delta,mpk,gam,jt,covs,b,f,allr,ik,reti,tpar,neg_m,lgt,m,w} =
       calloptc(alpha,delta,mpk,gam,
                    selpar,invest,rassetx,xtrfact,inst,rf,W,howto,b);
    if ergrf[1]; gosub makgraf; endif;

    alpha1 = alpha;                          /* save first stage parameters */
    delta1 = delta;
    mpk1 = mpk;
    b1 = b;
    gam1 = gam;

   /* ----------------------------------------------- */
   /*  5. Print first stage parameter estimates, se's;*/
   /*     call finds0.g to construct weighting matrix */
   /* ----------------------------------------------- */

    "";
    "2. First stage GMM (Identity weighting matrix). ";
        gosub outpar;

    "";
    format /RD 4,0;
    "Using " nwlag "newey-west lags";
    s = findS0(f,nwlag);
    gosub calcse;
    gosub outse;
    gosub chisq;
    gosub overid;


   /* ------------------------------------------------------------------ */
   /* 6. Second-stage GMM estimation using the returned weighting matrix */
   /* ------------------------------------------------------------------ */

    firststg = 0;   /* flag to tell se, jt computations to use second stg */

    if step2;

if 0;
      /* ITERATED GMM CODE. CHOOSE THIS OR TWO STAGE ONLY THAT FOLLOWS */
      "";
      "3. Iterating GMM (optimal weighting matrix) ";

      lastw = 0;                      /*was eye(cols(f)); */
      lastw =    invpd(findS0(f,nwlag));
      rhow = .5;  /* w is updated by w = rho* new w + (1-rho) last w */

      format /RD 4,0;
      "Using " nwlag "newey-west lags";
      "";

      iteridx = 1;
      do while iteridx <= 100;

          w = rhow*invpd(findS0(f,nwlag)) + (1-rhow)*lastw;
          lastw = w;

          {alpha,delta,mpk,gam,jt,covs,b,f,allr,ik,reti,tpar,neg_m,lgt,m,w} =
              calloptc(alpha,delta,mpk,gam,selpar,
                  invest,rassetx,xtrfact,inst,rf,W,howto,b);

          format /RD 10,4;
          "Iterating GMM for " iteridx "th time; JT = " jt
              "%ch" 100*(jt-oldchi[1])/jt;

          if iteridx == 1;
              oldchi = jt;
          else;
              oldchi = jt|oldchi;
          endif;

          if iteridx > 5;
              if (maxc(abs(oldchi[1:4]-oldchi[2:5])./jt) < 0.05) ;
                  break;
              endif;
          endif;

          iteridx = iteridx + 1;
          if iteridx == 200; "iterated GMM did not converge"; endif;
      endo;

      /* END OF ITERATED CODE */

else;
      /* TWO STAGE ONLY       */
      "";
      "3. second stage gmm ";

      w =invpd(findS0(f,nwlag));

      {alpha,delta,mpk,gam,jt,covs,b,f,allr,ik,reti,tpar,neg_m,lgt,m,w}=
          calloptc(alpha,delta,mpk,gam,selpar,
              invest,rassetx,xtrfact,inst,rf,W,howto,b);
       /* END TWO STAGE ONLY CODE   */
endif;

      gosub outpar;    /* parameters */
      gosub calcse;
      gosub outse;
      gosub chisq;
      gosub overid;/* chi squared tests on overid, including comp.s*/
      if ergrf[2];
          gosub makgraf;
      endif;

    endif;

/* -------------------------------------------------------- */
/* 7.  OPTIONAL: Perform OLS of retunrs on the instruments  */
/* -------------------------------------------------------- */
/*
    if (inst /= 0) and (howto[2] == 3);
        n = rows(inst);                   /* no need to redo with scaled fact*/
        k = cols(inst);
        y = rassetx~reti;
        if xtrcode == 2;
            y = y~xtrfact;
        endif;
        const = ones(n,1);
        xx = inv((const~inst)'(const~inst));
        b2 = xx*((const~inst)'(y));
        resid = y-(const~inst)*b2;
        s2 = sumc(resid^2)/(n-k-1) ;
        se = (diag(xx)*s2')^(1/2);
        t = b2./se;
        R2 = 1 - (stdc(resid)^2)./(stdc(y)^2);

        format /RD 5,0;
        "Regressions of returns on instruments (investment returns last)";
        cols(rassetx) "asset excess returns";
        cols(invest) "investment returns (not excess)";
        if xtrcode == 2;
            cols(xtrfact) "extra factors" ;
        endif;
        format /RD 5,2;
        "Coefficients"; b2;
        "t"; t;
        "R2"; R2';
        "% p val " ;
        i = 1;olspval = zeros(rows(r2),1);
        do while i<=rows(r2);
             olspval[i] = 100*cdffc( R2[i]/(1-R2[i])*(n-k-1)/k , k ,n-k-1);
             i=i+1;
        endo;
        olspval';


     endif;
*/

    retp(alpha1,delta1,mpk1,gam1,b1,neg_m,preder,erx);

                              /* note returns FIRST stage GMM est's */


/* ---------------------------------- */
/* Part III: Subroutines              */
/* ---------------------------------- */

/* ----------------------------------------------------- */
/* A. MAKGRAF: Sets up the publication graphics settings */
/*             to construct E(r) vs. pred E(r) graphs    */
/* ----------------------------------------------------- */

makgraf:
    if howto[2] /= 6;
        @preder = - covs'*b ; old @
         preder = (meanc(m.*allr)-_price)/meanc(m) ;
         /*
         " Mean excess returns %    pricing error % " ;
         (100*meanc(allr))~(100*meanc(m.*allr)/meanc(m)) ;
         */
    endif;
    if howto[2] == 6;
        preder =
        - 1/rows(allr)*(xtrfact^gam-meanc(xtrfact^gam))'(allr-meanc(allr)') /
        meanc(xtrfact^gam);
        preder = preder';
    endif;

     Erx = 0|meanc(allr[.,2:cols(allr)]) ; /* assumes 1st col is rf itself */
     preder = erx+preder;

     graphset;
     if not ergrf[3];
        _Pscreen = 0;               /* just for production runs */
     endif;
     _protate = 1;
     margin(0,0,0,3);
     _Plctrl = -1|0;
     _pltype = 6;
    _plegctl = 1;
   @ _plegctl = { 2 4 2 0 };@
     ylabel("Mean Excess Return. % per qtr.");
     xlabel("Predicted Mean Excess Return, % per qtr.");
/*
     xtics(0,3,.5,0);
     ytics(0,3,.5,0);
*/
     /* option to allow error bars in graph  */
     if 0;
        stderr = stdc(allr)/(rows(allr)^.5);
        _perrbar = 100*(preder*ones(1,3)~Erx~(Erx-stderr)~(Erx+stderr))~
                    ((6~2~1).*ones(rows(preder),3));
     endif;
     if grphcode;
        _Plctrl = -1|-1|0;
        _pstype = 8|3;
        _psymsiz = 3|5;
        _plegstr = "Nonscaled Returns\000Scaled Returns\00045[o] Line";
        xy(100*fill(fill(preder2,preder),(0|maxc(erx))),
           100*fill(fill(erx2,erx),(0|maxc(erx))));
     endif;
     if not grphcode;
        _pstype = 3;
        _psymsiz = 3;
        _plegctl = 0;
        xy(100*fill(preder,(0|maxc(Erx))) , 100*fill(Erx,(0|maxc(Erx))));
     endif;

return;

/* --------------------------------------------------------------------- */
/* B. OUTPAR: Generate and format the output of the estimated parameters */
/* --------------------------------------------------------------------- */

outpar:

    /* set up numbers of things for later use */

    nb = rows(b);

    ninv = cols(invest);
    nxtr = 0;
    if maxc(howto[2] .== 3|5|7);
        nxtr = cols(xtrfact);        /* scaled use their own */
    endif;

   "ninst ninv nxtr nb"; ninst  ninv nxtr nb;


  /* coefficient estimates */

  if not maxc(howto[2] .== (6|7|8));    /* gave up on formatting cons model; */

    "Parameter estimates" ;
    let label = "alpha" "delta" "mpk";
    let constlbl = "const";
    let blank = " ";
                                       /* top label:    mpk mpk b b b b b */
    label = blank | (selif(label,selpar).*.ones(ninv,1)) |
                ("b" .* ones(nb,1));
                                      /* 2 line: NR R  NR R [xtr] NR R [xtr] */
                                      /*           mpk   level    scaled      */
    if howto[2] == 2;
        label = label ~ (blank | invlbl | constlbl | invlbl ) ;
    endif;
    if howto[2] == 3;
        label = label ~ (blank | invlbl | constlbl | invlbl | xtrlbl ) ;
    endif;
    if howto[2] == 4;
        label = label @~ (blank | invlbl | (ones(cols(inst)+1,1).*.invlbl))@;
    endif;
    if howto[2] == 5;
        label = @label ~ @ (blank | invlbl |
                                (ones(cols(inst)+1,1) .*. (invlbl|xtrlbl) ) );
    endif;

    parvec = selif(alpha|delta|mpk,selpar.*.ones(ninv,1));

    table = (label')|
            ((blank|(100*parvec)|b)');         /* mpk multiplied by 100 */
    mask = zeros(cols(label),rows(label))|
           (0~ones(1,rows((parvec|b))));
    fmt = "*.*lf"~12~3;
    call printfm(table,mask,fmt);
  endif;

  if maxc(howto[2] .== 6|7|8);  /* consumption model */
    "parameter estimates";
    if maxc(howto[2] .== 7|8);  /* consumption model */
        parvec = selif(alpha|delta|mpk,selpar.*.ones(ninv,1))|gam;
        table = (100*parvec[1:rows(parvec)-1]|parvec[rows(parvec)]|b)';
        mask =  ones(1,rows(parvec|b));
    endif;
    if howto[2] == 6;
        parvec = gam;
        table = (parvec)';
        mask =  ones(1,rows(parvec));
    endif;
    fmt = "*.*lf"~12~2;
    call printfm(table,mask,fmt);
    "";
  endif;

return;

/* ---------------------------------------------------------- */
/* C. OUTSE: Setup the tables for the se's & t-values         */
/* ---------------------------------------------------------- */

outse:
   if not howto[2] == 6;
    table = ("s.e."~((100*se[1:cols(invest)])|se[cols(invest)+1:rows(se)])' )|
            ("   t "~ ((parvec|b)./se)');
    mask = (0~ones(1,rows(se)))|
           (0~ones(1,rows(se)));
    fmt = "*.*lf"~8~2;
    call printfm(table,mask,fmt);

   endif;

   if howto[2] == 6;
     format /RD 8,2;
     "se" se;
     " t" tpar/se;
   endif;

return;

/* ---------------------------------------------------------- */
/* D. CALCSE: Construct the Chi-Square tests for joint b=0    */
/*            for all selected sets of the factors            */
/* ---------------------------------------------------------- */

calcse:

    /* ---------------------------------------- */
    /* 1. calculate standard errors a la Hansen */
    /* ---------------------------------------- */

    tpar = encodec(alpha,delta,mpk,gam,selpar,maxc(howto[2].==6|7|8));
    call passtof(invest,rassetx,xtrfact,inst,
                    rf,ik,W,selpar,2|howto[2:3],alpha,delta,mpk,b);
                                            /* note howto[1]=2: hold b fixed */
    D0 = gradp(&findgtc,tpar);
    call passbakc;  /* clear globals, save m etc. from actual GMM not grad */

    if not howto[2] == 6;
                                             /* add deriv wrt b */
        D0 = D0~covs';
    endif;
    T = rows(f);
                                 /* note vcov of transformed par's */
                                 /* for optimal w, give same result */
                                 /* done explicitly in step2 for comfort */
    if firststg;
        varmat = 1/T*inv(D0'W*D0)*D0'W*S*W*D0*inv(D0'W*D0);
    else;
        varmat = 1/T*inv(D0'W*D0);
    endif;
    se = diag(varmat)^(1/2);


    /* ----------------------------------------------- */
    /* 2.  pvalues for tests on joint parameter vector */
    /* ----------------------------------------------- */

    ninst = 0;
    if (howto[2] == 4 or howto[2] == 5);
        ninst = cols(inst);
    elseif cols(inst) > 1;
        "Error in joint be tests; thinks no instruments but there are";
    endif;

    if howto[2] == 6;                  /* no b's in c model , no tests */
        goto endjb;
    endif;

    /* ---- joint all production coeffs ---- */

    let label = "coeffs" ;
    doff = ninv;
    trash = varmat[1:ninv,1:ninv];
    if not (abs(det(trash)) < 1E-30) ;
        chi2 = parvec'inv(trash)*parvec;
        pval = 100*cdfchic(chi2,doff);
    else;
        "variance matrix singular";
        chi2 = -99;
        pval = -99;
    endif;

    /* ---- joint all b ---- */

    label = label | "all b";

    select = seqb(ninv+1,ninv+nb,1);
    "all b" select';
    trash = varmat[select,select];
    doff = doff|nb;
    if not (abs(det(trash)) < 1E-30);
        chi2  = chi2 | b'inv(trash)*b;
        pval  = pval | 100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]);
    else;
        "variance matrix singular";
        chi2 = chi2 | -99;
        pval = pval | -99;
    endif;

    /* -- b on scaled constants -- */

    if ninst > 0;
        label = label | "insts";
        select = seqb(ninv+2,ninv+1+ninst,1);
        "scaled const" select';
        trash = varmat[select,select];
        doff = doff|rows(select);
        if not (abs(det(trash)) < 1E-30);
            chi2  = chi2 | b[select-ninv]'inv(trash)*b[select-ninv];
            pval  = pval | 100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]);
        else;
            "variance matrix singular";
            chi2 = chi2 | -99;
            pval = pval | -99;
        endif;
    endif;

    /* nonscaled extra factors, extra cons.*/

    if maxc(howto[2] .== 3|7);

        /* ---- investment b ---- */

        label = label | "inv. b";
        select = seqb(ninv+1+1,ninv+1+ninv,1);
        "investment b" select';
        doff = doff | rows(select);
        trash = varmat[select,select];
        if not (abs(det(trash)) < 1E-30);
            chi2 = chi2 | b[select-ninv]'inv(trash)*b[select-ninv];
            pval  = pval | 100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]);
        else;
            "variance matrix singular";
            chi2 = chi2 | -99;
            pval = pval | -99;
        endif;

        /* ---- extra b ---- */

        label = label | "extra b";
        select = seqb(ninv+1+ninv+1,ninv+1+ninv+nxtr,1);
        "extra b" select';
        doff = doff|rows(select);
        trash = varmat[select,select];
        if not (abs(det(trash)) < 1E-30);
           chi2 = chi2 | b[select-ninv]'inv(trash)*b[select-ninv];
            pval  = pval | 100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]);
        else;
            "variance matrix singular";
            chi2 = chi2 | -99;
            pval = pval | -99;
        endif;
    endif;

    /* tests on scaled investment factors */

    if (howto[2] == 4);

        /* unscaled  b */

        label = label | "unsc. b";
        select = seqb(ninv+1+ninst+1,ninv+1+ninst+ninv,1);
        "unscaled b, no const" select';
        trash = varmat[select,select];
        doff = doff|rows(select);
        if not (abs(det(trash)) < 1E-30);
            chi2 = chi2 | b[select-ninv]'inv(trash)*b[select-ninv];
            pval  = pval | 100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]);
        else;
            "variance matrix singular";
            chi2 = chi2 | -99;
            pval = pval | -99;
        endif;

        /* scaled b */

        label = label | "scaled b";
        select = seqb(ninv+1+ninst+ninv+1,ninv+nb,1);
        "scaled b" select';
        doff = doff|rows(select);
        trash = varmat[select,select];
        if not (abs(det(trash)) < 1E-30);
            chi2 = chi2 | b[select-ninv]'inv(trash)*b[select-ninv];
            pval  = pval | 100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]);
        else;
            "variance matrix singular";
            chi2 = chi2 | -99;
            pval = pval | -99;
        endif;

        /* bNR  bR, etc. */

        j1 = 1;
        do while j1 <= cols(invest);
            label = label | invlbl[j1];
            select = seqa(ninv+1+ninst+j1,ninv,ninst+1);
            "nr, r" select';
            doff = doff|rows(select);
            chi2 = chi2 | b[select-ninv]'inv(varmat[select,select])*
                           b[select-ninv];
            pval = pval | 100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]);
            j1 = j1 + 1;
        endo;

    endif;


    if maxc(howto[2] .== (5|8));

        /* ---- investment b ---- */
        ninst = cols(inst);
        label = label | "inv. b";
        doff = doff | (ninv*(ninst+1));
        select = seqb(ninv+1,ninv+ninv,1);
        k = 1;
        do while k <= ninst;
           select = select|seqb(ninv+k*(ninv+nxtr)+1,ninv+k*(ninv+nxtr)+ninv,1);
           k = k+1;
        endo;
        trash = varmat[select,select];
        if not (abs(det(trash)) < 1E-30);
            chi2 = chi2 | b[select-ninv]'inv(trash)*b[select-ninv];
            pval  = pval | 100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]);
        else;
            "variance matrix singular";
            chi2 = chi2 | -99;
            pval = pval | -99;
        endif;


        /* ---- extra b ---- */
        label = label | "extra b";
        if howto[2] == 5;
            doff = doff| (nxtr*(ninst+1));
            select = seqb(ninv+ninv+1,ninv+ninv+nxtr,1);
            k = 1;
            do while k <= ninst;
                select = select|seqb(ninv+k*(ninv+nxtr)+ninv+1,
                                ninv+(k+1)*(ninv+nxtr),1);
                k = k+1;
            endo;
        endif;
        if howto[2] == 8;
            doff = doff| 1;
            select = npar+ninv+1;
        endif;
        trash = varmat[select,select];
        if not (abs(det(trash)) < 1E-30);
            chi2 = chi2 | b[select-ninv]'inv(trash)*b[select-ninv];
            pval  = pval | 100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]);
        else;
            "variance matrix singular";
            chi2 = chi2 | -99;
            pval = pval | -99;
        endif;
    endif;
 endjb:
 return;

/*----------------------------------------------------------- */
/* D. CHISQ: Chi-Square tests of the b's                      */
/*----------------------------------------------------------- */

chisq:
    if not howto[2] == 6;
    " ";
    "Chi-squared tests on coefficients";
    label = blank | label;
    format /RD 12,2;
    $label';
    let label = "chi^2";
    $label;; chi2';
    format /RD 12,0;
    let label = "DF";
    $label;; doff';
    format /RD 12,3;
    let label = "% p val" ;
    $label;; pval';
  endif;

    /*
    "Notes:"
    "mean excess investment returns (%)" 100*meanc(reti-rf)';
    "correlation of estimated paramters";
    varmat./(se*se');
     */
return;




return;

/* ---------------------------------------------------------- */
/* D. OVERID: Construct the overidentifying restrictions test */
/* ---------------------------------------------------------- */

overid:

    /* ------------------------------- */
    /* 1. overidentifying restrictions */
    /* ------------------------------- */

    "";
    label = "Unrestr";
    if not howto[2] == 6;
        doff = cols(allr)-rows(tpar)-rows(b);
    endif;
    if howto[2] == 6;
        doff = cols(allr)-1;
    endif;
    if doff > 0;
        if firststg ; /* use formula for fixed weighting mx. This gives same
                          result for second stage. */
            i = eye(rows(s));
            term = (i-d0*inv(d0'*w*d0)*d0'*w);
            varg = term*s*term';
            {va,ve} = real(eighv(varg));
            npar = rows(tpar)+rows(b);
            if maxc(abs(va[1:npar])) > 1E-6;
                "Warning: large eigenvectors in pseudo-inverse. Check va";
                10000*va';
            endif;
            va[1:npar]=zeros(npar,1);   /* set eig should = 0 to 0 exactly */
            va[npar+1:rows(va)] = 1./va[npar+1:rows(va)];
            lam = diagrv(i,va);
            chi2 = T*lgt'(ve*lam*ve')*lgt;
            pval2 = 100*cdfchic(chi2,doff);
        else;     /* simpler, numerically equivalent 2nd stage formulas */
                  /* included just to reassure that it's the usual */
            chi2 = jt;
            pval2 = 100*cdfchic(chi2,doff);
        endif;
    else;
        doff "degrees of freedom, no overidentifying restrictions tests";
        goto endovid;
    endif;

    /* mean squared pricing errors */

    "RMSE pricing error (not including constant)" ;
    "100 [E(mR)'E(mR)]^(1/2)/E(m)*#rets ";
    preder = meanc(m.*allr[.,2:cols(allr)])/meanc(m);
    100*(preder'preder)^(1/2)/(cols(allr)-1);
     "Mean(m), (1-mean(m))*100";
     meanc(m) (1-meanc(m))*100;


   /* --------------------------------------------------------- */
   /* 2. difference of chisquared tests for which factors count */
   /* --------------------------------------------------------- */

   if maxc(howto[2] .== 3|5|7|8) ;

       /* ------------------ */
       /* a. investment only */
       /* ------------------ */

        jtur = jt;       /* save for later */
        oldhow = howto;  /* save for later */
        if maxc(oldhow[2] .== 3|7);
            howto[2] = 2;    /* just do investment return */
        endif;
        if maxc(oldhow[2] .== 5|8);
            howto[2] = 4;    /* just do scaled investment ret. */
        endif;
    /*  {alpha,delta,mpk,gam,jt,covs,b , f,allr,ik,reti,tpar,junk,lgt3} =  */
        {   T1,   T2, T3, T4,jt,  T5,bT,T7,allrT,T9, T10, tpT, T12,lgt3,m,w} =
            calloptc(alpha,delta,mpk,gam,selpar,
            invest,rassetx,xtrfact,inst,rf,W,howto,b);

        /*  "* Estimates of investment only system using unrestricted w *";
            gosub outpar;
            gosub outse; */

        label = label | "inv.only";
        doff = doff | (cols(allrT)-rows(tpT)-rows(bT));
        chi2 = chi2 | jt;
        pval2 = pval2 | (100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]));

        label = label | "inv-unr";
        if (jt-jtur)<0;
            chi2 = chi2 | 0;
            doff = doff | 0;
        else;
            doff = doff | (doff[rows(doff)]-doff[1]);
            chi2 = chi2 | (jt-jtur);
        endif;
        pval2 = pval2 | (100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]));

        howto = oldhow;                  /* restore howto */

        /* ----------------------- */
        /* b. extra factors  only  */
        /* ----------------------- */

        if maxc(howto[2] .== 3|5);
            let blank = " ";
            cfact = 0;
            if howto[2] == 5;
                cfact = 1;          /* include scaled factors in restr. test? */
            endif;
            {jt,dofobs,junk,pj,ej,lgt5} = obsfact(rassetx~(reti-rf),xtrfact,
                                0,inst,nwlag,0,0,w,blank,blank,cfact,
                                howto[3],0,0,0,0,2);
        endif;
        if maxc(howto[2] .== 7|8);
        /*  {alpha,delta,mpk,gam,jt,covs, b, f,allr,ik,reti,tpar,junk,lgt7} =*/
         {   T1,   T2, T3, T4,jt,  T5,T6,T7,allrT,T8, T9, T10, T11,lgt7,m,w} =
                calloptc(alpha,delta,mpk,-1,0|0|0,
                invest,rassetx,xtrfact,inst,rf,W,howto[1]|9|howto[3],b);
      /* make sure you're not overwriting an interesting w */
            dofobs = cols(allrT)-2;
        endif;

        label = label | "xtraonly";
        doff = doff | dofobs;
        chi2 = chi2 | jt;
        pval2 = pval2 | (100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]));

        label = label | "xtr-unr";
        if (jt-jtur)<0;
            chi2 = chi2 | 0;
            doff = doff | 0;
        else;
            doff = doff | (doff[rows(doff)]-doff[1]);
            chi2 = chi2 | (jt-jtur);
        endif;
        pval2 = pval2 | (100*cdfchic(chi2[rows(chi2)],doff[rows(doff)]));

    endif;

    "Chi squared tests on JT for overidentifying restrictions" ;
    format /RD 10,4;
    label = blank | label;
    $label';
    format /RD 10,2;
    "        JT" chi2';
    format /RD 10,0;
    "        DF" doff';
    format /RD 10,3;
    " % p value" pval2';

    endovid:

return;

endp;