* * 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)