@ program to simulate AR(1) bias with dummy varibles @

/*
IKE  JFK  LBJ  GRF  JEC  RWR
IKET  JFKT  LBJT  RMNT  GRFT  JECT  RWRT
JFK1  LBJ1  RMN1  GRF1  JEC1  RWR1
EVENTS  WATERGATE  IRAN  ASSASSIN  TROOPS
*/

load pres[141,24]=c:\eviews\presdum2.dat;

n=rows(pres)+1;                 @ number of observations @
t=seqa(1,1,n);                  @ counter for time @
y0=60;                          @ mean for y_0 @
rho=.95;                        @ AR(1) parameter @
b=3|0;                          @ a constant (and slope)       @
repls=1;
let stackr[1,2]= 0 0;
let stacks[1,2]= 0 0;
let stackz[1,1]= 0;
z = recserar(rndn(n-1,1),0,.9);   @ pseudo independent variable @



do while repls <= 500;
    x=ones(n,1)~rndn(n,1);          @ an x variable if you want it @
    e=rndn(n,1);                    @ a N(0,1) disturbance @

    truey = recserar(x*b+e,y0,rho);

@-------------------------------------------------------@
@     add some measurement error if desired             @
@-------------------------------------------------------@

    y = truey + 0*rndn(n,1);

@-------------------------------------------------------@
@     graph data if desired                             @
@-------------------------------------------------------@

    /*
    library pgraph;
    graphset;
    xy(t,y);
    */

    lagy=trimr(y,0,1);
    y=trimr(y,1,0);


@   z=trimr(z,1,0);@
    tdum=pres;

screen off;
    { vnam,m,b1,stdb,vc,stderr,sig1,cx,rsq1,resid,dwstat }=ols(0,y,lagy);
    { vnam,m,b2,stdb,vc,stderr,sig2,cx,rsq2,resid,dwstat }=ols(0,y,lagy~tdum~z);
screen on;
/*
    print " compare AR(1) model with AR(1)+dummies model";
    print " AR(1) estimate, R2, and sigma ";
    print (b1[2,1]|rsq1|sig1)~(b2[2,1]|rsq2|sig2);
*/
    stackr=stackr|(b1[2,1]~b2[2,1]);
    stacks=stacks|(sig1~sig2);
    stackz=stackz|b2[27,1];
    print repls;
    repls=repls+1;
endo;

stackr=trimr(stackr,1,0);
stacks=trimr(stacks,1,0);
stackz=trimr(stackz,1,0);
print "Lag-y coefficients: univariate vs. dummy variable specifications";
meanc(stackr)';
print "Empirical standard errors: univariate vs. dummy variable";
stdc(stackr)';
print "Mean estimates of sigma-sq";
meanc(stacks)';
print "Empirical standard errors of sigma-sq";
stdc(stacks)';

@-------------------------------------------------------@
@     graph simulation results                          @
@-------------------------------------------------------@

library pgraph;
graphset;
begwind;
window(2,1,0);
setwind(1);
title("Empirical Distribution of Lag-Y estimates w/o Controls");
xtics(0,1,.1,1);
{ b,m,freq } = histp(stackr[.,1],20);
nextwind;
title("Empirical Distribution of Lag-Y estimates with 24 Controls");
{ b,m,freq } = histp(stackr[.,2],20);
endwind;