*
* MIXTURE.RPF
* Example of a mixture model (not-Markovian) demonstrating the different
* estimation strategies.
*
* RATS User's Guide, Section 11.7.1.
*
open data ger4ind.prn
calendar(q) 1961
data(format=prn,org=columns) 1961:1 1986:4 ipgrowth
set x = 100.0*ipgrowth
*
* We'll use a two lag AR which will fully switch between regimes. In
* practice, this is probably too "loose" a model; there are almost
* certainly many modes for the likelihood function, since AR(2)'s can
* model many types of behavior.
*
equation xeq x
# constant x{1 to 2}
*
dec vect[vect] phi(2)
nonlin(parmset=regparms) phi sigsq
nonlin(parmset=msparms)  p
*
* Set up guess values. The coefficient vectors are separated by giving
* the first a larger intercept.
*
linreg(equation=xeq)
compute rstart=%regstart(),rend=%regend()
compute phi(1)=%beta
compute phi(2)=%beta
compute phi(1)(1)+=sqrt(%seesq)
compute phi(2)(1)-=sqrt(%seesq)
compute sigsq=%seesq
*
* Save the guess values so we can use them when demonstrating different
* techniques. (This step isn't generally necessary).
*
compute guesses=%parmspeek(regparms)
*
* This is the function which returns the regime likelihoods.
*
function RegimeF time
type vector  RegimeF
type integer time
local integer i
*
dim RegimeF(2)
ewise RegimeF(i)=exp(%logdensity(sigsq,%eqnrvalue(xeq,time,phi(i))))
end
*
* Maximum likelihood
*
frml logl = f=RegimeF(t),log(p*f(1)+(1-p)*f(2))
*
compute p=.3
maximize(parmset=regparms+msparms,pmethod=simplex,piters=5) logl rstart rend
*
* Compute the probability of the regimes at each time period.
*
set pt_t = f=RegimeF(t),p*f(1)/(p*f(1)+(1-p)*f(2))
graph(footer="Filtered Probability and Data",overlay=line) 2
# pt_t
# x
*
* EM
*
frml emlogl = f=RegimeF(t),pt_t*log(p*f(1))+(1-pt_t)*log((1-p)*f(2))
*
function EStep
set pt_t = f=RegimeF(t),p*f(1)/(p*f(1)+(1-p)*f(2))
end
*
* Reset the guess values.
*
compute %parmspoke(regparms,guesses)
compute p=.3
*
* Do 50 EM iterations, then switch over to ML. (EM itself doesn't
* converge in even 500 iterations). Note that the EM log likelihood
* isn't comparable to the ML log likelihood.
*
maximize(estep=Estep(),method=bhhh,parmset=regparms+msparms,iters=50) emlogl rstart rend
maximize(method=bfgs,parmset=regparms+msparms,hessian=%xx) logl rstart rend
*
* Bayesian MCMC
* Note that this has a severe "labeling" problem---the model doesn't
* really want to separate the regimes in the way that we intend.
*
compute %parmspoke(regparms,guesses)
compute p=.3
*
compute r1=r2=5
compute pdf=1,pscale=0.0
compute nburn=1000,ndraws=1000
dec series[integer] s
dec vect[vect] bphi(2)
dec vect[symm] hphi(2)
compute bphi(1)=phi(1)
compute bphi(2)=phi(2)
compute hphi(1)=%diag(||0.0,0.0,0.0||)
compute hphi(2)=%diag(||0.0,0.0,0.0||)
*
* Smallest regime size
*
compute minsize=5
*
set pstar rstart rend = 0.0
infobox(action=define,progress,lower=-nburn,upper=ndraws) "Gibbs sampling"
do draw=-nburn,ndraws
   *
   * Draw regimes given regression parameters and p. If the size of either regime is too small, reject and redraw.
   *
:redrawregimes
   gset s rstart rend = f=RegimeF(t),%ranbranch(||p*f(1),(1-p)*f(2)||)
   sstats rstart rend (s==1)>>count1 (s==2)>>count2
   if count1<minsize.or.count2<minsize {
      disp "Draw" draw "Redrawing regimes with regime of size" $
         %min(count1,count2)
      goto redrawregimes
   }
   *
   * Draw p given s
   *
   sstats rstart rend s==1>>count1
   compute p=%ranbeta(count1+r1,%nobs-count1+r2)
   *
   * Draw phi's given s, sigsq
   *
   cmom(smpl=(s==1),equation=xeq) rstart rend
   compute phi(1)=%ranmvpostcmom(%cmom,1.0/sigsq,hphi(1),bphi(1))
   cmom(smpl=(s==2),equation=xeq) rstart rend
   compute phi(2)=%ranmvpostcmom(%cmom,1.0/sigsq,hphi(2),bphi(2))
   *
   * Relabel if necessary
   *
   if phi(1)(1)<phi(2)(1) {
      disp "Relabel"
      compute temp=phi(1)
      compute phi(1)=phi(2)
      compute phi(2)=phi(1)
      gset s = 3-s
   }
   *
   * Draw sigsq given phi's and s
   *
   sstats rstart rend %eqnrvalue(xeq,t,phi(s))^2>>%rss
   compute sigsq=(pscale*pdf+%rss)/%ranchisqr(%nobs+pdf)
   infobox(current=draw)
   if draw>0
      set pstar rstart rend = pstar+(s==1)
end do draw
infobox(action=remove)
set pstar rstart rend = pstar/ndraws
graph(footer="Smoothed Probability of the Regime 1")
# pstar