/* ------------------------------------------------------------------------- */
/*                          NEWMONTE.GPG                                     */
/* Dickey-fuller style bootstrap.
  This one uses ONLY PAST C in y equation
  Generates data by
    Dc(t) = 0 (y(t-1)-c(t-1)) +            a2 Dc(t-1) +..+ b2 Dy(t-1) +..+ e(t)
    Dy(t) = 0 (y(t-1)-c(t-1)) + d1 Dc(t) + d2 Dc(t-1) +..+ f2 Dy(t-1) +..+ v(t)
  generates sampling distribution of coefficient on lagged y-c.
  (regressions include a constant)
  Under this null, y-c is a unit root, so mean should be zero, but this
  is much like an augmented Dickey-fuller test so the distribution will
  be nonstandard.

  Usage:
    {bc,by}  =  montec(a,b,d,f,se,sample,draws)

  Arguments:
    a b d f   coefficients on lagged differences. d starts with c[t], all
              others start with [t-1].

    sigma   = covariance matrix of shocks
    sample  = sample size
    draws   = nuber of replications

  Returns:
    bc = ( dc on  const ~ lagged ci ~    dc-1 dy-1 dc-2 dy-2 ..)
    by = ( dy on  const ~ lagged ci ~ dc dc-1 dy-1 dc-2 dy-2 .. );
        all draws of the coefficient
*/
/* ------------------------------------------------------------------------- */
library pgraph;
graphset;

clearg test; /* = 1? only go thru once to get sizes right */
clearg bc,by; /* returned as globals to save space */
clearg bccoig, bycoig; /* defined in print subroutine */
proc(3) = regress(dc,dy,clesy,lags,force);retp(0,0,0);endp;

output file=montecar.out reset;

grf     = 0;       /* make graphs of fake data? */
getdata = 1;       /* get data?                                              */
null    = 0;       /* switch in montecarlo. =0 simultes data forcing coeff
                        on lagged ci vector to zero =1 does not force        */
force   = 1;       /* switch in initial call to regress: =1 use parameters for
                        sim from regression excluding lagged ci vector
                        (reestimate under null)                              */
replict = 0;       /* = 1: replicate ols by using original errors, coeffs    */
if replict;
    force = 0;
    null  = 1;
endif;

draws   = 2000;         /* number of monte carlo draws? */
stdates = 20;         /* number of startup dates */
savem   = 0;         /* save bccoig, bycoig ?*/
seed = 2691;
rndseed(seed);
datav   = 1|2|3|4;       /* 1: GNP82 and Cn&s.
                             2: Labor income and CN&S
                             3: GNP-GGE and CN&S
                             4: P and d from nyse
                             5: postwar P and D only */

"DRAWS" draws;
#include getdat.g;

index = 1;
do while index <= rows(datav);

    dataset = datav[index];

    if dataset == 1; freq = 1; stpnt = 47; endif; /* copied from aerppr.gpg   */
    if dataset == 2; freq = 1; stpnt = 47; endif; /* if any of that changes   */
    if dataset == 3; freq = 1; stpnt = 47; endif; /* make sure to change here */
    if dataset == 4; freq = 4; stpnt = 27; endif; /* too!                     */
    if dataset == 5; freq = 4; stpnt = 47; endif;

    if dataset == 1;""; "GNP and Consumption";          "" ; endif;
    if dataset == 2;""; "Labor Income and Consumption"; "" ; endif;
    if dataset == 3;""; "GNP-GGE and Consumption"     ; "" ; endif;
    if dataset == 4;""; "Prices and diviends";          "" ; endif;
    if dataset == 5;""; "Postwar Prices and dividends" ;"" ; endif;

    if getdata;
        call getdat(stpnt,dataset);
    endif;
    #include datsetup.gpg

    /* ----------------------- */
    /* redo var for comparison */
    /* ----------------------- */

    {bmat,sigma,r,dxir,xir,statcom,e1} = var(x,a,lags,horiz,names,0,prnt);
    sampbc = bmat[2,1];
    sampby = bmat[2,2];


    /* ------------------------------ */
    /* get, print starting parameters */
    /* ------------------------------ */

    dx = x[2:rows(x),.]-x[1:rows(x)-1,.];
    clesy = x[2:rows(x),1]-x[2:rows(x),2];
    dc = dx[.,1];
    dy = dx[.,2];
    {truebc,trueby,e} = regress(dc,dy,clesy,lags,force);
    "";
    save cerr_1 = e;  @ for stockmc.gpg @
    @ save truebc;    for stockmc.gpg    @
    @ save trueby;                       @

    format /RD 12,6;
    "parameters for simulations";
    label = "const"~"c-y(-1)"~"dc(-1)"~"dc(-2)"~"dy(-1)"~"dy(-2)";
    $label;
    " bc" truebc';
    " by" trueby';
    "var/cov matrix (1/T)" ;
    e'e/rows(e);

    call montec(truebc~trueby,draws,x,e,null,stdates,replict,grf);
    gosub printit;
    if savem;
        if dataset == 1; save bccoig1 = bccoig; save bycoig1 = bycoig; endif;
        if dataset == 2; save bccoig2 = bccoig; save bycoig2 = bycoig; endif;
        if dataset == 3; save bccoig3 = bccoig; save bycoig3 = bycoig; endif;
        if dataset == 4; save bccoig4 = bccoig; save bycoig4 = bycoig; endif;
    endif;
    index = index + 1;
endo;
output off;
end;

/* ------------- */
/* print results */
/* ------------- */

printit:

format /RZ 12,3;
label = "const"~"c-y(-1)"~"dc(-1)"~"dc(-2)"~"dy(-1)"~"dy(-2)";

"mean OLS coefficients" ;
$label;

