The RATS Software Forum

Code: Select all

*
* 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)
diff(center) x / xc
*
* The YuleLags procedure does a quick, efficient examination of a range of AR
* models for stationary data.
*
@yulelags(max=20) x
*
* The tests shown are done on the demeaned data
*
source startest.src
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
compute scalef=1.0/sqrt(%variance)
nonlin(parmset=starparms) gamma c
frml flstar = %logistic(1.8*gamma*(x{3}-c),1.0)
compute c=%mean,gamma=2.0
equation standard x
# constant x{1 to 11}
equation transit x
# constant x{1 to 11}
*
frml(equation=standard,vector=phi1) phi1f
frml(equation=transit ,vector=phi2) phi2f
frml star x = f=flstar,phi1f+f*phi2f
*
nonlin(parmset=regparms) phi1 phi2
nonlin(parmset=starparms) gamma c
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} (The article shows lag 9 rather than 10, but this
* specification fits quite a bit better). This is now estimated with all the
* parameters.
*
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
nlls(parmset=regparms,frml=star) x
nlls(parmset=regparms+starparms,frml=star) x
*
* This will do the model shown in the article
*
equation standard x
# x{1}
equation transit x
# x{2 3 4 9 11}
frml(equation=standard,vector=phi1) phi1f
frml(equation=transit ,vector=phi2) phi2f
nlls(parmset=regparms,frml=star) x
nlls(parmset=regparms+starparms,frml=star) x
*
stats %resids
*
* Analysis of roots of the polynomials
*
associate standard phi1
associate transit  phi2
compute upper=standard+transit
*
* The equation <<upper>> is created by adding the right hand sides of the standard
* and transit polynomials.
*
@LagPolyRoots %eqnlagpoly(upper,x)
@LagPolyRoots %eqnlagpoly(standard,x)
*
* Compute eventual forecasting function
*
group starmodel star>>sforecast
forecast(model=starmodel,steps=100,from=1935:1)
graph 2
# sforecast
# x 1920:1 1934:1
*
set base 1935:1 1935:1+99 = sforecast
forecast(model=starmodel,steps=100,from=1935:1,shocks=1.0)
set irf1 1935:1 1935:1+99 = sforecast-base
forecast(model=starmodel,steps=100,from=1935:1,shocks=3.0)
set irf3 1935:1 1935:1+99 = sforecast-base
graph(header="Responses to shocks",key=upright,klabels=||"1.0 Shock","3.0 Shock"||) 2
# irf1
# irf3

Code: Select all

set shocks 1935:1 1935:1+99 = %if(t==1935:1,delta,%ran(mystddev))
forecast(model=starmodel,steps=100,from=1935:1,paths)
# shocks
set irf1 1935:1 1935:1+99 = sforecast-base

Code: Select all

set shocks 1935:1 1935:1+99 = %if(t==1935:1,0.1,%ran(1))

Code: Select all

* This will do the model shown in the article
*
equation standard dlEX_CA
# Constant dlneer{4 8 12} dlppi{3 10 11} dlEX_CA{1 2 3 4 12}

equation transit dlEX_CA
# Constant dlneer{0 8 12} dlppi{1 2 5 6 7} dlEX_CA{1 2 3 4 9}

frml(equation=standard,vector=phi1) phi1f
frml(equation=transit ,vector=phi2) phi2f
nlls(parmset=regparms,frml=star) dlEX_CA
nlls(parmset=regparms+starparms,frml=star) dlEX_CA
*
stats %resids
*
* Analysis of roots of the polynomials
*
associate standard phi1
associate transit  phi2
compute upper=standard+transit
*
* The equation <<upper>> is created by adding the right hand sides of the standard
* and transit polynomials.
*
@LagPolyRoots %eqnlagpoly(upper,dlEX_CA)
@LagPolyRoots %eqnlagpoly(standard,dlEX_CA)
**************************************************************
* Compute eventual forecasting function
group starmodel star>>sforecast
forecast(model=starmodel,steps=100,from=2012:4)
graph 2
# sforecast
# dlEX_CA 2003:1 2012:3

set base 2012:4 2012:4+99 = sforecast
forecast(model=starmodel,steps=100,from=2012:4,shocks=1.0)
set irf1 2012:4 2012:4+99 = sforecast-base
forecast(model=starmodel,steps=100,from=2012:4,shocks=3.0)
set irf3 2012:4 2012:4+99 = sforecast-base
graph(header="Responses to shocks",key=upright,klabels=||"1.0 Shock","3.0 Shock"||) 2
# irf1
# irf3

Code: Select all

## SR14. Empty Range on Series SFORECAST.
## I2. Expected Instruction Here
>>>>#<<<<
## SR14. Empty Range on Series IRF1.
## I2. Expected Instruction Here
>>>>#<<<<