cal 1952 1 4
all 1998:4
open data drisampl.rat
data(format=rats) / gdpq
set g = 100*(gdpq/gdpq{1}-1)

* EM Algorithmn

* Starting Values for Estimated paramters from the summary statistics of data

statistics g


* note that in the way we are setting at truth


compute a1=%mean+1
compute a2=%mean-1


* get initial estimate of phi from AR(1) coefficient

linreg g * * coeff
# g{1}

compute phi=%beta(1)

* Get initial estimate of sigma from variance of residuals

compute sigma=%seesq**.5


* Starting probabilities, states labeled ij for s_t=j,s_t-1=i
* states are 11 - st=1 st-1=1
*            21 - st=1 st-1=2
*            12 - st=2 st-1=1
*            22 - st=2 st-1=2
* So, we know that 1-p1 = p2 and 1-p3=p4

comp p11 = .7
comp p21 = 1-p11
comp p22 = .7
comp p12 = 1-p22


* ACTUAL CONVERGENCE ALGORITHM

compute  conv1=0.0


while conv1 < 5 {

* evaluate unconditional probabilities of each regime
* label pr11-pr22
* evaluated as P{s_t|s_t-1}*P{s_t=1}

* State 11
comp pr11 = p11*(1-p22)/(2-p11-p22)

* State 12
comp pr12 = p12*(1-p11)/(2-p11-p22)

* State 21
comp pr21 = p21*(1-p22)/(2-p11-p22)

* State 22
comp pr22 = p22*(1-p11)/(2-p11-p22)


*evaluate the density of g|state for state 11,21,12,22, ie the likelihood that g was
*generated by the dgp for that state

set f11 = %density(((g-a1)-phi*(g{1}-a1))/sigma)
set f21 = %density(((g-a1)-phi*(g{1}-a2))/sigma)
set f12 = %density(((g-a2)-phi*(g{1}-a1))/sigma)
set f22 = %density(((g-a2)-phi*(g{1}-a2))/sigma)


*evaluate numerator of P(st|yt) for states 1-4 ie [22.3.7]

set rp11 / =f11*pr11
set rp21 / =f21*pr21
set rp12 / =f12*pr12
set rp22 / =f22*pr22


*evaluate denominator of [22.3.7]

set fy = rp11 + rp12 + rp21 + rp22

* Evaluate expression in 22.3.7 in Hamilton for P{st=j|y;param}

set hp11 = rp11/fy
set hp21 = rp21/fy
set hp12 = rp12/fy
set hp22 = rp22/fy

*set test = hp11 + hp12 + hp21 + hp22
*statistics test
* Using the expression in 22.3.10, the estimated probabilities should be the mean of the likelihoods of each obs
* being from a given state, and the sums of the state probabilities are stored as sum1-sum4 for later use.
* notationally, p1n represents the new estimate of p1, etc.  Updating done at the end.


statistics(print) hp11
compute pr11n =%mean
*compute pr21n= 1- pr11n
comp sum11=%nobs*%mean

statistics(print) hp21
comp sum21=%nobs*%mean
compute pr21n=%mean

statistics(print) hp12
comp sum12=%nobs*%mean
compute pr12n=%mean

statistics(print) hp22
comp sum22=%nobs*%mean
compute pr22n= %mean

* Now, using new regime estimates, work backwards to get new estimates for transition
* matrix
* State 11
comp p11n = pr11n/(1-p22)*(2-p11-p22)


* State 21
comp p21n = 1-p11n

* State 22
comp p22n = pr22n/(1-p11)*(2-p11-p22)

* State 12
comp p12n = 1-p22n


dis p11n p12n
dis p21n p22n



*Mean Estimator
set m11 = hp11*(g-phi*(g{1}-a1))
set m21 = hp21*(g-phi*(g{1}-a2))
set m12 = hp12*(g-phi*(g{1}-a2))
set m22 = hp22*(g-a2-phi*(g{1}-a2))

statistics m11
comp s11=%mean*%nobs
statistics m21
comp s21=%mean*%nobs
statistics m12
comp s12=%mean*%nobs
statistics m22
comp s22=%mean*%nobs

comp a1n = (pr11/(pr11+pr21))*(s11/sum11)+(pr21/(pr11+pr21))*(s21/sum21)
comp a2n = (pr22/(pr22+pr12))*(s22/sum22)+(pr12/(pr22+pr12))*(s12/sum12)
dis a1n a2n


* Variance Estimator
* Weighted average of squared residual for each state
set res11 = hp11*(g-a1-phi*(g{1}-a1))**2
set res21 = hp21*(g-a1-phi*(g{1}-a2))**2
set res12 = hp12*(g-a2-phi*(g{1}-a2))**2
set res22 = hp22*(g-a2-phi*(g{1}-a2))**2

statistics res11
comp v11=%mean*%nobs
statistics res21
comp v21=%mean*%nobs
statistics res12
comp v12=%mean*%nobs
statistics res22
comp v22=%mean*%nobs
comp var = pr11*(v11/sum11)+pr21*(v21/sum21)+pr22*(v22/sum22)+pr12*(v12/sum12)
comp sigman = var**(0.5)

* These would be appropriate if there was a different variance for regime 1 and 2
* comp var1n = (pr11/(pr11+pr21))*(v11/sum11)+(pr21/(pr11+pr21))*(v21/sum21)
* comp var2n = (pr22/(pr22+pr12))*(v22/sum22)+(pr12/(pr22+pr12))*(v12/sum12)


* Phi Estimator
set reg1 = g-a1
set reg2 = g-a2

linreg reg1 * * coeff
# reg1{1}
compute phi11=%beta(1)

linreg reg1 * * coeff
# reg2{1}
compute phi21=%beta(1)

linreg reg2 * * coeff
# reg1{1}
compute phi12=%beta(1)

linreg reg2 * * coeff
# reg2{1}
compute phi22=%beta(1)

comp phin=pr11*phi11+pr12*phi12+pr21*phi21+pr22*phi22


* Updating - use newly calculted estimates to update

comp p11 = p11n
comp p21 = p21n
comp p12 = p12n
comp p22 = p22n
comp a1=a1n
comp a2=a2n
comp phi=phin

comp sigma=sigman

dis "Probabilities"

dis p11 p12
dis p21 p22

dis "Standard Error"
dis sigma
dis "Means"
dis a1
dis a2

dis "phi"
dis phi

*compute conv = %if(%abs(diff1)+%abs(diff2)<tol,1.0,0.0)
*compute diff=%abs(diff1)

compute conv1=conv1+1
}
end while
statistics g
*
*  Hamilton's GDP 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=1
compute nstates=2**(nlags+1)
dec vect mu(2) rhoh(nlags)
dec vect tp(5)
*
nonlin mu rhoh 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
*
*  This is just an extension to a large number of states of the simpler
*  formulas from the manual, making heavy use of the %do function. The %(...)
*  function allows complex expressions to be embedded within a function call.
*
frml ggrowth = $
   %do(i,1,nstates,reg(i)=$
  %(u=g-mu(lagstate(i,1)),%do(j,1,nlags,u=u-rhoh(j)*(g{j}-mu(lagstate(i,j+1)))),u)),$
   %do(i,1,nstates,f(i)=%density(reg(i)/sigma)/sigma),$
   %do(i,1,nstates,p(i)=%(sp=0.0,%do(j,1,nstates,sp=sp+tp(transit(i,j))*pstar(j)),sp)),$
   fpt=0.0,%do(i,1,nstates,fpt=fpt+(fp(i)=f(i)*p(i))),$
   %do(i,1,nstates,pstar(i)=fp(i)/fpt),prob1{0}=%dot(mask,pstar),log(fpt)
*
*  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.
*  With this much data, the initial distribution is unlikely to matter much, so
*  we just give them all an equal probability.
*
frml initial = tp=||0.0,p11,1-p22,1-p11,p22||,%do(i,1,nstates,pstar(i)=1.0/nstates)
*
*  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 mu(1)=a1,mu(2)=a2
linreg g
# constant g{1}
compute rhoh=%xsubvec(%beta,2,2)
maximize(start=initial,method=bfgs,trace) ggrowth nlags+2 *
graph(style=polygon,overlay=line) 2
# g * %regend()
# prob1 %regstart() %regend()