@ program to to kf assuming unevenly spaced observations @ a0=rndn(1,1); @ voter's best guess of mean: prior @ p0=100; @ uncertainty of prior RELATIVE concerning differential @ q= .05; @ disturbance variance: in underlying party differential @ h= 1; @ sampling variance in party differential at each point @ gam=.9; @ autoregressive coefficient on underlying differential @ cases=100; format /rd 9,4; @---------------------------------------------------@ @ 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); days=seqa(1,1,cases); */ alpha=recserar(sqrt(q)*rndn(3*cases,1),0,gam); y=alpha+sqrt(h)*rndn(3*cases,1); days=seqa(1,1,3*cases); selectr=round(-.1+rndu(3*cases,1)); alpha=selif(alpha,selectr .eq 1); days=selif(days,selectr .eq 1); y=selif(y,selectr .eq 1); alpha=alpha[1:100]; days=days[1:100]; y=y[1:100]; @--------------------------------------------------------@ @ 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; qp=q; gamp=gam; i=0;j=0;prevday=0; do while i < cases-1; elapse = days[i+1,1]-prevday; gam = gamp^elapse; qsum=0; do while j < elapse; qsum=qsum+gamp^(2*j); j=j+1; endo; q = qsum*qp; k=(gam^2*p1+q)./(gam^2*p1+q+h); a1=gam*a1+k*(y[2+i,1]-gam*a1); print elapse~a1~k~gam~q~qsum; p1=h*k; ahat=ahat|a1; stackk=stackk|k; prevday=days[i+1,1]; i=i+1;j=0; endo; library pgraph; graphset; @ display the progression of weights over observations @ xy(days,stackk[1:cases,.]); @ display observed, true, and fitted values @ _plegctl=1; _plctrl={ 1,1,1 }; _plegstr="observed Y\000Filtered Alpha Estimate\000True Alpha"; xy(days,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];