@ program to simulate ww70 with distinct individual means @
@ and autoregressive shocks that many individuals simultaneously @
library pgraph;
graphset;

let stacky[1,3]= . . . ;

@ number of time points @
t=200;          
@ initialize counters   @                       
sumy=zeros(t,1);
k=1;
@ number of individuals @
n=1000;
@ create autoregressive disturbances to the individual time series @
rho=.9;
e=recserar(rndn(t,1),0,rho);

@ create individual time series @
do while k < n;
    @ b=autoregressive parameter at the individual level @
	b=0;
	a=ones(t,1)*6*rndu(1,1)*(1-b);
    @ each individual's disturbance has a common (e) and unique component
@
	u=(rndn(t,1)+e)*.2;
	y = recserar(a+u,0,b);
    @ create panel data, using the 80th, 90th, and 100th time points @
	stacky=stacky|(y[80]~y[90]~y[100]);
	sumy=sumy+y;

	k=k+1;
endo;

dat=trimr(stacky,1,0);

print "Correlation matrix: WW70 Input ";
corrx(dat);

vc=vcx(dat);
bhat=vc[3,1]/vc[2,1];
print;
print "WW70 Estimates of time2's influence on time3:---> " bhat;

@ create aggregate time series by summing individual series @
sumy=sumy/(k-1);
shorty =sumy[101:t];
shorty1=sumy[100:t-1];
print "
       lag-1 autocorrelation in aggregate time series";
corr1=corrx(shorty~shorty1);
corr1[2,1];
print;
print " Simulation assumes that rho= " rho;
print "                     and   b= " b;
print;
months=seqa(1,1,t);
title("Aggregate Time Series");
ylabel("Mean Partisanship");
xy(months[100:t],sumy[100:t]);

/* OUTPUT

Correlation matrix: WW70 Input 

       1.0000000       0.98789963       0.98756389 
      0.98789963        1.0000000       0.98759295 
      0.98756389       0.98759295        1.0000000 

WW70 Estimates of time2's influence on time3:--->       0.99992860 

       lag-1 autocorrelation in aggregate time series
      0.83126605 

 Simulation assumes that rho=       0.90000000 
                     and   b=       0.00000000 
*/