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