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

    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

                   /* eliminated in favor of howto[4] */
             /*      4: factors = scaled investment returs   */
             /*      5: factors = scaled investment + extrafact  */

                   /* consumption options. c = xtrafact */
                   6: factors = beta consumption gr ^ gamma
                   7: factors = inv. return + beta cons gr ^ gamma
                   8: factors = scaled inv. return + beta cons gr ^ gamma

                  10: use xtrafact only (for inv. growth, capm, crr models )

        howto[3] = 0: use unconditional moments (even if scaled factors)
                   1: use scaled returns

        howto[4] = 0: use unscaled factors
                   1: use scaled factors
                   2: scale investment returns but not extra factors

        howto[5] = 0: price investment returns too
                   1: do not price investment returns

        howto[6] = 0: do not scale constant by instrument
                   1: do scale constant by instrument

        howto[7] = method.
                /*    1: Include E(m) or e(mrf) as a moment, factor levels    */
                /*          2: Old way. demean factors, m=1+b'(f-Ef)          */
                /*          3: prespecify b0 = 1                              */
                /*          4: just like 1; calculated by trans to cov Rf     */

    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,select);

     clearg _invest,_rassetx,_xtrfact, _inst,
            _rf,_ik,_w,_selpar,_howto,_alpha,_delta,_mpk,_b,_select;
    _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;
    _select = select;
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,constf,method;
    local w11,w12,w21,w22,fs,rs,rbar,covs,invtrm,b0str,b1,b0;
    clearg _reti,_m,_f,_lgt,_covs,_allr,_factors,_gam,_price,_b0;

    _b0 = 1;        /* for use with method 3; m = b0 + b'f           */
    method = _howto[7];

    if 0;
        format /RD 10,6; "findjtc: parameters" tparms';
    endif;

    T = rows(_rassetx);

    ccode = 0;
    if maxc(_howto[2] .== 6|7|8);
        ccode = 1;
    endif;
    if maxc(_howto[2] .== 7|8) and maxc(maxc(_select .== 4)) ;
        ccode = 0;                             /* no c if doing chi2 diff test*/
    endif;

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

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

    if (_howto[2] /= 6) and (not (_howto[2] == 10 and _howto[5] == 1));
                                 /* don't do Ri if only consumption model */
                                 /* or extra only, not pricing Ri */
        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 you use it            */
    /* ------------------------- */

    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 maxc(_howto[2] .== 4|5);
        "Findjtc.g: howto[2] = " _howto[2] " This value is no longer used";
    endif;

    if _howto[2] == 6;                      /* consumption alone  */
        _factors = _xtrfact^_gam;
    endif;

    if (_howto[2] == 7) or (_howto[2] == 8);  /* ri + consumption ^gamma */
       _factors = _reti;
    endif;

    if _howto[2] == 10;
        _factors = _xtrfact;
    endif;

    if (_howto[4] > 0) or (_howto[2] == 8); /* add scaled factors if required */
        trash = _factors;
        if _howto[4] == 2;
            trash = _reti;      /* in this case only scale ri by insts */
        endif;
        j = 1;
        do while j <= cols(_inst);
             _factors = _factors~(trash.*_inst[.,j]);
            j = j+1;
        endo;
    endif;

    if (_howto[2] == 7) ;                          /* add consumption to  */
        _factors = _factors~(_xtrfact^_gam);       /* ri + cons model */
    endif;
    if (_howto[2] == 8);
        _factors =
            _factors[.,1:2]~(_xtrfact^_gam)~_factors[.,3:cols(_factors)];
    endif;


/* ------------------------------------- */
/* add constant, possibly demean factors */
/* ------------------------------------- */

constf = ones(rows(_factors),1);
if (_howto[6] == 1) and (_howto[4] == 1);
                                                /* if your'e scaling factors*/
                                                  /* z times 1 as factors */
    j = 1;
    do while j <= cols(_inst);                    /* this produces lower jts, */
        constf = constf~_inst[.,j];               /* but are scaled consts id?*/
        j = j+1;
    endo;
endif;

if ((method == 1) or (method == 4)) and (_howto[2] /= 6 );
    _factors = constf~_factors;
endif;
if method == 2;
   _factors = _factors-meanc(_factors)';
endif;

/* methods not mentioned leave factors as is */

