* ARMA_Rail.PRG
* Box Jenkins Procedures (Example)

********* 1. Reading the data file

calendar 1965 1 4
allocate 1978:4

open data rail.txt
data(format=prn,org=cols) / rail

* use the name, y
set y = rail
statistics y

graph 1
# y

********* 2. Checking ACFs and plot of the data
cor(partial=pacf, qstats,number=14) y

source(noecho) bjident.src
@bjident y

* Checking ACFs of the differenced data (with option, diff=1)
source(noecho) bjident.src
@bjident(diff=1) y

********* 3. Estimating an AR(1) model
boxjenk(define=eq1,ar=1,ma=0,constant) y / res1

* Checking ACFs of the residual
cor(partial=pacf, qstats,number=14) res1

* Forecasting 4 values
forecast(print) 1 4 1979:1
# eq1 y_for

* Alternatively, forecast using the code, bjfore.src
source(noecho) bjfore.src
@BJFORE(DIFFS=0,ARS=1, MAS= 1) y 1979:1 1979:4 y_for1 res1

print 1979:1 1979:4 y_for1

********* 4. Estimating an ARMA(1,1) model
boxjenk(define=eq2,ar=1,ma=1,constant) y / res2

* Checking ACFs of the residual
cor(partial=pacf, qstats,number=14) res2

* Forecasting 4 values
forecast(print) 1 4 1979:1
# eq2 y_for

********* 5. Estimating an ARIMA(0,1,1) model
* difference the data first, and call it dy.
diff y / dy

* Checking ACFs of the differenced data
cor(partial=pacf, qstats,number=14) dy

boxjenk(define=eq3,ar=0,ma=1,constant) dy / res3

* Checking ACFs of the residual
cor(partial=pacf, qstats,number=14) res3

* Forecasting 4 values
forecast(print) 1 4 1979:1
# eq3 y_for

***** Do-loop procedure with y to find p,q values with a min SBC

com sbcmin = 100000000., sbcp = 0, sbcq = 0
do i=0,3
do j=0,3
 boxjenk(ar=i,ma=j,constant,noprint) y
 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(1,1) model
boxjenk(define=eq1,ar=1,ma=1, constant) y / resids
compute sbc = log(%rss/%nobs) + %nreg*log(%nobs)/(%nobs)
display  '  SBC = ' sbc

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

* Forecasting 4 values
forecast(print) 1 4 1979:1
# eq1 y_for