options double crt nodate;
name bjfpac 'Box-Jenkins series F -- Partial AutoCorrelations';
?
? See McCullough, B.D., "Algorithm choice for (partial) autocorrelation
? functions," Journal of Economic and Social Measurement 24, 1998,
? pp.265-278.
?  by Clint Cummins  (updated 10/1999)

? The BJIDENT command in TSP was revised on 3/21/1998 to make the
? BURG method the default (instead of Yule-Walker), and also to
? offer the OLS method.  The code below can be run on either
? newer versions of TSP, or pre-3/98 versions of TSP, by setting
? the value of the flag  pre398  .
const pre398,0;  ? set to 1 if you have TSP dated before 3/21/98
const maxlag,20;

load;

title 'PAC - various methods, in rough order of accuracy';

title 'PAC - Yule-Walker method';
? Yule-Walker method (least accurate; subject to numerical problems
? when the condition number of the autocorrelation matrix is high)
? (was default method in TSP until 3/21/98; replaced by Burg method)
? Lag  McCullough     TSP
?  1     -.39      -.38988
?  2      .18       .17971
?  3      .00       .0022645
?  4     -.04      -.044277
?  5     -.07      -.069406
?  6     -.12      -.12062
?  7      .02       .019680
?  8      .01       .0048884
?  9     -.06      -.056497
? 10      .00       .0037055
? 11      .14       .14280
? 12     -.01      -.0094066
? 13      .09       .091964
? 14      .17       .16693
? 15     -.00      -.0012949     
? 16      .22       .22069     
? 17      .05       .052808     
? 18     -.10      -.10519     
? 19      .04       .042236     
? 20      .05       .049782     
if pre398; then;
   bjident(nlag=maxlag,nlagp=maxlag,nocumpl,noplot) f;
else;
   bjident(pacmeth=yulewalk,nlag=maxlag,nlagp=maxlag,nocumpl,noplot) f;
copy @pac pacyw;

title 'PAC - OLS method';
? OLS - more accurate method (pick off highest AR coefficient,
? same as conditional ML estimation -- does not impose stationarity)
?    Note:  it is also possible to run a "reverse" OLS estimation
? on leads instead of lags; the estimates will be different.
? Both ways are consistent.
?      McCullough     TSP
?  1     -.42      -.42490
?  2      .19       .19040
?  3      .01       .0083800
?  4     -.07      -.065703
?  5     -.07      -.069805
?  6     -.15      -.14537
?  7      .05       .045633
?  8     -.00      -.0020036
?  9     -.10      -.098443
? 10      .05       .049197
? 11      .18       .18271
? 12     -.06      -.054537
? 13      .09       .093309
? 14      .18       .18338
? 15      .01       .014066
? 16      .43       .43275
? 17      .01       .010086
? 18     -.14      -.14126
? 19      .11       .11339
? 20      .18       .17729
if pre398; then; do;
  do p=1,maxlag;
    set mp = -p;
    olsq(silent) f f(-1)-f(mp) c;
    set @pac(p) = @coef(p);
  enddo;
enddo;
else;
  bjident(pacmeth=ols,nlag=maxlag,nlagp=maxlag,nocumpl,noplot) f;
copy @pac pacols;

title 'PAC - Burg method';
?      McCullough     TSP
?  1     -.41      -.40667
?  2      .18       .18301
?  3      .01       .010302
?  4     -.07      -.066470
?  5     -.07      -.065889
?  6     -.14      -.14488
?  7      .05       .045059
?  8      .01       .0054400
?  9     -.10      -.10481
? 10      .05       .053789
? 11      .16       .16216
? 12     -.03      -.030592
? 13      .07       .072159
? 14      .16       .16477
? 15      .02       .018646
? 16      .32       .32302
? 17      .00      -.0036412
? 18     -.18      -.17796
? 19      .11       .10838
? 20      .08       .075237
if pre398; then; do;
   title 'Burg method not implemented before 3/98';
   mform(nrow=maxlag,ncol=1) @pac = @na;
enddo;
else;
   bjident(pacmeth=burg,nlag=maxlag,nlagp=maxlag,nocumpl,noplot) f;
copy @pac pacburg;

title 'PAC - Exact ML method';
?  The exact ML method requires a lot of computation, so it
?  is not very practical.  Discussed only briefly by McCullough.
?     Box-Jenkins    TSP/BJEST
?     (uncond SSR)   (exact ML)  ifconv
?  1     -.40        -.41908       1
?  2      .19         .18734       1
?  3      .01         .010527      0
?  4     -.07        -.072289      1
?  5     -.07        -.069309      1
?  6     -.15        -.15332       0
?  7      .05         .046334      0
?  8      .00         .0042786     1
?  9     -.10        -.11898       1
? 10      .05         .065427      0
? 11      .18         .16419       0
? 12     -.05        -.034254      1
? 13      .09         .077616      0
? 14      .18         .17602       0
? 15      .01         .016050      0
? 16       .          .35492       1
? 17       .         -.012877      0
? 18       .         -.19804       0
? 19       .          .12764       0
? 20       .          .084537      0
? (p.66 in Box-Jenkins(1976)).
mmake pacussr -.40 .19 .01 -.07 -.07 -.15 .05 .00 -.10 .05
               .18 -.05 .09 .18 .01 @na @na @na @na @na;
? most accurate method (pick off highest AR coefficient,
? using exact ML estimation)
copy @pac paceml;
do p=1,maxlag;
  bjest(nar=p,const,exactml,silent,
        nocumpl,noplot,tol=1e-7,grad=c4,maxit=50) f;
  print p,@ifconv;
  set paceml(p) = @coef(p);
enddo;
mmake pacs pacyw pacols pacburg paceml;
smpl 1,maxlag;
unmake pacs yulewalk ols burg exactml;
print       yulewalk ols burg exactml;
end;

noprint;
smpl 1,70;
read f;  ? p.530 - yields from batch chemical process
47 64 23 71 38 64 55 41 59 48 71 35 57 40 58
44 80 55 37 74 51 57 50 60 45 57 50 45 25 59
50 71 56 74 50 58 45 54 36 54 48 55 45 57 50
62 44 64 43 52 38 59 55 41 53 49 34 35 54 45
68 38 50 60 39 59 40 57 54 23;