* ARMA4_ex.PRG
* Box Jenkins Procedures (Example)

calendar 1960 1 4
allocate 1991:4

open data us.prn
data(format=prn,org=cols) / M1 TBILL r3 r10 GDPDEF GDP85 GOVT

set y = gdp85

statistics y

graph 1
# y

* Checking ACFs and plot
cor(partial=pacf) y

source(noecho) bjident.src
@bjident y

** Taking the first difference of logged data
log  y / ly
diff ly / gw

cor(partial=pacf) gw

source(noecho) bjident.src
@bjident gw

** Do-loop procedure

com sbcmin = 100000000., sbcp = 0, sbcq = 0
do i=0,5
do j=0,5
 boxjenk(ar=i,ma=j,constant,noprint) gw
 compute aic = log(%rss/%nobs) + 2*(%nreg)/(%nobs)
 compute sbc = log(%rss/%nobs) + %nreg*log(%nobs)/(%nobs)
 display 'p ='  i  ' q =' j  '  AIC = ' aic     '  SBC = ' sbc
 if sbc < sbcmin
 {
  compute sbcp = i
  compute sbcq = j
  compute sbcmin = sbc
  }
endo j
endo i
disp
disp 'Selected values of p and q using SBC = ' sbcp  sbcq
disp

* Now, Re-Estimating the ARMA(3,2) model
boxjenk(define=eq1,ar=3,ma=2, constant) gw / resids
compute sbc = log(%rss/%nobs) + %nreg*log(%nobs)/(%nobs)
display  '  SBC = ' sbc

* Chekcing residuals
cor(partial=pacf,qstats,number=24,span=8,dfc=%nreg) resids
source(noecho) bjident.src
@bjident resids

* Forecasting 1-16 values starting from the 101th observation and save it
forecast(print) 1 16 1992:1
# eq1 dy_for

print 1991:1 1991:4 ly
*set 1992:1 1995:4 ly_for = ly(-1) + dy_for
*errors
set y_for = exp(ly_for)
print 1992:1 1995:4 y  ly_for  dy_for  y_for

source(noecho) bjfore.src
@bjfore(diff=1, constant, ars=3, mas=2) ly 1992:1 1995:4 ly_for
* The above regression uses the un-differenced loggged series.

*set 1992:1 1995:4 y_for = exp(ly_for)
*errors
print 1992:1 1995:4 ly_for y_for