repls=1000;
j=0;
let stackp[1,4]= . . . .;
do while j<repls;


a0= 0;       @ voter's best guess of mean: prior                      @
p0=25;       @ uncertainty of prior RELATIVE concerning differential  @
q=  4;       @ disturbance variance: in underlying party differential @
h= 12.5;     @ sampling variance in party differential at each point  @
gam=1;       @ autoregressive coefficient on underlying differential  @

cases=10;

@---------------------------------------------------@
@ alpha represents the true underlying differential @
@ y is what is observed at each time point          @
@---------------------------------------------------@
alpha=recserar(sqrt(q)*rndn(cases,1),0,gam);
y=alpha+sqrt(h)*rndn(cases,1);

@--------------------------------------------------------@
@ set up the Kalman Filter weighting scheme:             @
@   Y_t will be weighted by k and alpha_hat[t-1] by (1-k)@
@   to produce the guess about the true alpha_t          @
@--------------------------------------------------------@
k=(gam^2*p0+q)./(gam^2*p0+q+h);
a1=gam*a0+k*(y[1,1]-gam*a0);
p1=h*k;

@ run the KF iteratively for all observations @

ahat=a1;
stackk=k;
i=0;
do while i < cases-1;

	k=(gam^2*p1+q)./(gam^2*p1+q+h);
	a1=gam*a1+k*(y[2+i,1]-gam*a1);
	p1=h*k;

	ahat=ahat|a1;
	stackk=stackk|k;

	i=i+1;
endo;

/*
library pgraph;
graphset;
@ display the progression of weights over observations @
xy(seqa(1,1,cases),stackk[1:cases,.]);
@ display observed, true, and fitted values  @
_plegctl=1;
_plegstr="observed Y\000Filtered Alpha Estimate\000True Alpha";
xy(seqa(1,1,rows(ahat)),y~ahat~alpha);

@ display table of input and output @
format /rd 10,5;
s1="y";s2="ahat";s3="alpha";s4="k";
print $s1~s2~s3~s4;
let blank[1,1]= . ;
print blank~a0~blank;
print y[1:10]~ahat[1:10]~alpha[1:10]~stackk[1:10];
print;

*/
stackp=stackp|(y[5]~ahat[5]~alpha[5]~meanc(y[1:5]));

j=j+1;
endo;
stackp=trimr(stackp,1,0);
format /rd 13,5;
print "means, standard deviations, mad, RMSE";
s1="y";s2="ahat";s3="alpha";s4="meanc(y1:5)";
print $s1~s2~s3~s4;
meanc(stackp)';
stdc(stackp)';
meanc(abs(stackp[.,3]-stackp))';
sqrt(meanc((stackp[.,3]-stackp)^2)');