* ARMA_simulation_hw.PRG
* Box Jenkins Procedures

* Reading the data and plot
all 1000
open data arma_simulation.dat
data(format=prn,org=obs) / Y1 y2 y3 y4

set y = y1

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

* Estimating an AR(2) model
boxjenk(constant, ar=2) y / resids

* Checking residuals
cor(partial=pacf,qstats,number=24,span=8,dfc=%nreg) resids / ACF

* AIC and BIC
compute aic = log(%rss/%nobs) + 2*(%nreg)/(%nobs)
compute BIC = log(%rss/%nobs) + %nreg*log(%nobs)/(%nobs)
display 'AIC = ' aic     '  BIC = ' BIC

* Estimating ARMA(1,1) model
boxjenk(noprint,ar=1,ma=1, constant) y
compute aic = log(%rss/%nobs) + 2*(%nreg)/(%nobs)
compute BIC = log(%rss/%nobs) + %nreg*log(%nobs)/(%nobs)
display 'AIC = ' aic     '  BIC = ' BIC

* Re-estimating ARMA(1,1) model, while definfing it as eq1 for forecasting
boxjenk(define=eq1,constant,ar=1,ma=1) y / resids
compute aic = log(%rss/%nobs) + 2*(%nreg)/(%nobs)
compute BIC = log(%rss/%nobs) + %nreg*log(%nobs)/(%nobs)
display 'AIC = ' aic     '  BIC = ' BIC
  * Chekcing residuals
cor(partial=pacf,qstats,number=24,span=8,dfc=%nreg) resids / ACF

* The above forecasting can be done by using BJFORE.SRC
source(noecho) bjfore.src
@bjfore(ars=1,mas=1,constant) y 101 108 fores
print / fores

** Do-loop procedure
do i=0,5
do j=0,5
 boxjenk(ar=i,ma=j,constant,noprint) y
 compute aic = log(%rss/%nobs) + 2*(%nreg)/(%nobs)
 compute BIC = log(%rss/%nobs) + %nreg*log(%nobs)/(%nobs)
 display 'p ='  i  ' q =' j  '  AIC = ' aic     '  BIC = ' BIC
endo j
endo i

** Alternatively,
com BICmin = 100000000., BICp = 0, BICq = 0
do i=0,5
do j=0,5
 boxjenk(ar=i,ma=j,constant,noprint) y
 compute aic = log(%rss/%nobs) + 2*(%nreg)/(%nobs)
 compute BIC = log(%rss/%nobs) + %nreg*log(%nobs)/(%nobs)
 display 'p ='  i  ' q =' j  '  AIC = ' aic     '  BIC = ' BIC
 if BIC < BICmin
 {
  compute BICp = i
  compute BICq = j
  compute BICmin = BIC
  }
endo j
endo i
disp
disp 'Selected values of p and q using BIC = ' BICp  BICq
disp

* Now, Re-Estimating the ARMA(1,0) model
boxjenk(define=eq2,ar=1,ma=0, constant) y / resids
compute BIC = log(%rss/%nobs) + %nreg*log(%nobs)/(%nobs)
display  '  BIC = ' BIC

* Chekcing residuals
cor(partial=pacf,qstats,number=24,span=8,dfc=%nreg) resids / ACF

* Forecasting 1-15 values starting from the 101th observation and save it
forecast(print) 1 15 1001 1015
# eq2 y_for

smpl 1 1015
graph(style=line) 2
# y
# y_for