*
* ARIMA.RPF
* RATS Version 8, User's Guide, Example from Section 9.4
* From Enders, Applied Econometric Time Series, 3rd edition
*
open data quarterly.xls
calendar(q) 1960:1
data(format=xls,org=columns) / tbill r10
graph(header="T-Bill and 10-year Bond Rates",key=upleft) 2
# tbill
# r10
*
* Compute spread and first difference of spread:
*
set spread = r10 - tbill
diff spread / dspread

spgraph(vfields=2,footer="FIGURE 2.5 Time Path of Interest Rate Spread")
 graph(header="Panel (a): The interest rate spread")
 # spread
 graph(header="Panel (b): First difference of the spread")
 # dspread
spgraph(done)
*
* Compute and graph autocorrelations of the spread itself
*
corr(results=corrs,partial=pcorrs,number=12) spread
graph(footer="FIGURE 2.6 ACF and PACF of the Spread",key=below,$
   style=bar,nodates,min=-1.0,max=1.0,number=0) 2
# corrs
# pcorrs
*
* Compute and graph autocorrelations for the first difference
*
corr(results=dcorrs,partial=dpcorrs,number=12) dspread
graph(header="Correlations of the first difference of spread",key=below,$
   style=bar,nodates,min=-1.0,max=1.0,number=0) 2
# dcorrs
# dpcorrs
*
* Estimating various candidate models. To ensure comparability, they are
* all run over the period from 1961:4 to the end of data. 1961:4 is the
* earliest period for which an AR(7) model (the longest AR considered)
* can be run using conditional least squares.
*
boxjenk(constant,ar=7) spread 1961:4 *
correlate(results=rescorrs,number=12,span=4,qstats,dfc=%narma) %resids
graph(header="AR(7)",style=bar,nodates,min=-1.0,max=1.0,number=1)
# rescorrs 2 *
*
@regcorrs(dfc=%narma,number=20,qstats,report,$
   method=burg,title="AR(7) model diagnostics")
*
* AR(2)
*
boxjenk(constant,ar=2) spread 1961:4 2008:1
@regcorrs(dfc=%narma,number=20,qstats,report,$
   method=burg,title="AR(2) model diagnostics")
*
* AR({1,2,7})
*
boxjenk(constant,ar=||1,2,7||) spread 1961:4 2008:1
@regcorrs(dfc=%narma,number=20,qstats,report,$
   method=burg,title="AR({1,2,7}) model diagnostics")
*
* ARMA(1,1)
*
boxjenk(constant,ar=1,ma=1) spread 1961:4 2008:1
@regcorrs(dfc=%narma,number=20,qstats,report,$
   method=burg,title="ARMA(1,1) model diagnostics")
*
* ARMA(2,(1,7))
*
boxjenk(constant,ar=2,ma=||1,7||) spread 1961:4 2008:1
@regcorrs(dfc=%narma,number=20,qstats,report,$
   method=burg,title="ARMA(2,(1,7)) model diagnostics")
*
* Estimate two candidate models and compare forecasting performance:
*
* Use all available data through 1995:3 to estimate each model (use *
* for "start" parameter to let RATS determine start date)
*
boxjenk(constant,define=ar7eq,ar=7)                spread * 1995:3
boxjenk(constant,define=ar2ma17eq,ar=2,ma=||1,7||) spread * 1995:3
*
* Compute and print one-step forecast (and corresponding actual value)
* for 1995:4.
*
uforecast(equation=ar7eq,static,print) onestep_ar7 1995:4 1995:4
uforecast(equation=ar2ma17eq,static,print) onestep_ar2ma17 1995:4 1995:4
*
* Now, compute and graph dynamic forecasts for 1995:4 through 2008:1.
*
uforecast(equation=ar7eq,print) fore_ar7 1995:4 2008:1
uforecast(equation=ar2ma17eq,print) fore_ar2ma17 1995:4 2008:1
graph(header="Forecasts vs Actuals",key=upleft) 3
# spread 1995:1 *
# fore_ar7
# fore_ar2ma17 / 5
*
* Next, compute one-step forecasts for all 50 periods, re-estimating at
* each period.
*
do time=1995:4,2008:1
   boxjenk(noprint,constant,define=ar7eq,ar=7) spread * time-1
   boxjenk(noprint,constant,define=ar2ma17eq,ar=2,ma=||1,7||) spread * time-1
   uforecast(equation=ar7eq,static) forecast_ar7 time time
   uforecast(equation=ar2ma17eq,static) forecast_ar2ma17 time time
