/*                               MONTEC.PRG

Does monte carlo of linear factor pricing models, gmm tests. TO answer
questions about distribution of b, Jt, and spectral density estimates

*/

/*
Regression model for cross section of asset excess returns

  R_t = a + bet f_t + e_t   E(e_i f) = 0; E(e_t) = 0.

If correctly priced, E(R) = bet * lam. Thus,

  bet_i  lam = a_i + bet_i E(f)

hence,
  a_i = bet_i ( lam - E(f) ).

If f is priced (CAPM) then, lam = E(f), a_i = 0. If not, not nec. true.
For example, if f is demeaned, then E(f) = 0, a_i = bet_i lam = E(R_i).
*/

library pgraph, gmm;
Sampl = 1000;                                /* sample size */
maxit = 1;                                  /* nubmer of monte carlo draws*/
rndseed 1234;
verbose = 1;                                /* give stats on each run */

Ef = 0.07;                                  /* mean of factors */
sigf = 0.15;                                /* factor vcv matrix */
nf = ones(rows(ef),1);

lam = Ef;
betas = 0.5|0.75|1.0|1.25|1.5;
nr = rows(betas);
sigma = 0.20*eye(rows(betas));              /* residual vcv matrix */
alphas = betas*(lam-Ef);

"assumed alphas and betas";
alphas~betas;

b1t = -inv(sigf)*lam;                      /* true b's */
b0t = 1-b1t'Ef;
"Assumed b " b0t~b1t';

draw = 1;
do while draw <= maxit;
    f = ef' + rndn(sampl,nf)*chol(sigf);            /* draw factors */
    r = alphas'.*ones(sampl,nr) +
           f*betas' + rndn(sampl,nr)*chol(sigma);   /* draw returns */
        /* check:
                sigma = (1~2)|(2~5);
                test = rndn(1000,2)*(chol(sigma));
                (test'test)/1000 ;
        */
    if verbose;
        olsbeta = r/(ones(sampl,1)~f);
        format /RD 10,4;
        "OLS  betas" olsbeta';
     endif;


    W = eye(nr+1);                              /* start with id matrix */

    w11 = w[1,1];
    w12 = w[1,2:cols(w)];
    w21 = w[2:rows(w),1];
    w22 = w[2:rows(w), 2:cols(w)];

    rbar = meanc(r);
    covs = 1/sampl*(f-meanc(f)')'(r - 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(f)'b1;
    b = b0|b1;

    if verbose;
        "b " b';
    endif;

    bigcovs = (1~meanc(r)')|(meanc(f)~covs);


    draw = draw+1;
endo;