meanbc = meanc(bc);
meanbc';
meanby = meanc(by);
meanby';
"";

if draws > 1;

    "standard deviation of OLS coefficients";
    $label;
    sdbc = stdc(bc);
    sdbc';
    sdby = stdc(by);
    sdby';
     "";
    "c and y cointegration t ratios: (sample b - mean bx)/stddev(bx)";
     (sampbc-meanbc[2])/sdbc[2];
     (sampby-meanby[2])/sdby[2];
endif;

if draws >= 100 ;
    bccoig = sortc(bc[.,2],1);
    bycoig = sortc(by[.,2],1);
    label = "5%"~"10%"~"50%"~"10%"~"5%";
    points = draws.*(1/100~1/20~1/10~1/2~9/10~19/20~99/100);
    "distribution of coefficient on lagged ci vector under 0 null";
    $label;
    "bc"
    bccoig[points]';
    "by"
    bycoig[points]';
    "";
    if draws > 1000;
        clear bc,by;
    endif;
    "p-values of sample bc, by (order of closest monte carlo bx)" ;
    pval = indexcat(bccoig,
           selif(bccoig,abs(sampbc-bccoig).== minc(abs(sampbc-bccoig))))
           /draws*100;
    pval;
    pval = indexcat(bycoig,
           selif(bycoig,abs(sampby-bycoig).== minc(abs(sampby-bycoig))))
           /draws*100;
    pval;
    "top and bottom 10 coefficients on cointegrating vector";
    "  bc      by ";
    ((bccoig[1:10]~bycoig[1:10])|
   (bccoig[rows(bccoig)-10:rows(bccoig)]~bycoig[rows(bccoig)-10:rows(bccoig)]));
endif;

return;


/* -------------------------------------------------------------------------- */
/*                               PROCS                                        */
/* -------------------------------------------------------------------------- */


proc(3) = regress(dc,dy,clesy,lags,force);
    local rhv,j,sample,bc,by,e;
    sample = rows(dy)-lags;

    if force == 0;
       Rhv = ones(sample,1)~(clesy[2:sample+lags-1]);
    endif;
    if force == 1;
        Rhv = ones(sample,1);
    endif;

   if lags > 0;
      j = 1;
      do while j <= lags;
         Rhv = Rhv~dc[1+lags-j:sample+lags-j];
         j = j+1;
      endo;
      j = 1;
      do while j <= lags;
         Rhv = Rhv~dy[1+lags-j:sample+lags-j];
         j = j+1;
      endo;
    endif;

    bc = (dc[1+lags:sample+lags]/rhv);
    by = (dy[1+lags:sample+lags]/rhv);
    e  = (dc[1+lags:sample+lags]-rhv*bc) ~ (dy[1+lags:sample+lags]-rhv*by);

    if force == 1;
        bc = bc[1]|0|bc[2:rows(bc)];
        by = by[1]|0|by[2:rows(by)];
    endif;

retp(bc,by,e);
endp;


proc(0) = montec(bmat,draws,x,e,null,stdates,replict,grf);
clearg const,a,b,d,f,g,lags,sig2,trial,t,rhv,j,
      betac,betay,sample,trash1,trash2,trash3;
clearg dcl,dyl,clesyl,eandv,dxl;

      lags = (rows(bmat)-2)/2; "lags" lags;
      sample = rows(e);

      const = bmat[1,.];
      g = bmat[2,.];
      a = bmat[2+1:2+lags,1];
      b = bmat[2+lags+1:2+2*lags,1];
      d = bmat[2+1:2+lags,2];
      f = bmat[2+lags+1:2+2*lags,2];

      dxl = x[2:rows(x),.]-x[1:rows(x)-1,.];
      clesyl = x[2:rows(x),1]-x[2:rows(x),2];
      dcl = dxl[.,1];
      dyl = dxl[.,2];
      if stdates >0;
        dcl = dcl|zeros(stdates,1);
        dyl = dyl|zeros(stdates,1);
        clesyl = clesyl|zeros(stdates,1);
      endif;
      bc = zeros(draws,2+2*lags);
      by = zeros(draws,2+2*lags);

      format /RD 10,0;

      trial = 1;

      do while trial <= draws;

        if (trial == 1) or ((trunc(trial/100) == (trial/100))) or 0;
            output off;
            "draw #" trial "of" draws;
            output on;
        endif;

       /* --------------- */
       /* generate errors */
       /* --------------- */

       /* bootstrap */
        if not replict;
            eandv = submat(e,
                trunc( rows(e)*rndu(stdates+rows(e),1)+1),1|2);
        endif;

     /* test to see if it replicates ols */
        if replict;
           eandv = e;
        endif;

     /* ------------- */
     /* generate data */
     /* ------------- */

        t = lags+1;
        do while t <= sample+lags+stdates;
            dcl[t,.] = const[1,1] + null*g[1]*clesyl[t-1]
                      + a'rev(dcl[t-lags:t-1])
                      + b'rev(dyl[t-lags:t-1]) + eandv[t-lags,1];
            dyl[t,.] = const[1,2] + null*g[2]*clesyl[t-1]
                      + d'rev(dcl[t-lags:t-1])
                      + f'rev(dyl[t-lags:t-1]) + eandv[t-lags,2];
            clesyl[t] = clesyl[t-1] + dcl[t]-dyl[t];
            t = t+1;
        endo;

        if grf;
            title("log ratio");
            xy(seqa(1,1,rows(clesyl)),clesyl);
        endif;

        {trash1, trash2, trash3} = regress(
            dcl[1+stdates:rows(dcl)],
            dyl[1+stdates:rows(dyl)],
            clesyl[1+stdates:rows(clesyl)],lags,0);

        bc[trial,.] = trash1';
        by[trial,.] = trash2';

        trial = trial+1;
    endo;
endp;

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