* * 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>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)>%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