/*    Jim Stock Monte carlo

      My additions are contained in ------------- comments
*/

/* --------------------*/
library pgraph;
graphset;

smode   = 0;
smode2  = 0;
jcmode1 = 1;
replict = 0;     /* use actual errors to reproduce data */
graphs  = 0;

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

nst = 22;         /* number of startup values; must >= klags */
nmc  = 1;
n = 60;
rhot = 1;
seed = 2691;
klags = 2;

fboot = 1;       /*  bootstrap */
fmonte = 0;      /* monte carlo */

load ehat = cerr_1;
ehat = ehat[.,2]~ehat[.,1];   /* file has d~p */
load lddat;
load lpdat;
load truebcr;     /*  these are generated by newmonte using replict */
load truebyr;
load truebcf;     /*  these are generated by newmonte using force=1 */
load truebyf;

tcrit= 2.11;
rcrit = .225;
output file = cg1_3.out reset;

/* --------------------------------------------- */
if smode; "Stock program";           endif;
if jcmode1  ; "with estimated constants"; endif;
/* --------------------------------------------- */

format /RDN 10,3;
"                 NMC = " nmc;
"            n = nobs = " n;
"                seed = " seed;
"klags = no. var lags = " klags;
"";

let sig[2,2]=
    379.66   184.96
    184.96   205.05 ;

if fboot; "Computed by bootstrap";
else;
    "Monte Carlo using covariance matrix";
    sig;
    sigchol = chol(sig);
endif;

rho = 1;

    ntot=n+nst;
    rndseed(seed);
    imc=1;
    do until imc > nmc;
        if fboot;
            indx9 = trunc((rows(ehat)*rndu(ntot-2,1))+1);
            e = zeros(2,2)|ehat[indx9,.];
        endif;
        if fmonte;
            e = rndn(ntot,2);
            e = e*sigchol;
        endif;
        if replict;
            e = zeros(nst,2)|ehat;
        endif;

        dlp = zeros(rows(e),1);
        dld = dlp;
        ldp = dlp;

/* ----------------------------------------------- */
/* jc starts with first two observations from data */

        if jcmode1 or replict;
            dld[1] =  15.619237;
            dlp[1] =  28.743327;
            dld[2] =  12.536408;
            dlp[2] =  32.903990;
            ldp[1|2] = -300.708086 | -321.075667;
         endif;

/* jc has constants in regression  */
        if smode;
            constd = 0;
            constp = 0;
        endif;
        if jcmode1;
            constd =  truebcf[1];
            constp =  truebyf[1];
        endif;
/* ----------------------------------------------- */

        if smode2;
            dlp = e[.,1];
            ldp = recserar(e[.,2],e[1,2],rho);
            ld = ldp+cumsumc(dlp);
            dld = 0|(ld[2:rows(ld),1]-ld[1:(rows(ld)-1),1]);
            /* what's going on here ? */
        endif;

        if smode or jcmode1;
            i9 = 3;                         /* stock has 4 */
            do until i9 > rows(dlp);
                dlp[i9,1] =
                           constp +
                           truebyf[5]     * dlp[i9-1,1] +
                           truebyf[6]     * dlp[i9-2,1] +
                           truebyf[3]     * dld[i9-1,1] +
                           truebyf[4]     * dld[i9-2,1]  + e[i9,1];

                dld[i9,1] =
                            constd +
                           truebcf[5]      * dlp[i9-1,1] +
                           truebcf[6]      * dlp[i9-2,1] +
                           truebcf[3]      * dld[i9-1,1] +
                           truebcf[4]      * dld[i9-2,1]  + e[i9,2];

                i9 = i9+1;
            endo;
/* stock:
            ld = cumsumc(dld);
            lp = cumsumc(dlp);
            ldp = ld-lp;
*/
/* me: initialize to level*/
            ldp = recserar(dld-dlp,ldp[1],1);

        endif;

        if replict;
            i9 = 3;                         /* stock has 4 */

            do until i9 > rows(dlp);
                dlp[i9,1] =
                           truebyr[1] +
                           truebyr[2]      * ldp[i9-1,1] +
                           truebyr[5]      * dlp[i9-1,1] +
                           truebyr[6]      * dlp[i9-2,1] +
                           truebyr[3]      * dld[i9-1,1] +
                           truebyr[4]      * dld[i9-2,1]  + e[i9,1];

                dld[i9,1] =
                            truebcr[1] +
                            truebcr[2]      * ldp[i9-1,1] +
                            truebcr[5]      * dlp[i9-1,1] +
                            truebcr[6]      * dlp[i9-2,1] +
                            truebcr[3]      * dld[i9-1,1] +
                            truebcr[4]      * dld[i9-2,1]  + e[i9,2];

                ldp[i9] = ldp[i9-1] + dld[i9] - dlp[i9];
                i9 = i9+1;
            endo;
            ld = cumsumc(dld);
            lp = cumsumc(dlp);
        endif;

        if graphs;
            xy(seqa(1,1,rows(lp)),lp~ld);
            xy(seqa(1,1,rows(ldp)),ldp);
        endif;

        y = dlp[(nst+1):rows(dlp),1];
        x = ones(n,1)~ldp[nst:(rows(ldp)-1),1];
        k = 1;
        do until k > klags;
            x = x~dlp[(nst+1-k):(rows(dlp)-k),1]~dld[(nst+1-k):(rows(dld)-k),1];
            k = k+1;
        endo;
        xxi = inv(x'x);
        beta = xxi*(x'y);
        ee = y-(x*beta);
        seesq = (ee'ee)/(rows(ee)-cols(x));
        rdfi = beta[2,1];
        tdfi = beta[2,1]/sqrt(seesq*xxi[2,2]);

        if imc == 1;
            tdf = tdfi;
            rdf = rdfi;
        else;
            tdf = tdf|tdfi;
            rdf = rdf|rdfi;
        endif;

        imc = imc+1;
    endo;

/* --------------------------------------------------------- */
    format /RDN 8,3;
    "distribution of r and t ";
    rdfs = sortc(rdf,1);
    tdfs = sortc(tdf,1);
    table = seqa(1,1,rows(rdfs))~rdfs~tdfs;
    if nmc < 100;
        table;
    endif;
    if nmc >= 100;
        table[rows(table)-30:rows(table),.];
    endif;
    "percentage greater than .225";
    sampby = .225;
    /* closest:
    pval = indexcat(rdfs,
           selif(rdfs,abs(sampby-rdfs).== minc(abs(sampby-rdfs))))
           /nmc*100;
     */
     /* # greater */
    pval = (indexcat(rdfs,minc(selif(rdfs,rdfs .>= sampby)))-1)/nmc*100;
    100-pval;

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

output off;