*
* Application from section 4.5, pp 194-198
* (Data on web site shorter than described)
*
open data q-unemrate.txt
calendar(q) 1948
data(format=free,org=columns) 1948:1 1991:2 unemp
*
* Estimate the seasonal ARIMA model
*
boxjenk(maxl,diffs=1,ar=1,sar=1,sma=1) unemp
*
diff unemp / dx
dofor delay = 1 2 3 4
   set z = dx{delay}
   @tsaytest(threshold=z) dx
   # dx{1 to 5}
end dofor delay
*
* Ori-F test
*
dec symm[series] m(5,5)
do i=1,5
   do j=1,i
      set m(i,j) = dx{i}*dx{j}
   end do j
end do i
linreg dx
# dx{1 to 5} m
exclude(title="Ori-F test")
# m
*
* LST modification
*
dec vect[series] lst(5)
do i=1,5
   set lst(i) = dx{i}**3
end do i
linreg dx
# dx{1 to 5} m lst
exclude(title="LST test")
# m lst
*
* Estimation of the threshold model (4.53)
*
linreg(smpl=dx{2}<=.1) dx
# constant dx{1 2}
@lagpolyroots %eqnlagpoly(0,dx)
linreg(smpl=dx{2}>.1) dx
# constant dx{1 2}
@lagpolyroots %eqnlagpoly(0,dx)
*
* Markov-Switching model
*
* This is much simpler than the Hamilton switching model given in the RATS manual.
* The complication there is that the state of the economy at each lag matters
* because the autoregression isn't on the observed data, but on the difference
* between the observed data and the unobserved state-dependent mean.
*
source markov.src
*
compute nstates=2
dec vector pstar(nstates)
dec rect p(nstates-1,nstates)
dec series[vect] pt_t pt_t1 psmooth
gset pt_t    = %zeros(nstates,1)
gset pt_t1   = %zeros(nstates,1)
gset psmooth = %zeros(nstates,1)
nonlin(parmset=msparms) p
*
* Application-specific
*
compute nlags=2
dec rect phi(nlags+1,nstates)
dec vect sigma(nstates)
nonlin(parmset=apparms) phi sigma
*
* Note that the states aren't globally identified. We would get exactly the same
* likelihood if state 1 were contraction and 2 were expansion. The choice of
* initial values for phi and the optimization method can alter whether, in the
* process of estimation, the states "trade places." In this case, we start with
* the OLS coefficients for state 1, and their sign-flipped values for state 2.
*
linreg dx
# constant dx{1 to nlags}
ewise phi(i,j)=%if(j==1,%beta(j),-%beta(j))
compute sigma=%fill(2,1,sqrt(%seesq))
*
* ARStateF returns the vector of likelihoods of the two states at period "time"
* given the parameters of the autoregression.
*
function ARStateF time
type vector ARStateF
type integer time
*
local integer i j
local real u
local vect reg
*
dim reg(nstates)
*
do i=1,nstates
   compute u=dx(time)-phi(1,i)
   do j=1,nlags
      compute u=u-phi(j+1,i)*dx(time-j)
   end do j
   compute reg(i)=u
end do i
*
dim ARStateF(nstates)
ewise ARStateF(i)=%density(reg(i)/sigma(i))/sigma(i)
end
*
* The log likelihood function recursively calculates the vector PSTAR of estimated
* state probabilities. The likelihoods in the states are computed by the ARStateF
* function, and then %mcstate and %msupdate update the PSTAR vector and compute
* the likelihood. We use a START formula on MAXIMIZE to set PSTAR to the ergodic
* probabilities.
*
compute p=||.90,.20||
*
frml logl = f=ARStateF(t),pt_t1=%mcstate(p,pstar),pt_t=pstar=%msupdate(f,pt_t1,fpt),log(fpt)
maximize(start=(pstar=%mcergodic(p)),parmset=apparms+msparms,$
 method=bhhh,pmethod=simplex,piters=5) logl nlags+2 *
*
do i=1,2
   disp "Conditional Mean for State" i phi(1,i)/(1-phi(2,i)-phi(3,i))
end do i