*
*  ARMAX.PRG
*  Example 9.4 (see also ARIMA.PRG)
*
* Combined Regression/Time Series Model (ARMAX)
*
* See Section 19.6 in Pindyck and Rubinfeld (4th Edition, 1998)
* This builds on OLS and ARIMA work demonstrated in BASICS.PRG
* and ARIMA.PRG.
* (P&R use data through March, 1996, ours ends in February, 1996)
*********************************************************************

*********************************************************************
*** Set the CALENDAR and ALLOCATE range. CAL to start in January, 1960:
*** and read in the Data:
**********************************************************************

calendar 1960 1 12
allocate 1999:9
open data basics.wks
data(format=wks,organization=observation) / rate m1 m2 ip ppi

**********************************************************************
*** Generate new series
**********************************************************************
*  M2(t) - M2(t-1) divided by M2(t-1), and
*  PPI(t) - PPI(t-1) divided by PPI(t-1):
set grm2 = (m2 - m2{1})/m2{1}
set grppi = (ppi - ppi{1})/ppi{1}

display 'Simple OLS model:'
linreg(print) rate * 1995:8 olsres
# constant ip grm2 grppi{1}

correlate(partial=olspcorrs,number=35) olsres / olscorrs
graph(key=below,style=bar,nodates,min=-1.0,max=1.0,number=1) 2
# olscorrs 2 25
# olspcorrs 2 25

* After experimenting with other models, P&R settle on ARIMA(8,0,2):
boxjenk(ar=8,ma=2,constant,iters=100,define=armax) olsres / armares
correlate(partial=armapcorrs,number=24) armares / armacorrs
graph(window='ARMA(8,0,2) Correlations',key=below,style=bar,nodates,min=-1.0,max=1.0,number=1) 2
# armacorrs 2 25
# armapcorrs 2 25

* Estimate P&R's ARMA model with regressors as inputs:
display "P&R's ARMAX BOXJENK:"
nlpar(sub=100)
boxjenk(inputs=3,ar=8,ma=2,constant,iters=100,define=armax) rate / armaxres
# ip 0 0 0
# grm2 0 0 0
# grppi 0 0 1

* "Forecasts" (actually a series of 1-step forecasts) in levels:
smpl 1961:11 1995:8
steps 1
# armax arma_fore
graph(key=below,header='Interest Rate One-Step Forecasts, Actuals') 2
# rate
# arma_fore
smpl

* Now try a much simpler model:
boxjenk(ar=2,ma=1,constant,iters=100,define=simple) olsres / armares
correlate(partial=armapcorrs,number=24) armares / armacorrs
graph(window='ARMA(2,0,0) Correlations',key=below,style=bar,nodates,min=-1.0,max=1.0,number=1) 2
# armacorrs 2 25
# armapcorrs 2 25

display "Smaller ARMAX model:"
boxjenk(inputs=3,ar=2,ma=1,constant,iters=100,define=armaxsmall) rate / armaxsres
# ip 0 0 0
# grm2 0 0 0
# grppi 0 0 1

* "Forecasts" (actually a series of 1-step forecasts) in levels:
smpl 1961:11 1995:8
steps 1
# armaxsmall armasmall_fore
smpl

smpl 1990:1 1995:8
graph(key=below,header='Interest Rate One-Step Forecasts, Actuals') 3
# rate
# arma_fore
# armasmall_fore

* Now compare actual multi-step forecasts:
smpl 1994:1 1995:8
forecast 1
# armaxsmall armasmall_f
forecat 1
# armax arma_f
graph(key=below,header='Interest Rate 20-Step Forecasts, Actuals') 3
# rate
# arma_f
# armasmall_f

* Theil U's:
Display "Theil U's for P&R Model:"
theil(setup) 1 24 1995:8
# armax
do time=1961:1,1995:8
 theil time
end
theil(dump)

Display "Theil U's for Small Model:"
theil(setup) 1 24 1995:8
# armaxsmall
do time=1961:1,1995:8
 theil time
end
theil(dump)

*********************************************************************
*** Try capping estimation at 1985:12 and producing Theil U's for
***  1986:1 through 1995:8. Will need to generate the out of sample
***  residuals manually using STEPS.
**********************************************************************

display "Limiting estimation range to 1990:12"
compute regend = 1990:12
smpl * regend

* Estimate P&R's ARMA model with regressors as inputs:
display "P&R's ARMAX BOXJENK:"
nlpar(sub=100)
boxjenk(inputs=3,ar=8,ma=2,constant,iters=100,define=armax) rate / armaxres
# ip 0 0 0
# grm2 0 0 0
# grppi 0 0 1

display "Smaller ARMAX model:"
boxjenk(inputs=3,ar=2,ma=1,constant,iters=100,define=armaxsmall) rate / armaxsres
# ip 0 0 0
# grm2 0 0 0
# grppi 0 0 1

* Generate out of sample residuals:
smpl regend+1 1995:8

steps 1
# armax armaxfstep * armaxres
steps 1
# armaxsmall armaxfsteps * armaxsres

* Theil U's:
Display "Theil U's for P&R Model, 1991:1-- 1995:8:"
compute start = regend+1, end = 1995:8

theil(setup) 1 24 end
# armax
do time=start,end
 theil time
end
theil(dump)

Display "Theil U's for Small Model, 1991:1-- 1995:8:"
theil(setup) 1 24 end
# armaxsmall
do time=start,end
 theil time
end
theil(dump)