options crt;  ? 6e6
in kopcke;
freq q;
set m = 7;
dot e s;
set m = m+1;  ? 8 for equip, 9 for structures
smpl 52:1,86:4;
x. = (1-q)*k.lag;  ? (a)
xp. = q*k.lag;  ? for (d)
smpl 56:1,79:4;
ar1 i. c x.(4,m,none) k.lag;  ? (b)
rename @logl loglu;
?
ar1 i. c x.(4,m,far) k.lag;
set lrtest = 2*(loglu-@logl);  ? (c)
cdf(chisq,df=1) lrtest;
?
ar1 i. c xp.(4,m,none) k.lag;  ? (d)
rename @logl loglxp;
print loglu loglxp;
?
?  these do loops estimate all PDLs for up to 9 lags
?
smpl 56:1,86:4;
const minsbic 999999 nl 0 mincase 0;
supres @coef;  ? turn off all output
do nlags=4,9;
   ar1(silent) i. c x.(4,nlags,none) k.lag;
   sbmin 1 1;
   ar1(silent) i. c x.(4,nlags,far) k.lag;
   sbmin 2 0;
   ar1(silent) i. c xp.(4,nlags,none) k.lag;
   sbmin 3 1;
   ar1(silent) i. c xp.(4,nlags,far) k.lag;
   sbmin 4 0;
enddo;
supres;
print minsbic nl mincase;  ? (d) -- min of @SBIC
?  print selected model
goto mincase;   ? this is like a case statement in other languages
1 ar1 i. c x.(4,nl,none) k.lag;
  goto 10;
2 ar1 i. c x.(4,nl,far) k.lag;
  goto 10;
3 ar1 i. c xp.(4,nl,none) k.lag;
  goto 10;
4 ar1 i. c xp.(4,nl,far) k.lag;
10 enddot;
?
proc sbmin case nofar;
  ? SBIC includes C, RHO, K.LAG, free params of PDL
  set @sbic = -@logl + (3+4-nofar)*log(@nob)/2;
  if (@sbic < minsbic); then; do;
    set minsbic = @sbic;
    set nl = nlags;
    set mincase = case;
  enddo;
endproc;