new;
library pgraph;

rhoin=0.94;
n=100;
t=100;
rho={};
rhof={};
rhoinc={};

i=1;
do until i>=1001;

rhoinb=rhoin+uni(0.059,-0.059,n);             @sample from a uniform distribution@                                                            @vector of heterogeneous coeffients drawn from a uniform distribution@
eta=rndn(n,1);                                 @fixed effects@                                                                         @ N(0,1) fixed effects@
dat=average(rhoinb,t,1,eta);        @this produces an N,T panel of Within transformed data@
dat1=average(rhoinb,t,1,zeros(n,1));    @dat1  no fixed effects@                      
dat=packr(dat);
y=dat[.,2];
x=dat[.,3];
b=inv(x'x)*(x'y);                      @fe estimates@                                                                                                           @Fixed  estimates@
rho=rho|b[1];								@save persistence estimates@
dat1=packr(dat1);
y1=dat1[.,2];
x1=dat1[.,3];
bf=inv(x1'x1)*(x1'y1);
rhof=rhof|bf[1];
rhoinc=rhoinc|meanc(rhoinb);			 @save true value@


"iter";;i;
i=i+1;
endo;
 
"Estimated persistence fixed effects (FE in DGP)";;meanc(trimr(rho,100,0));
 "Estimated persistence fixed effects (no FE in DGP)";;meanc(trimr(rhof,100,0));
"Mean Persistence";;meanc(trimr(rhoinc,100,0));


 

proc(1)=average(rho,t,sigma,eta);		 @rho=N by 1 vector of ar parameters, t=T,sigma=error variance,eta=fixed effects@
local y,i,yy,index,xx,x;
y={};
x={};
index={};
i=1;
do until i>=rows(rho)+1;
yy=ar1(rho[i],t+50,sigma,eta[i]);				 @call procedure to generate AR series@
y=y|(trimr(yy,50,0)-meanc(trimr(yy,50,0)));		 @save de-meaned dependent variable after deleting 1st 50@
xx=lag1(trimr(yy,50,0));						 @create lagged series@
x=x|(xx-meanc(packr(xx)));						 @save its within transform@
index=index|ones(rows(trimr(yy,50,0)),1).*i;	 @create an index (not used here but useful for RCM etc)@
i=i+1;
endo;
retp(index~y~x);
endp;

proc ar1(b,t,sigma,cont);							   @ note b=AR parameter,t=length,sigma=error variance,cont=constant term@
local y,i,v;
y=zeros(t,1);
v=norm(t,0,sigma);										@normal(0,sigma) error@

i=2;
do until i>=t+1;
y[i]=cont+(y[i-1]*b)+v[i];
i=i+1;
endo;
retp((miss(y,0)));
endp;

proc(1)=norm(nn,nor_mean,nor_var);
 local st_norm,normal;
  st_norm = rndn(nn,1) ;   
  normal = (rndn(nn,1) * sqrt(nor_var)) + nor_mean;
  retp(normal);
endp;

proc(1)=uni(max,min,n);
local st_uni,uniform;
  st_uni = rndu(n,1) ;   @U(0,1)@
  uniform = (rndu(n,1) * (max - min)) + min;   @U(min,max)@
retp(uniform);
endp;