IRF for STAR models

[Javi]

Dear Rats users:

I am looking for a code to estimate IRF for STAR models. Is there anybody who has a code to estimate it?

Thanks in advance.



Best regards,

[TomDoan]

Because of the non-linearities, STAR models don't really have an IRF like a VAR does; the response to a shock depends on the
current value and doesn't scale up linearly. The "eventual forecasting function" is commonly used to give some idea of the
long-term behavior. An example of that is in the Terasvirta 1994 replication file LYNX.PRG.

If you want something closer to an actual IRF, you can do the eventual forecast as shown in that example, then repeat with a
SHOCKS option on the FORECAST. The difference between the two would be the analog of an IRF. Note, however, that this can be quite sensitive to the size of the shock, that is, (unlike linear models) the response to a "1.0" shock won't be half the response to a "2.0" shock.

This adds the calculation of the responses to shocks of size 1.0 and 3.0 to the end of the lynx example. You can see how the behavior is quite far from being scales of each other.

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

[luxu1983]

forecast(model=starmodel,steps=100,from=1935:1,shocks=1.0)
what do the shocks stand for? can stand for standard deviation

[TomDoan]

No. That's 1.0 and 3.0. If you want one standard deviation, use SHOCKS=||sqrt(%seesq)||.

[luxu1983]

dear sir
can this method (using forcecast)calculateing the irf extend to multivariable nonlinear model ?

[TomDoan]

It's not "the IRF", but yes.

[luxu1983]

It's not "the IRF", but yes.

dear tom
can i call it GIRF?
i find G Koop, MH Pesaran, and SM Potter(1996) extend linear irf to the nonlinear case by defining a generalized impulse response function that is very similar to your method
but in their method includes " simulation"

may you provide the example that is nonlinear case using the generalized impulse response function ?

[TomDoan]

What they're doing is taking the base forecast, then, instead of doing one shock of a certain size and taking the difference as an estimate of an IRF, they do one shock and then draw a realization of the error process beyond that and average. Since

1. the average isn't necessarily a good measure of central tendency in a highly non-linear function,
2. the average depends upon the distribution from which the other shocks are drawn,
3. the results are dependent upon the initial values

the whole thing still makes little sense. If you want to pursue it, read chapter 13 of the User's Guide, particularly sections 13.1 and 13.2. Those will explain how to get the averages. A single draw would be done with something like:

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

where delta is the initial shock and mystddev is the (mean zero) standard deviation for the Normal density for the future shocks.

[luxu1983]

dear Tom
thank you very very much
i modify it follow you suggestion adding shocks
but it doesn't work well

group starmodel star>>sforecast
forecast(model=starmodel,steps=100,from=1935:1)
set base 1935:1 1935:1+99 = sforecast
set shocks = %if(t==1935:1,0.1,%ran(1))
forecast(model=starmodel,steps=100,from=1935:1,path)
# shocks
set irf1 1935:1 1935:1+99 = sforecast-base
graph(header="Responses to shocks",key=upright)
# irf1

[TomDoan]

dear Tom
thank you very very much
i modify it follow you suggestion adding shocks
but it doesn't work well

group starmodel star>>sforecast
forecast(model=starmodel,steps=100,from=1935:1)
set base 1935:1 1935:1+99 = sforecast
set shocks = %if(t==1935:1,0.1,%ran(1))
forecast(model=starmodel,steps=100,from=1935:1,path)
# shocks
set irf1 1935:1 1935:1+99 = sforecast-base
graph(header="Responses to shocks",key=upright)
# irf1

Sorry. The set shocks should have read

Code: Select all

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

I've also corrected that in the original post. However, %ran(1) is way out of bounds for this series. As I mentioned in point 2 above, the average across draws will depend upon the distribution you use there. When you consistently hit the process with shocks which are something like six standard deviations, you don't know what's going to happen.

[guo]

Dear Tom,

Thank you so much for your previous advice and I really appreciate your suggestions.

I am doing ESTAR model, and I want to get nonlinear "IRF" in my analysis.

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



However, the error emerges,

Code: Select all

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



Under this circumstance, what should I do to get the "IRF", because I have some exogenous variables (dlneer and dlppi) in the ESTAR model.

Thank you very much for your precious time and attention.

Your help is greatly appreciated.

Kind Regards

guo

[TomDoan]

Without the complete program and data set, I can't tell much about it.