*
* TARMODELS.RPF
* Example of estimation of a STAR Model
* RATS User's Guide, Example from Section 11.5.
*
* This is based on example 3 from Terasvirta(1994), "Specification,
* Estimation and Evaluation of Smooth Transition Autoregressive Models",
* JASA, vol 89, pp 208-218.
*
cal(a) 1821
open data lynx.dat
data(org=cols) 1821:1 1934:1 lynx
set x = log(lynx)/log(10)
*
* The YuleLags procedure does a quick, efficient examination of a range
* of AR models for stationary data. Because this is for a standard AR,
* with no threshold effects, it might not give a good estimate of the
* lag length in the presence of a threshold.
*
@yulelags(max=20) x
*
* Standard practice is to pick the delay by choosing the one which gives
* the most significant test statistic. The tests are done on the
* demeaned data using 11 lags and up to a 9 period delay.
*
diff(center) x / xc
do d=1,9
@StarTest(p=11,d=d) xc
end do d
*
* This initially fixes the c and gamma at the mean and two times the
* reciprocal of the standard error, respectively. (The logistic exponent
* is rescaled by 1/sigma as shown in the article), so gamma is, with
* that parameterization, put at 2. This gives a fairly sharp split
* between larger and smaller values, so two branches should have
* decidedly different coefficients if there is evidence of STAR. With c
* and gamma fixed, NLLS is actually just OLS.
*
stats x
*
* In practice, you would use the following. The author used a rounded value.
*
compute scalef=1.0/sqrt(%variance)
nonlin(parmset=starparms) gamma c
frml glstar = %logistic(scalef*gamma*(x{3}-c),1.0)
compute c=%mean,gamma=2.0
*
* These are the two branches. Putting these in as linear equations makes
* specification more flexible---you can easily change the lag lists and
* the rest goes through as is.
*
equation standard x
# constant x{1 to 11}
equation transit x
# constant x{1 to 11}
*
* Convert the linear equations into FRML's with phi1 and phi2 as the
* coefficient vectors. Put them together with the transition function to
* make the STAR formula.
*
frml(equation=standard,vector=phi1) phi1f
frml(equation=transit ,vector=phi2) phi2f
frml star x = g=glstar,phi1f+g*phi2f
*
nonlin(parmset=regparms) phi1 phi2
nonlin(parmset=starparms) gamma c
*
* Estimate the model with the STARPARMS left out.
*
compute phi1=%const(0.0),phi2=%const(0.0)
nlls(parmset=regparms,frml=star) x
*
* Based upon the initial results, the standard equation is trimmed to
* just x{1} and transit to x{2 3 4 10 11}.
*
equation standard x
# x{1}
equation transit x
# x{2 3 4 10 11}
frml(equation=standard,vector=phi1) phi1f
frml(equation=transit ,vector=phi2) phi2f
*
* The new specification is re-estimated with just the regression
* parameters, then including the STAR parameters.
*
* Because the PHI1 and PHI2 get re-dimensioned when we use them for
* the new equations, the previous values are lost (and should be,
* since the role of each element has likely changed).
*
compute phi1=%const(0.0),phi2=%const(0.0)
nlls(parmset=regparms,frml=star) x
nlls(parmset=regparms+starparms,frml=star) x
*
* Analysis of roots of the polynomials
*
compute %eqnsetcoeffs(standard,phi1)
compute %eqnsetcoeffs(transit,phi2)
compute upper=standard+transit
*
* The equation <> is created by adding the right hand sides of
* the standard and transit polynomials.
*
@LagPolyRoots(title="Above Threshold Polynomial") %eqnlagpoly(upper,x)
@LagPolyRoots(title="Below Threshold Polynomial") %eqnlagpoly(standard,x)
*
* To forecast the non-linear equation, we need to create a MODEL with
* the one formula. This uses SFORECAST for the results for all
* forecasting calculations, so we need to pull information out that we
* need as we do forecasts.
*
group starmodel star>>sforecast
*
* Compute eventual forecasting function starting in 1925:1.
*
forecast(model=starmodel,steps=40,from=1925:1)
set eforecast 1925:1 1964:1 = sforecast
*
* Computes forecasts for 40 steps beginning in 1925 by taking
* the mean of simulated values over 10000 replications.
*
set meanf 1925:1 1964:1 = 0.0
compute nreps=10000
do reps=1,nreps
simulate(model=starmodel,from=1925:1,to=1964:1)
set meanf 1925:1 1964:1 = meanf+sforecast
end do reps
set meanf 1925:1 1964:1 = meanf/nreps
*
* Graphs the eventual forecast, the mean forecast and the last few years
* of actual data. While the process has a fairly regular cycle of
* roughly eight years, shocks will shift the timing enough that even 16
* years out, the probability is fairly high that it will be in a
* somewhat difference phase than it is today. The mean forecasts thus
* eventually converge to the series mean.
*
graph(footer="Forecasts of Lynx Data with LSTAR Model",$
key=loright,klabels=||"Eventual Forecast","Mean Forecast","Actual"||) 3
# eforecast
# meanf
# x 1920:1 1934:1