/* -------------------------------------------------------------------------- */
/*                           FINDJTC.G

Saved 2/29 in case new modifications were a disaster
    version of findjt adapted to include consumption model

    Contains findjt and supporting procs.  Findjt takes parameters gives
JT statistic. Since findjt must be a function of only searched parameters,
procs passtof and passback simulate passing other variables by creating globals
reading them and destroying them.

Calls to these procs must look like the following:


[ set up all right hand variables to passtof, in particular
    ik = makeik(delta,invest);
and
    tpar = encode(alpha,delta,mpk,selpar);
b may be 0 if howto[1] /= 2 .   ]

    call passtof(invest,rassetx,xtrfact,rf,ik,W,selpar,howto,alpha,delta,mpk,b);
then
    Jt = findjt(tpar)
or
    call to optmum(&findjt, tpar)
or
    findgt(tpar);
or
    gradwb(tpar);
then
    {reti,ik,m,f,lgt,b,covs,allr} = passback;


Vectors are ordered as follows:

rassetx, etc. :
asset returns~investment returns~scaled asset return1~scaled inv ret. 1,...

factors:
ri1~r12~x1 ri1z1~ri2z1~x1x1 ...

*/
/* -------------------------------------------------------------------------- */

library pgraph,optmum,gmm;


/* -------------------------------------------------------------------------- */
/*                              proc passtof

    Passes globals to findjt.
Usage:
    call passtof(invest,rassetx,xtrfact,inst,
                      rf,ik,W,selpar,howto,alpha,delta,mpk,b);
Inputs:
    invest
        Investment series
    rassetx
        excess returns to test model. adds investment returns automatically
    xtrfact
        extra observable factors. Market, CRR, etc.
    inst
        instruments, 0 if only constant.
    rf
        risk-free rate for making excess returns
    ik
        i/k series.  if selpar[2] /= 1, uses this one rather than recreate
    W
        GMM weighting matrix. 0 if you want I, will create right size
    selpar
        3x1 vector, tells you which of alpha delta mpk t9o search over.
    howto
        howto[1] = 1: estimate b
                   2: use given b, hold constant for gradients, std. errors
        howto[2] = 1: factors = excess investment returns
                   2: factors = investment returns (not much diff)
                   3: factors = investment returns + xtrfact
                   4: factors = scaled investment returs
                   5: factors = scaled investment + extrafact
                   6: factors = consumption gr ^ gamma
                   7: factors = inv. return + cons gr ^ gamma
                   8: factors = scaled inv. return + cons gr ^ gamma
                      (consumption growth is passed in xtrfact)
                   9: factors = cons.gr ^ gamma, price i ret too (used for
                       chi2 difference tests only )
        howto[3] = 0: use unconditional moments, even if scaled factors
                   1: use scaled returns
    alpha delta mpk

    b
        if howto[1] = 2, uses this one rather than estimate.

*/
/* -------------------------------------------------------------------------- */

proc(0) =
    passtof(invest,rassetx,xtrfact,inst,rf,ik,W,selpar,howto,alpha,delta,mpk,b);

     clearg _invest,_rassetx,_xtrfact, _inst,
            _rf,_ik,_w,_selpar,_howto,_alpha,_delta,_mpk,_b;
    _invest = invest;
    _rassetx = rassetx;
    _xtrfact = xtrfact;
    _inst = inst;
    _rf = rf;
    _ik = ik;
    _w = w;
    _selpar = selpar;
    _howto = howto;
    _alpha = alpha;
    _delta = delta;
    _mpk = mpk;
    _b = b;
endp;


/* -------------------------------------------------------------------------- */
/*                          Proc findjt

Takes Parameters + globals defined above -> T*JT, globals in passback.
Designed for use with optmum.

Usage:
    Jt = findjt(tparms)

Notes:
    _allr gives all excess returns to be tested.  This is
        rassetx|inv.ret|(rassetx x z1)|(inv. ret x z1)| ....

*/
/* -------------------------------------------------------------------------- */