if maxc(maxc(_select)) /= 0;                          /* delete some factors */
    if rows(_select) /= cols(_factors);                /* if any left */
        trash = zeros(cols(_factors),1);
        j = 1;
        do while j <= rows(trash);
            if maxc(maxc(_select.==j));
                trash[j] = 1;
            endif;
            j = j+1;
        endo;
        _factors = delif(_factors',trash)';
    else;
       _m = zeros(rows(_factors),1);  /* if all b = 0, incl constant, then m=0*/
       goto mdefined;
    endif;
endif;

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

    if (_howto[2] == 6) or (_howto[5] == 1) ;
        trash = _rassetx;
    elseif _howto[5] == 0;
        trash = _rassetx~rix;
    endif;

    _allr  =  trash;
    if (_howto[3] == 1);                          /* add scaled returns */
        j = 1;
        do while j <= cols(_inst);
            _allr = _allr ~ ((trash).*_inst[.,j]);
          j = j+1;
        endo;
    endif;


    if (method == 1) or (method == 4 );
        _price = 1|zeros(cols(_allr),1);   /* for rf or 1 in first asset */
@        _allr = ones(rows(_rf),1)~_allr;@   /* with R=1, sets 1 = E(m) */
        _allr = _rf~_allr;                 /* option to use tbill rate */
    endif;
    if method == 2;
        _price = zeros(cols(_allr),1);    /* for old way */
    endif;
    if (method == 3) ;
        _price = zeros(cols(_allr),1);
    endif;

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

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

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

/*       Here's what we're trying to achieve, with QR because more stable
         price not zero version
            _b =  inv(_covs*_W*_covs')*_covs*_W*_price;
         price zero versiom
            _b = - inv(_covs*_W*_covs')*_covs*_W*meanc(_allr)*b0;
*/
        _covs = 1/T*(_factors)'(_allr);

        clearg _cw,_ctil,_ptil;
        local q,r;

        if method == 1;

        /* standard way: */

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

        if (method == 2) or (method == 3);

            /* version for 0 price throughout, old way */

            _cw = chol(_w);                  /*  c'c=w */
            _ctil = _cw*_covs';
            _ptil = _cw*meanc(_allr) ;
            {q,r} = qqr(_ctil);
            _b = -inv(R)*Q'_ptil;

            /* exactly same as _b = olsqr(_ptil,_ctil); */
            if method == 3;
                _b = _b*_b0;
            endif;
        endif;

        if method == 4;          /* same as method 1, calculated w covar */
            w11 = _w[1,1];
            w12 = _w[1,2:cols(_w)];
            w21 = _w[2:rows(_w),1];
            w22 = _w[2:rows(_w), 2:cols(_w)];

            if cols(_factors) > 1;
                fs     = _factors[.,2:cols(_factors)];
                rs     = _allr[.,2:cols(_allr)];
                rbar   = meanc(rs);
                covs   = 1/T*(fs-meanc(fs)')'(rs - rbar');
                invtrm = inv(covs*w22*covs');
                b0str  = w11 + rbar'w21 -
                        (w12 + rbar'w22)*covs'*invtrm*covs*w21;
                b0str  = b0str/(w11 + w12*rbar + rbar'w21 + rbar'w22*rbar -
                      (w12 + rbar'w22)*covs'*invtrm*covs*(w22*rbar + w21) );
                b1 = invtrm*covs*(w21 - (w22*rbar + w21)*b0str);
                b0 = b0str - meanc(fs)'b1;
                _b = b0|b1;
                _covs = 1/T*(_factors)'(_allr);
                                                 /* included so that   */
                                                 /*  D0 can be calculated */
            endif;
            if cols(_factors) == 1;          /* assume 1 is only factor */
                if not (_factors == ones(rows(_factors),1));
                    "Findjtc: error. assumed factors = 1, they aren't ";
                endif;

                _covs = 1/T*(_factors)'(_allr);
                _b =  inv(_covs*_W*_covs')*_covs*_W*_price;
                 /* covs = ER'. b = (ER'WER)^-1 ER'WP; -1 is scalar */
            endif;
        endif;
    endif;

    if (method == 1) or (method == 4) ;
        _m = _factors*_b;
    endif;
    if method == 2;
        _m = 1 + _factors*_b;
    endif;
    if method == 3;
        _m = _b0 + _factors*_b;
    endif;

    mdefined:

    _f   = _m .* _allr;
    _lgt = meanc(_f) - _price;
    JT   = T*_lgT'_W*_lgT;

    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;

/* -------------------------------------------------------------------------- */
/*                            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(10) = passbakc;          /* globals cleared & defined by findjt above */
    local reti,ik,m,f,lgt,covs,b,allr,w,price;
    reti = _reti;
    ik   = _ik;
    m    = _m;
    f    = _f;
    lgt  = _lgt;
    b    = _b;
    covs = _covs;
    allr = _allr;
    w =_w;
    price = _price;
                                /* done to enforce no other procs use these */

    clearg _reti,_m,_f,_lgt,_b,_covs,_factors,_allr;
    clearg _invest,_rassetx,_xtrfact,_inst,_price,
                 _rf,_ik,_w,_selpar,_howto,_alpha,_delta,_mpk,_gam,_select;

    retp(reti,ik,m,f,lgt,b,covs,allr,w,price);
endp;