calendar 1947 1 4
allocate 1991:3
*allocate 1988:3

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

set y = gnp82

statistics y

* Plot of the data
graph(key=upleft)
# y

print / y

* Compute and graph autocorrelations:
correlate(partial=PACF, QSTATS) y / ACF
graph(key=below,style=bar,nodates,min=-1.0,max=1.0,number=1) 2
# ACF
# PACF

*** Take log if needed (ALL values are positive and big.)
*log  y / log_y

*** Take the first difference if non-stationary
diff y / gw
*diff log_y / gw

*** Set: "mydata" (no change), "log_data" (logged), or "gw" (FD)
set y = gw

* Searching for a minimum of SBC (also called, SIC)

com sbcmin = 100000000., sbcp = 0, sbcq = 0
com aicmin = 100000000., aicp = 0, aicq = 0

do i=0,3
do j=0,3

boxjenk(ar=i,ma=j,constant,noprint)  y  / resids
compute sbc = log(%rss/%nobs) + %nreg*log(%nobs)/(%nobs)
compute aic = log(%rss/%nobs) + 2*(%nreg)/%nobs
*compute aic = %nobs*log(%rss) + 2*%nreg
*compute sbc = %nobs*log(%rss) + %nreg*log(%nobs)

 display 'p ='  i  ' q =' j  '  AIC = ' aic     '  SBC = ' sbc
 if sbc < sbcmin
 {
  compute sbcp = i
  compute sbcq = j
  compute sbcmin = sbc
  }
 if aic < aicmin
 {
  compute aicp = i
  compute aicq = j
  compute aicmin = aic
  }
endo j
endo i
disp
disp 'Selected values of p and q using SBC = ' sbcp  sbcq
disp 'Selected values of p and q using AIC = ' aicp  aicq

disp

* Now, Estimating the ARMA(p,q) model
***  using the selected values of p (AR=p) and q (MA=q)
disp '***** Using the selected orders using SBC'

boxjenk(define=eq1, AR=aicp, MA=aicq, constant) y / resids
compute sbc = log(%rss/%nobs) + %nreg*log(%nobs)/(%nobs)
compute aic = log(%rss/%nobs) + 2*(%nreg)/%nobs
*compute aic = %nobs*log(%rss) + 2*%nreg
*compute sbc = %nobs*log(%rss) + %nreg*log(%nobs)

display  'aic = ' aic
* Checking the residuals
cor(partial=PACF, qstats, number=36, span=12, dfc=%nreg) resids /  ACF

graph(key=below,style=bar,nodates,min=-1.0,max=1.0,number=1, $
 header = 'ACF and PACF of Residuals') 2
# ACF
# PACF

* Forecasting 1-12 values starting from the 1001-th observation and save it
forecast(print) 1 12 1988:4
# eq1 y_for

smpl 1947:1 1991:3
graph(style=line) 2
# y
# y_for