proc(1) = findjtc(tparms);

    local ccode,rix,erix,sigin,f,JT,T,j,trash;
    clearg _reti,_m,_f,_lgt,_covs,_allr,_factors,_gam,_price;

    histry = histry|(tparms'~-99);


    T = rows(_rassetx);
    ccode = 0;
    if maxc(_howto[2] .== 6|7|8|9);
        ccode = 1;
    endif;

    /* ------------------------------ */
    /* form investment returns, etc.  */
    /* ------------------------------ */

    {_alpha,_delta,_mpk,_gam} =decodec(tparms,_selpar,ccode,_alpha,_delta,_mpk,
                                            _gam);

    if _howto[2] /= 6;             /* don't do Ri if only consumption model */
        if _selpar[2] == 1;        /* if choosing delta, form new i/k series */
            _ik = makeik(_delta,_invest);
        endif;

        {_reti,_ik} =  invret(_alpha,_delta,_mpk,_ik);
        rix = _reti - _rf;
        Erix = meanc(rix);
    endif;

    /* ------------------------- */
    /* set up factors            */
    /* NEW: no demeaning factors */
    /* Make sure to echo any changes to gradwb! */
    /* ------------------------- */

    if _howto[2] == 1;                     /* excess investment returns */
        _factors = rix;
    endif;
    if _howto[2] == 2;                     /* investment returns */
        _factors = _reti;
    endif;
    if _howto[2] == 3;                     /* comparison tests w. other facts*/
        _factors = _reti~_xtrfact;
    endif;
    if _howto[2] == 4;                     /* scaled investment returns */
        _factors = _reti;
        j = 1;
        do while j <= cols(_inst);
             _factors = _factors~(_reti.*_inst[.,j]); /* No scaling 1 by z */
        @    _factors = _factors~( (ones(rows(_reti),1)~_reti).*_inst[.,j]);@
                                                       /* also scale 1 by z */
            j = j+1;
        endo;
    endif;
    if _howto[2] == 5;
        _factors = (_reti~_xtrfact);      /* scaled and comparison */
        j = 1;
        do while j <= cols(_inst);        /* no scaling 1 by z!!! */
            _factors  = _factors ~ ( (_reti~_xtrfact).*_inst[.,j]);
            j = j+1;
        endo;
    endif;

  /*   if maxc(_howto[2] .== 6|9;)    m = cgr^gam directly below  endif;   */

    if _howto[2] == 7;
        _factors = _reti~(_xtrfact^_gam);
    endif;

    if _howto[2] == 8;
        _factors = _reti;
        j = 1;                                 /* is it ok to demean cons? */
        do while j <= cols(_inst);
        /*_factors = _factors~((_reti.*_inst[.,j])-meanc(_reti.*_inst[.,j])');*/
            _factors = _factors~(_reti.*_inst[.,j]);
            j = j+1;
        endo;
        _factors = _factors~(_xtrfact^_gam);
        /* _factors = _factors~(_xtrfact^_gam - meanc(_xtrfact^_gam));  */
    endif;

    if (levels[1]==1) or (levels[1]==2);
        _factors = _rassetx[.,1]~_factors;
        @_factors = ones(rows(_factors),1)~_factors;@ /* adds constant */
    endif;
    if (levels[1] == 0);
        _factors = _factors-meanc(_factors);
    endif;

    /* ------------------------------------------ */
    /* set up vector of asset returns to test     */
    /* _allr is  asset~investment .*. instruments */
    /* ------------------------------------------ */

    if _howto[2] == 6;
        trash = _rassetx;
    else;
        if levels[1]==1;
            trash = _rassetx~_reti;
            _price = ones(cols(trash),1);
        elseif levels[1] == 2 ;
            trash = _rassetx~rix;
            _price = 1|zeros(cols(trash)-1,1);
        elseif levels[1] == 0 ;
            trash = _rassetx~rix;
        endif;
    endif;

    _allr = trash;

    if _howto[3] == 1;
        j = 1;
        do while j <= cols(_inst);
          if levels[2] == 0;      /* scale rf by inst too */
            _allr = _allr ~ ((trash).*_inst[.,j]);
            if levels[1] == 1;
                _price = _price|ones(cols(trash),1);
            elseif levels[1] == 2;
                _price = _price|meanc(inst[.,j])|zeros(cols(trash)-1,1);
            endif;
          endif;
          if levels[2] == 1;      /* do not scale rf by inst */
            _allr = _allr ~ ((trash[.,2:cols(trash)]).*_inst[.,j]);
            if levels[1]==1;
                _price = _price|ones(cols(trash)-1,1);
            elseif levels[1] == 2;
                _price = _price|zeros(cols(trash)-1,1);
            endif;
          endif;
          j = j+1;
        endo;
    endif;

     /* test for overid: fix this up and get a return code */

     if cols(_allr) < cols(_factors)+cols(_invest);
        "System underidentified, will bomb";
     endif;


     /* ------------------------------ */
     /* set up W, covs, find b, do GMM */
     /* ------------------------------ */

    if _W == 0;
        _W = eye(cols(_allr));
    endif;

    if not maxc(_howto[2] .== 6|9);
        _covs = 1/T*(_factors)'(_allr);
        if levels[1]==0;
            _covs = _covs - meanc(_factors)*(meanc(_allr)');
        endif;

        clearg b0;
        b0 = 1;
        if _howto[1] /= 2;            /* if _howto[1] == 2, use given b */

/*
          if levels[1] == 1;
            _b =  inv(_covs*_W*_covs')*_covs*_W*_price;
          elseif levels[1]==2;
            _b =  inv(_covs*_W*_covs')*_covs*_W*_price;
          elseif levels[1]==0;
            _b = - inv(_covs*_W*_covs')*_covs*_W*meanc(_allr)*b0;
           endif;
*/
          clearg _cw,_ctil,_ptil;
          local q,r;

          _cw = chol(_w);                  /*  c'c=w */
          _ctil = _cw*_covs';
          _ptil = _cw*_price ;
          if (levels[1] == 1) or (levels[1] == 2);
            {q,r} = qqr(_ctil);
            _b = inv(R)*Q'_ptil;
            /* exactly same as _b = olsqr(_ptil,_ctil); */
          endif;
        endif;

        if (levels[1]==1) or (levels[1]==2);
            _m = _factors*_b;
        elseif levels[1]==0;
            @_m = b0+ (_factors)*_b;@
            _m = b0+ (_factors-meanc(_factors)')*_b;
        endif;
    endif;
    if maxc(_howto[2] .== 6|9);
        _m = _xtrfact^(_gam);
    endif;

     _f   = (_m .* (_allr));
     if levels[1]==1;
        _lgt = meanc(_f) - 1;
     elseif levels[1]==2;
        _lgt = meanc(_f) - _price;
     else;
        _lgT = meanc(_f);
     endif;
    JT   = T*_lgT'_W*_lgT;

    histry[rows(histry),cols(histry)]= jt;

    retp(JT);
endp;

/* ---------------------------------------------------------------------------*/
/*                                 Proc findgtc

    Calls findjt, returns gt.  Use to take grad of gt with respect to pars
for GMM standard errrors.

*/
/* ---------------------------------------------------------------------------*/

proc(1) = findgtc(tpars);
    call findjtc(tpars);
    retp(_lgt);
endp;


/* ---------------------------------------------------------------------------*/
/*                          new proc gradwbc
   does analytical gradients of Jt wrt mpk's.

   Merge with findjtc!  gradwb then just sets a switch and calls.

*/
/* ---------------------------------------------------------------------------*/

proc(1) = gradwbc(tpars);

    local jt, ccode,drdmp,projmat,djt,i,drdmpi,dfdmpi,j,trash,dradmpi,
    dcdmpi,dctil;

    Jt = findjtc(tpars) ; /* included to make sure globals ok */

    if maxc(_howto[2] .== 6|7|8|9);
        ccode = 1; "Error: gradwbc not implemented for consumption model ";
    endif;

    /* find d RIi / d mpk i  */

    if _selpar[3] == 1;
        drdmp = (1-_delta')./(1 + _alpha'.* _ik[1:rows(_reti),.]);
    else;
        "Error: Gradwbc only implemented for mpk. ";
    endif;

    /*projmat = _ctil*inv(_ctil'_ctil)*_ctil';  */     /* compute this once */
    local q,r;
    {q,r} = qqr(_ctil);
    projmat = Q*Q';
    projmat = eye(rows(projmat))-projmat;

    dJt = zeros(1,rows(tpars));                /* will hold answer */

    /* have to do this item by item, since formula for dJT/dpar doesn't
       easily vectorize */

    i = 1;
    do while i <= rows(tpars);

        /* ------------------------- */
        /* set up df/dmpki           */
        /* ------------------------- */

        drdmpi = zeros(rows(drdmp),cols(drdmp));
        drdmpi[.,i] = drdmp[.,i];

        if _howto[2] == 1;                     /* excess investment returns */
            "Error: gradwb only implemented for howto = 2,4";
        endif;
        if _howto[2] == 2;                     /* investment returns */
            dfdmpi = drdmpi;
        endif;
        if _howto[2] == 3;                /* comparison tests w. other facts*/
            dfdmpi = drdmpi~zeros(rows(_xtrfact),cols(_xtrfact));
        endif;
        if _howto[2] == 4;                     /* scaled investment returns */
            dfdmpi = drdmpi;
            j = 1;
            do while j <= cols(_inst);
                dfdmpi = dfdmpi~( (ones(rows(_reti),1)~drdmpi).*_inst[.,j]);
                                                         /* also scale 1 by z */
                j = j+1;
            endo;
        endif;
        if _howto[2] >= 5;
             " Gradwb: howto >= 5 not yet implemented";
        endif;

        if (levels[1]==1) or (levels[1]==2);
            dfdmpi = zeros(rows(dfdmpi),1)~dfdmpi;
        endif;
        if (levels[1] == 0);
            "gradwp not implemented for levels[1] = 0 ";
        endif;

        /* ------------------------------------------ */
        /* set up vector of drasst/dmpi                  */
        /* ------------------------------------------ */

        trash = zeros(rows(_rassetx),cols(_rassetx));

        if _howto[2] == 6;
        else;
            if levels[1]==1;
                trash = trash~drdmpi;
            elseif levels[1] == 2 ;
                trash = trash~drdmpi;
            elseif levels[1] == 0 ;
                trash = trash~drdmpi;
            endif;
        endif;

        dradmpi = trash;

        if _howto[3] == 1;
            j = 1;
            do while j <= cols(_inst);
              if levels[2] == 0;      /* scale rf by inst too */
                dradmpi = dradmpi ~ ((trash).*_inst[.,j]);
              endif;
              if levels[2] == 1;      /* do not scale rf by inst */
                dradmpi = dradmpi ~ ((trash[.,2:cols(trash)]).*_inst[.,j]);
              endif;
              j = j+1;
            endo;
        endif;

        /* calculate c/theta[i] */

       if not maxc(_howto[2] .== 6|9);
             dcdmpi = 1/rows(dradmpi)*
                       (dradmpi'_factors + _allr'dfdmpi);   /* d/dmpi E Rf'*/
       endif;

       dctil = _cw*dcdmpi;

       djt[i] = -2*rows(_allr)*_ptil'projmat*dctil*_b;

       /* "FIx following that olsqr only accepts vector y */
       /* djt[i] = -2*rows(_allr)*_ptil'(dctil - _ctil*olsqr(dctil,_ctil))*_b;*/


       i = i+1;
    endo;

    retp(djt);
endp;



/* ---------------------------------------------------------------------------*/

/* ---------------------------------------------------------------------------*/
/*               old    proc gradwb

  Analitycal grad of Jt wrt a,d,mpk, using optimal b's. use to speed up searches
using findjt.  Uses globals from passto and defined in a first
call to findjt. clear when done with a call to  passback.
DO NOT use when searching over delta, since it does notcalculate change in ik
when you recreate series.
Note, though it seems
to be right, there are significant differences with results of gradp findjt.
Why?

example call and comparison to gradp:

    alpha = 3|3;          /* original starting values */
    delta = .1|.1;
    mpk =   .2|.2;
    invest  = gin82~gir82;
    W = eye(cols(rassetx)+cols(inv));
    ik = makeik(delta,invest);
    selpar = 1|0|0;
    tpar = encode(alpha,delta,mpk,selpar);
    call passtof(invest,rassetx,xtrfact,rf,ik,W,selpar,howto,alpha,delta,mpk,b);
    format /RD 12,5; gradwb(tpar);
    call passback ;
    call passtof(invest,rassetx,xtrfact,rf,ik,W,selpar,howto,alpha,delta,mpk,b);
    gradp(&findjt,tpar);

this proc needs to be cleaned!
*/
/* ---------------------------------------------------------------------------*/
/*
proc(1) = gradwbc(tpars);
    local c,tik,drda,er,re,ca,ca1,rix,np,nf,nre,j,nxtr,
            edrda1,edrda2,dreda,dreda1,drida,drida1,
            invtrm,xa,ga,djt,counter,g,gterm,Tl,edrda,edrterm,k;

    call findjtc(tpars);
    C = _covs';                       /* C = E(Rf') */
    Er = meanc(_allr);
    Re = (_allr)-Er';
    invtrm = inv(C'_W*C);

    Tl = rows(_reti);
    nf = cols(C);          /* number of factors */
    np = cols(_invest);    /* number of investment returns */
    nre = cols(Re);        /* number of assets, incl scaled */
    nxtr = 0;
    if (_howto[2] == 3) or (_howto[2] == 5);
        nxtr = cols(_xtrfact);
    endif;                 /* number of extra factors, not incl scaled */


    gterm = (eye(rows(C)) - C*invtrm*C'_W);
    g = gterm*Er;

    djt = zeros(sumc(_selpar)*rows(_alpha),1);
                            /* this is where the answer goes */

    counter = 1;

    /* ---------------------------------------------------------*/
    /* calculate dRI/dparameter.  note this is a matrix:        */
    /* drda =  dR1/da1 ~ dR2/da2 ....                           */
    /* ---------------------------------------------------------*/

    if _selpar[1] == 1;
        drda = (1-_delta') .* (_ik[2:Tl+1,.] + 1/2*_ik[2:Tl+1,.]^2) ./
                            (1 + _alpha' .* _ik[1:Tl,.]) -
                _reti .* _ik[1:Tl,.] ./ (1 + _alpha' .* _ik[1:Tl,.] );
        gosub update;
    endif;

    if _selpar[2] == 1;           /* warning: does not include effect on ik */
        "gradwb: Using analytical grad on delta, does not include ik effect";
        drda = _reti./(1-_delta');
        gosub update;
    endif;

    if _selpar[3] == 1;
        drda = (1-_delta')./(1 + _alpha'.* _ik[1:Tl,.]);
        gosub update;
    endif;


    retp(Tl*djt');


update:

    j = 1;                                  /* loops through parameters   */
    do while j <= np;

        /* ---------------------------------- */
        /* form dRIda =                       */
        /*  d(RI1 RI2 Rx zRI1 zRI2 zRx ..)/da */
        /* ---------------------------------- */

        dRIda = zeros(Tl,np+nxtr);
        dRIda[.,j] = drda[.,j];
        if (_howto[2] == 4) or (_howto[2] == 5);
            drIda1 = drIda;
            k = 1;
            do while k <= cols(_inst);
                drIda1 = dRIda1~(_inst[.,k].*dRIda);
                k = k+1;
            endo;
            drIda = drida1;
        endif;

        /* ---------------------------- */
        /* form dReda =                 */
        /* d(Ra RI1 RI2 zRa zRI1 zRI2..)*/
        /*    /da                       */
        /* ---------------------------- */

        dreda = zeros(Tl,cols(_rassetx)+cols(_invest));
        dReda[.,cols(_rassetx)+j]=drda[.,j];
        if _howto[3] == 1;
            dreda1 = dreda;
            k = 1;
            do while k <= cols(_inst);
                dreda1 = dreda1~(dreda.*_inst[.,k]);
                k = k+1;
            endo;
            dreda = dreda1;
        endif;

        /* ------- */
        /* form Ca */
        /* ------- */

        Ca = 1/Tl*((Re'drida)+(dreda'_factors));

        /* -------------------------------------------- */
        /* form E(dR/da) =                              */
        /* (0 for assets)|(0 0 dRI/da 0 0 )|            */
        /* (0 for scaled assets)|(0 0 zdRIda 0 0 )| ... */
        /* -------------------------------------------- */

        edrda1 = meanc(drda);
        edrda = 0*edrda1;
        edrda[j] = edrda1[j];
        edrda = zeros(cols(_rassetx),1)|edrda;

        if _howto[3] == 1;
            k = 1;
            do while k <= cols(_inst);
                edrda2 = meanc(_inst[.,k].*drda);
                edrda1 = 0*edrda2;
                edrda1[j] = edrda2[j];
                edrda = edrda|(zeros(cols(_rassetx),1)|edrda1);
                k = k+1;
            endo;
        endif;

        Xa  = (eye(rows(C)) - C*invtrm*C'_W )*Ca*invtrm*C';
        ga = - (Xa+Xa') * _W * Er;

        ga = ga + gterm*edrda;

        dJt[counter] = (2*ga'_w*g);
        counter = counter+1;
        j = j+1;
     endo;

return;
endp;
*/
/* -------------------------------------------------------------------------- */
/*                            proc passback

    Returns globals defined by findjt, then clears them so you don't get tempted
to use them.

Usage:
    {reti,ik,m,f,lgt,b,covs,allr} = passback;

Returns:
    reti = investment return
    ik   = inv/capital ratio
    m    = m
    f    = m*asset return~investment return
    lgt  = E(f), Hansens' gT
    b    = coefficients in m = 1 + (f-Ef)'b
    covs = cov(F,ex. asset return ~ ex invest. return)
    allr = rassetx | (reti-rf) | rassetx .* z1 |  (reti-rf) .* z1 | ....
*/
/* -------------------------------------------------------------------------- */


proc(9) = passbakc;          /* globals cleared & defined by findjt above */
    local reti,ik,m,f,lgt,covs,b,allr,neg_m;
    reti = _reti;
    ik   = _ik;
    neg_m = counts(_m,0)/rows(_m);
    m    = _m;
    f    = _f;
    lgt  = _lgt;
    b    = _b;
    covs = _covs;
    allr = _allr;
                                /* done to enforce no other procs use these */
    clearg _reti,_m,_f,_lgt,_b,_covs,_factors;
    clearg _invest,_rassetx,_xtrfact,_inst,
                 _rf,_ik,_w,_selpar,_howto,_alpha,_delta,_mpk,_gam;
    retp(reti,ik,m,f,lgt,b,covs,allr,neg_m);
endp;