/* ---------------------------------------------------------------------------*/
/*                          new proc gradwbc
   does analytical gradients of Jt wrt mpk's.
   this used to be appended to findjtc.g.
   I found numerical gradients no slower than analytical since system has
   only two variables to search over. Analytical also no more numerically
   stable than numerical: this suffers the same matrix inversion problems.


   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;

        dfdmpi = zeros(rows(dfdmpi),1)~dfdmpi;

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

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

        if _howto[2] == 6;
        else;
           trash = trash~drdmpi;
        endif;

        dradmpi = trash;

        if _howto[3] == 1;
            j = 1;
            do while j <= cols(_inst);
              dradmpi = dradmpi ~ ((trash[.,2:cols(trash)]).*_inst[.,j]);
              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;