*
* Example 12.6 from pp 591-597
* Caution: even though this uses a coarser grid for drawing the alpha's and takes
* fewer draws, it still takes a long time.
*
open data m-geln.txt
calendar(m) 1926
data(format=free,org=columns) 1926:1 1999:12 ge
*
* Estimate GARCH-M model. Because this uses the square root of the variance (not
* the variance itself), it can't be done with the GARCH instruction.
*
nonlin c a1 a2 b1 dh
dec series h u
frml uf = ge-dh*sqrt(h)
frml hf = c+a1*h{1}+a2*h{2}+b1*u{1}**2
*
linreg(noprint) ge
# constant
set u = %resids
set h = %seesq
compute c=%seesq/(1-.85),a1=a2=.4,b1=.05,dh=0.0
*
frml garchm = h=hf,u=uf,%logdensity(h,u)
maximize garchm 3 *
*
* Initialize e's and draw a random sample for the s series based upon the
* transition probabilities.
*
compute e1=e2=.1
set(first=%raninteger(1,2)) s = 1+(%uniform(0.0,1.0)>%if(s{1}==1,1-e1,e2))
*
* Prior for beta's
*
compute beta10=0.3,hbeta10=1.0/.09
compute beta20=1.3,hbeta20=1.0/.09
compute beta1=beta10,beta2=beta20
*
* Priors for alpha's
*
dec rect alpha(2,3)
compute alpha(1,1)=1.0,alpha(1,2)=0.6,alpha(1,3)=0.2
compute alpha(2,1)=2.0,alpha(2,2)=0.7,alpha(2,3)=0.1
*
* Formulas for the variances in the two branches
*
frml h1f = alpha(1,1)+alpha(1,2)*h{1}+alpha(1,3)*u{1}**2
frml h2f = alpha(2,1)+alpha(2,2)*h{1}+alpha(2,3)*u{1}**2
*
* Set up the grids for the alpha parameters
*
compute ngrid=100
dec rect[vect] agrid(2,3)
compute agrid(1,1)=agrid(2,1)=%seqa(0.0,6.0/ngrid,ngrid)
compute agrid(1,2)=agrid(2,2)=%seqa(0.0,1.0/ngrid,ngrid)
compute agrid(1,3)=agrid(2,3)=%seqa(0.0,0.5/ngrid,ngrid)
compute [vector] fgrid=%zeros(ngrid,1)
*
* Priors for e's
*
compute gamma11=5,gamma12=95
compute gamma21=5,gamma22=95
*
* nburn is the number of draws in the burn-in period
* nkeep is the number of draws to keep
*
compute nburn=100,nkeep=400
*
* We're keeping track of the full distribution on ten statistics
* (the beta's, alpha's and e's) and the time series for the states
*
dec vect[series] stats(10)
do i=1,10
   set stats(i) 1 nkeep = 0.0
end do i
*
set shat 1 1999:12 = 0.0
*
dec real lastlogl
*
infobox(action=define,progress,lower=-(nburn-1),upper=nkeep) "Gibbs Sampling"
do draws=-(nburn-1),nkeep
   infobox(current=draws)
   *
   * Calculate the residuals given the current parameters
   *
   set u 2 * = h=%if(s==1,h1f,h2f),ge-sqrt(h)*%if(s==1,beta1,beta2)
   *
   * Draw beta's given states
   *
   set rstd = ge/sqrt(h)
   sstats(smpl=(s==1)) 2 * rstd>>sum1
   compute vbeta1=1.0/(%nobs+hbeta10),mbeta1=vbeta1*(sum1+hbeta10*beta10),beta1=mbeta1+%ran(sqrt(vbeta1))
   sstats(smpl=(s==2)) 2 * rstd>>sum2
   compute vbeta2=1.0/(%nobs+hbeta20),mbeta2=vbeta2*(sum2+hbeta20*beta20),beta2=mbeta2+%ran(sqrt(vbeta2))
   *
   * Draw alpha's
   *
   do i=1,2
      do j=1,3
         do k=1,ngrid
            compute alpha(i,j)=agrid(i,j)(k)
            set u 2 * = h=%if(s==1,h1f,h2f),ge-sqrt(h)*%if(s==1,beta1,beta2)
            sstats(smpl=(s==i)) 2 * %logdensity(h,u)>>logl
            compute fgrid(k)=logl
            if (j==2.or.j==3).and.alpha(i,2)+alpha(i,3)>=1.0
               break
         end do k
         compute nused=k
         *
         * One detail in doing the draws from the grid is that you have to make sure that
         * you don't either overflow or underflow when converting from log densities. The
         * simplest way to do this is to subtract off the largest value. That way the
         * "density" function will range between 0 and 1.
         *
         compute biglogl=%maxvalue(%xsubvec(fgrid,1,nused))
         ewise fgrid(k)=%if(k<=nused,exp(fgrid(k)-biglogl),0.0)
         compute alpha(i,j)=%rangrid(%xsubvec(agrid(i,j),1,nused),%xsubvec(fgrid,1,nused))
      end do j
   end do i
   *
   * Draw e's
   *
   sstats(smpl=(s{1}==1)) 2 * s==2>>count
   compute e1=%ranbeta(gamma11+count,gamma12+%nobs-count)
   sstats(smpl=(s{1}==2)) 2 * s==1>>count
   compute e2=%ranbeta(gamma21+count,gamma22+%nobs-count)
   *
   * Draw s's. Since the first and last data points have only one neighbor, they are
   * covered by a different transition calculation. For each time period, we need to
   * compute the function value twice, one for S=1, one for S=2. Since we end up keeping
   * one of those two values, we've already done one of the calculations that we need
   * at the next stage. lastlogl is that value.
   *
   do time=1,1999:12
      compute sold=s(time)
      if sold==1.and.time>1
         compute logl1=lastlogl
      else {
         compute s(time)=1
         set u 2 * = h=%if(s==1,h1f,h2f),ge-sqrt(h)*%if(s==1,beta1,beta2)
         sstats 2 * %logdensity(h,u)>>logl1
      }
      if time==1
          compute p1g=%if(s(time+1)==1,1-e1,e1)
      else
      if time==1999:12
         compute p1g=%if(s(time-1)==1,1-e1,e2)
      else
         compute p1g=%if(s(time-1)==1,1-e1,e2)*%if(s(time+1)==1,1-e1,e1)

      if sold==2.and.time>1
         compute logl2=lastlogl
      else {
         compute s(time)=2
         set u 2 * = h=%if(s==1,h1f,h2f),ge-sqrt(h)*%if(s==1,beta1,beta2)
         sstats 2 * %logdensity(h,u)>>logl2
      }
      if time==1
         compute p2g=%if(s(time+1)==1,e2,1-e2)
      else
      if time==1999:12
         compute p2g=%if(s(time-1)==1,e1,1-e2)
      else
         compute p2g=%if(s(time-1)==1,e1,1-e2)*%if(s(time+1)==1,e2,1-e2)
      *
      * Again, you have to be careful to avoid over/underflows in converting to densities.
      * Since the probability depends only upon the ratios of the densities, we can
      * subtract the two log densities before taking exp's
      *
      compute s(time)=%ranbranch(||p1g*exp(logl1-logl2),p2g||)
      compute lastlogl=%if(s(time)==1,logl1,logl2)
   end do time
   if draws>0 {
      set shat 1 1999:12 = shat+s
      compute stats(1)(draws)=beta1
      compute stats(2)(draws)=beta2
      compute stats(3)(draws)=e1
      compute stats(4)(draws)=e2
      compute stats(5)(draws)=alpha(1,1)
      compute stats(6)(draws)=alpha(1,2)
      compute stats(7)(draws)=alpha(1,3)
      compute stats(8)(draws)=alpha(2,1)
      compute stats(9)(draws)=alpha(2,2)
      compute stats(10)(draws)=alpha(2,3)
    }
end do draws
infobox(action=remove)
*
* The probability of state 2 is the average of shat (the state) - 1.0
*
set pprob 1 1999:12 = (shat/nkeep)-1.0
*
spgraph(vfields=2,footer="Markov Switching GARCH-M model")
graph(header="(a) Monthly log returns")
# ge
graph(header="(b) Posterior Probability of state 2")
# pprob
spgraph(done)
*
dec vect[strings] vlabels(10)
compute vlabels=||"beta1","beta2","e1","e2","alpha10","alpha11","alpha12","alpha20","alpha21","alpha22"||
do i=1,10
   density(type=histogram,counts,maxgrid=12) stats(i) 1 nkeep xx fx
   scatter(style=bar,hlabel=vlabels(i))
   # xx fx
end do i