end do
*
* Graph the forecast with the actuals.
*
graph(header="Forecasts vs Actuals",key=upleft) 3
# spread 1995:1 *
# forecast_ar7
# forecast_ar2ma17eq / 5
*
* Compute means, variances:
*
set e1 = forecast_ar7 - spread
set e2 = forecast_ar2ma17 - spread
stats e1
stats e2
*
* Using UForeErrors procedure:
*
@uforeerrors spread forecast_ar7
@uforeerrors spread forecast_ar2ma17
*
* Check for bias:
*
linreg spread
# constant forecast_ar7
test(all)
# 0.0 1.0
linreg spread
# constant forecast_ar2ma17
test(all)
# 0.0 1.0
*
* Granger-Newbold test
* Computing it directly:
*
set x = e1 + e2
set z = e1 - e2
cmom(corr)
# x z
compute corrcoef = %cmom(2,1)
compute gnstat = corrcoef/sqrt((1-corrcoef^2)/(%nobs-1))
cdf(noprint) ttest gnstat %nobs-1
display "G-N statistic = " *.### gnstat "Signif. Level = " #.##### %signif/2
display
*
* Using the GNewbold procedure:
*
@gnewbold spread forecast_ar7 forecast_ar2ma17
*
* Diebold-Mariano statistic on 4th power
*
set d = e1^4 - e2^4
statistics(noprint) d
compute dmstat = sqrt((%nobs-1)/%variance)*%mean
cdf(noprint) ttest dmstat %nobs-1
display "D-M statistic on 4th Power= " *.### dmstat "Signif. Level = " #.##### %signif/2
display
*
* RATS also includes a procedure for doing D-M tests based on MSE or MAE:
*
@dmariano(criterion=mse) spread forecast_ar7 forecast_ar2ma17
@dmariano(criterion=mae) spread forecast_ar7 forecast_ar2ma17
*
* Using THEIL to generate forecast performance statistics for dynamic
* forecasts up to 8 steps ahead.
*
boxjenk(noprint,constant,ar=7,define=ar7eq) spread * 1995:3
theil(setup,steps=8,to=2008:1)
# ar7eq
do time=1995:4,2008:1
   theil time
   boxjenk(noprint,constant,ar=7,define=ar7eq) spread * time
end do time
theil(dump,title="Theil U results for AR(7) Model")
*
boxjenk(noprint,constant,define=ar2ma17eq,ar=2,ma=||1,7||) spread * 1995:3
theil(setup,steps=8,to=2008:1)
# ar7eq
do time=1995:4,2008:1
   theil time
   boxjenk(noprint,constant,define=ar2ma17eq,ar=2,ma=||1,7||) spread * time
end do time
theil(dump,title="Theil U results for ARMA(2,(1,7)) Model")
*
* Using the procedures
*
* We'll read in two additional series used to demonstrate additional
* features of some of the procedures:
*
open data quarterly.xls
calendar(q) 1960:1
data(format=xls,org=columns) / tbill r10 ppinsa m1nsa
set spread = r10 - tbill
*
* Compare log, square root, and no transformation:
*
@bjtrans ppinsa
@bjtrans m1nsa
*
* Enders suggests the possibility of transforming a version of SPREAD
* adjusted to have positive values for all entries. Here's how you could
* do that:
*
set spreadshift = spread + 2
@bjtrans spreadshift
*
* BJIDENT plots for SPREAD:
*
@bjident(diffs=1) spread
*
* BJIDENT plots for the log of M1NSA (via the TRANS=LOG option):
*
@bjident(diffs=1,sdiffs=1,trans=log,number=24,report,qstat) m1nsa
*
* Using BJFORE, first on the AR(7) model for SPREAD:
*
@bjfore(ars=7,constant) spread 1995:4 2008:1 fore
graph(header="Forecasts vs Actuals") 2
# fore
# spread 1993:1 *
@regcorrs(dfc=%narma,method=burg)
*
* And on an AR(1)SMA(1) model for M1NSA:
*
@bjfore(ars=1,smas=1,trans=log,diffs=1,sdiffs=1,constant) m1nsa 2008:2+1 2008:2+12 m1fore m1resids
graph(header="M1NSA Out of Sample Forecasts and Actuals") 2
# m1nsa 2000:1 *
# m1fore
*
* Do an "autofit" allowing up to 7 AR and 7 MA, using the Schwarz criterion
*
@bjautofit(constant,pmax=7,qmax=7,crit=sbc) spread

boxjenk(constant,define=ar2ma1eq,ar=%%autop,ma=%%autoq) spread
@regcorrs(dfc=%narma,method=burg)

do time=1995:4,2008:1
   boxjenk(noprint,constant,define=ar2ma1eq,ar=2,ma=1) spread * time-1
   uforecast(equation=ar2ma1eq,static) forecast_ar2ma1 time time
end do
@uforeerrors spread forecast_ar2ma1