* ARMA3.PRG
* Box Jenkins Procedures (Example)

calendar 1947 1 4
allocate 1988:3

open data gnp82.dat
data(format=prn,org=cols) / gnp82

set y = gnp82

statistics y

graph 1
# y

* Checking ACFs and plot
cor(partial=pacf) y / ACF
   * Modify the next line  ** No differencing will be made.

source(noecho) bjident.src
@bjident y

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

cor(partial=pacf) gw / ACF
   * Modify the next line  ** No differencing will be made.

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,5) model
boxjenk(define=eq1,ar=3,ma=5, constant) gw / resids
compute sbc = log(%rss) + %nreg*log(%nobs)/(%nobs)
display  '  SBC = ' sbc

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

* Forecasting 1-16 values starting from 1988:4th observation and save it
forecast(print) 1 16 1988:4
# eq1 gw_for

smpl 1947:1 1992:3
graph(style=line)   2
# gw
# gw_for