*
* Hamilton Markov chain model for US GNP from page 697
*
cal(q) 1951:1
open data gnpdata.prn
data(format=prn,org=columns) 1951:1 1984:4
*
set g = 100*log(gnp/gnp{1})
graph
# g
*
* Hamilton's GNP growth model has two possible states at each period (interpreted
* as expansion and contraction). With one lag plus the current, there are 4
* possible combinations of states when analyzing a given time period.
*
compute nlags=4
compute nstates=2**(nlags+1)
dec vect mu(2) phi(nlags)
dec vect tp(5)
*
nonlin mu phi p11 p22 sigma
*
* To speed matters up, create a lookup table which gives the mapping from the
* coding for the states into the binary choice for each lag. The current state is
* 1 for 1-16 and 2 for 17-32. The first lag is 1 for 1-8 and 17-24 and 2 for 9-16
* and 25-32, etc.
*
dec rect[integer] lagstate(nstates,nlags+1)
ewise lagstate(i,j)=1+%mod((i-1)/2**(nlags+1-j),2)
*
* We also create a mapping from state to state which will pick out the transition
* probability. transit will be a number from 1 to 5. 1 is the most common value,
* which represents a zero transition probability. For instance, if we're in
* (1,1,2,2,1), the only states which can be obtained are {1,1,1,2,2} and
* {2,1,1,2,2}. The transition probability to the first of these would be p11, and
* for the second it would be 1-p11. Given our coding, it's easy to tell whether a
* state can be obtained from another, since the final four spots of the new state
* have to be the same as the lead four spots of the old state.
*
dec rect[integer] transit(nstates,nstates)
ewise transit(i,j)=fix(%if(%mod(i-1,2**nlags)==(j-1)/2,2*lagstate(i,1)+lagstate(j,1)-1,1))
*
dec vect reg(nstates) f(nstates) p(nstates) fp(nstates) pstar(nstates)
*
* The mask series is used to sum up the probability of being in state 1 at any
* point in time. The prob1 series will hold the result of this calculation.
*
dec vect mask(nstates)
ewise mask(i)=(i<=nstates/2)
set prob1 = 0.0
*
* We create a special purpose function to do a one period update of the Markov
* chain. We can't use one of the standard functions described in the RATS manual
* because we don't have a full 32x32 transition matrix, instead using the lookup
* table into a smaller one.
*
function MarkovProb time
type real MarkovProb
type integer time
*
local integer i j
local real u
local real sp fpt

do i=1,nstates
   compute u=g(time)-mu(lagstate(i,1))
   do j=1,nlags
      compute u=u-phi(j)*(g(time-j)-mu(lagstate(i,j+1)))
   end do j
   compute reg(i)=u
end do i
*
ewise f(i)=%density(reg(i)/sigma)/sigma
do i=1,nstates
   compute sp=0.0
   do j=1,nstates
      compute sp=sp+tp(transit(i,j))*pstar(j)
   end do j
   compute p(i)=sp
end do i
compute fp=f.*p
compute fpt=%sum(fp)
compute pstar=fp*(1.0/fpt)
compute prob1(time)=%dot(mask,pstar)
compute MarkovProb=fpt
end
*
frml ggrowth = log(MarkovProb(t))
*
* pstar (the probability distribution of the states at the previous period) needs
* to be initialized with each function evaluation, as it gets updated at each
* period. We also copy the transition probabilities into the tp vector. The
* initial distribution is set equal to the ergodic probabilities of the states,
* given the current parameters. The solution method is shown on page 684 of
* Hamilton.
*
* The "A" matrix from 684 has 1-the transition probs in the top NxN, with a row
* of 1's below that.
*
function MarkovInit
type vector MarkovInit
*
local rect a
local integer i j
*
compute tp=||0.0,p11,1-p22,1-p11,p22||
dim a(nstates+1,nstates)
ewise a(i,j)=%if(i>nstates,1.0,(i==j)-tp(transit(i,j)))
compute MarkovInit=%xcol(inv(tr(a)*a)*tr(a),nstates+1)
end
*
* 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 mu and the optimization method can alter whether, in the
* process of estimation, the states "trade places."
*
compute p11=.85,p22=.70
compute mu(1)=1.0,mu(2)=0.0
linreg g
# constant g{1 to 4}
compute phi=%xsubvec(%beta,2,5)
compute sigma=%seesq
*
maximize(start=(pstar=MarkovInit()),method=bfgs) ggrowth nlags+2 *
*
* To create the shading marking the recessions, create a dummy series which is 1
* when the recessq series is 1, and 0 otherwise. (recessq is 1 for NBER
* recessions and -1 for expansions).
*
set contract = recessq==1
*
spgraph(vfields=2)
graph(header="Quarterly Growth Rate of US GNP",shade=contract)
# g %regstart() %regend()
set prob2 = 1-prob1
graph(style=polygon,header="Probability of Economy Being in Contraction",shade=contract)
# prob2 %regstart() %regend()
spgraph(done)