options crt double;
name kald 'Kalman with initially singular dummy variable and 2 breaks';
const t0,100 t1,200 t,300;
const g0,0 g1,1 g2,0;
smpl 1,t;
random(seedin=98421) e;
trend obs;
d = obs > t1;
y = g0*(obs <= t0) +
    g1*((t0 < obs) & (obs <= t1)) +
    g2*d + e;

? Default without prior  (error -- singular prior variance)
kalman(noetrans) y c d;

? with diffuse prior  (works, but not exactly the same as OLS)
mmake bpr 0 0;
mat vbpr = 10000*ident(2);
kalman(noetrans,bprior=bpr,vbprior=vbpr) y c d;
unmake @state alpha beta;
print alpha beta;
?plot(preview) alpha;
?plot(preview) beta;
cusum = @recres;
smpl 2,t;
cusum = cusum(-1) + @recres;
smpl 1,t;

plots cusum;  ? cusum fails in OLS
olsq y c d;
set ssrols = @ssr;

random f;
olsq(silent) y c f;  ? just to make OLSQ store bounds for cusum plot
?plot(preview,origin) cusum @csub5% @cslb5%;

? Use Proc Recres below to compute recursive residual by brute force
? (inverting X'X at each observation).
list xs c d;
recres er;
?print er;
set er(1) = 0;  ? set first g recursive residuals to zero; g = rank(X(1:k))
mat ssrr = er'er;
print ssrr ssrols;  ? SSR matches OLS to 5 digits

Proc Recres er;
?
?  Computes recursive residuals even if the first K observations
?  are singular.
er = 0;
copy y yt;
mmake x xs;
length xs nx;
mat xpx = 0*ident(nx);
copy xpx v;
mform(nrow=nx,ncol=1) b=0;
supres smpl;
set nob = @nob;
do i=1,nob;
  smpl i,i;
  mmake x xs;
? er(t) = [y(t) - x(t)*b(t-1)]/sqrt[x(t)*v(t-1)*x(t)'+1]
  mat e1 = yt - x*b;
  ?mat den = sqrt( x*yinv(xpx)*x' + 1 );
  mat den = sqrt( x*v*x' + 1 );
  set er(i) = e1/den;
  ?print i e1 den er;
  mat xpx = xpx + x'x;
  ?mat b = b + yinv(xpx)*x'e1;
  mat v = yinv(xpx);
  mat b = b + v*x'e1;
  ?print b xpx;
enddo;
smpl 1,nob;
endproc;