/* Econometrics Project: Markov Regime Switching Model
   Code by Jamus Jerome Lim, based on original code by James D. Hamilton
   October 25, 2003 */

output file=C:\gauss\moose\regime1.out reset;
format /m1 /ro 14,6;

/* RXREST

   This program is written in the GAUSS programming language.  It makes use of
   the procedure SMOOTH3 as well as the GAUSS procs HESSP and GRADP. It is
   modification of the EMEST program orginally written by James D. Hamilton.

   This program (RXREST) estimates a 2-state Markov switching model
   in which the mean and the variance of an nv-dimensional vector vary with
   the state.  Estimation is via the EM algorithm described in "Analysis of
   Time Series Subject to Changes in Regime," by James D. Hamilton
   (Journal of Econometrics, July 1990).

    Variable meanings:
            n = number of observations
            nv = dimension of vector (number of variables)
            kc = 1 suppresses output print
            kc = 2 prints probabilities and smoothed probabilities
            kb = 1 directs SMOOTH3 to just calculate objective function
            kb = 2 directs SMOOTH3 to calculate smoothed probabilities
            kb = 4 directs SMOOTH3 to iterate once on normal equations */

/* This section fits a scalar process for the real exchange rate, as in "Long Swings in the
    Exchange Rate: Are They in the Data and Do Markets Know It?" by Charles
    Engel and James D. Hamilton (AER, Sept. 1990).  */

   nv = 1;                     			@ nv is the number of variables analyzed @
   n = 94;                     	 		@ n is the number of observations on each variable @
   load y[n,35] = C:\gauss\moose\exrate.dat;	@ this loads exchange rate data @
   y = y[.,1];                 	 		@ this selects country 1 @ 
/* y = y[.,5];                 	 		@ this would select country 5 @ */

@  The next three paramaters (a,b, and c) determine the strength
       of the Bayesian prior imposed.  One is essentially proceeding as if
       one had c observations from each regime with mean zero
       and a observations from each regime with variance b/a.
       See the paper by James D. Hamilton, "A quasi- Bayesian approach 
       to estimating paramters for mixtures of normal distributions" ,
       Journal of Business and Econoomic Statistics, January 1991 @

   a = .2;
   b = 1.0;
   c = .1;

@  Next set initial guesses for parameters in the vector th, whose dimension
      should be 2*(nv+nv^2+1) by 1.  Parameters are listed in the following order:

          Mean (or mean vector) for state 1
          Mean (or mean vector) for state 2
          The Markov probability for going from state 1 to state 1
          The Markov probability for going from state 2 to state 2
          The (nv by nv) variance-covariance matrix for state 1
          The (nv by nv) variance-covariance matrix for state 2  @

   let th[6,1] =  5 -5 .5 .5 20 20;

   nvmu2 = nv + 1;             		@ nvmu2 is the index at which mu2 parameters begin @
   nvp = (2*nv) + 1;           		@ nvp is the index for p @
   nvq  = (2*nv)+2;            		@ nvq is the index for q @
   nvsig1 = (2*nv) + 3;        		@ nvsig1 is the index at which sig1 parameters begin @
   nvsig2 = (2*nv) + 3 + (nv^2); 	@ nvsig2 is the index at which sig2 parameters begin @
   nvend = 2*(nv+nv^2+1);      		@ nvend is the index for the last element @

   "Parameters for Bayesian prior (a, b, c)";
    a;b;c;"";

   kb = 4;
   kc = 1;
   it = 1;
   ierr = 0;             	@ a value of 1 indicates an 80287 exception encountered @
   dh = 1.e-8;         		@ definition of small step for gradient and hessian @
   deltalik = 10;      		@ change in likelihood function @
   deltach = zeros(6,1);	@ change in parameter vector over previous iteration @
   deltath  = 10;      		@ biggest element of deltach @
   thtol  = 1.e-8;     		@ when change in parameters is less than this, convergence is obtained @
   mxit = 400;         		@ maximum number of iterations allowed @
   vof = 0;

#include smooth3.gss;

     "";"";
          "starting parameter values";
          "   Means for state 1";
              th[1:(nvmu2-1),1]';
           "   Means for state 2";
               th[nvmu2:(nvp-1),1]';
           "   Prob from state 1 to state 1 (p)";
               th[nvp,1];
            "   Prob from state 2 to state 2 (q)";
                th[nvq,1];
            "   Variance-covariance matrix for state 1";;
                reshape(th[nvsig1:(nvsig2-1),1],nv,nv);
            "   Variance-covariance matrix for state 2";;
                reshape(th[nvsig2:nvend,1],nv,nv);"";

      kb =1; vof = ofn(th); kb = 4;

@ iterate on normal equations @

     it = 1;
     do until (it > mxit) or (deltath < thtol) or (ierr == 1);
           th = ofn(th);
     it = it+1;
     endo;

     "";"      Number of iterations required:";;it;
        "      Change in theta vector";;deltath;


@ print out final results @

     kb=1;"";
     "";"             FINAL RESULTS";"";

"Estimated theta vector:";th';"";

"Estimated parameters:";
   "   Means for state 1";
        th[1:(nvmu2-1),1]';
   "   Means for state 2";
        th[nvmu2:(nvp-1),1]';
   "   Prob from state 1 to state 1 (p)";
        th[nvp,1];
   "   Prob from state 2 to state 2 (q)";
        th[nvq,1];
   "   Variance-covariance matrix for state 1";;
        reshape(th[nvsig1:(nvsig2-1),1],nv,nv);
   "   Variance-covariance matrix for state 2";;
        reshape(th[nvsig2:nvend,1],nv,nv);"";
     "Value of log likelihood function:";;vof=ofn(th);vof;


@ To calculate standard errors and gradient, we need to reconfigure
  parameter vector to eliminate redundant terms when nv > 1 @

xk = zeros(nvq + (nv+1)*nv,1);   @ xk will include nonredundant terms @

If nv > 1;
   xk[1:nvq,1] = th[1:nvq,1];
   k1 = nvsig1;
   k2 = nvsig1;

   i = 1;
   do until i > nv;
      j = 1;
      do until j > nv;
         if j >= i;
            xk[k1,1] = th[k2,1];
            k1 = k1 + 1;
            k2 = k2+1;
         else;
            k2 = k2 + 1;
         endif;
      j = j+1;
      endo;
   i = i+1;
   endo;

   i = 1;
   do until i > nv;
      j = 1;
      do until j > nv;
         if j >= i;
            xk[k1,1] = th[k2,1];
            k1 = k1 + 1;
            k2 = k2+1;
         else;
            k2 = k2 + 1;
         endif;
      j = j+1;
      endo;
   i = i+1;
   endo;

   "Parameter vector has been reconfigured as follows:";
   xk';

endif;

@ ----------------------------------------------------------------  @

proc convert(xk);  			@ This proc converts to redundant parameterization @
   local qth,i,j,k1,k2,qsig;
   qth = zeros(nvend,1);
   qsig = zeros(nv,nv);
   qth[1:nvq,1] = xk[1:nvq,1];
   k1 = nvsig1;
   k2 = nvsig1;

   i = 1;
   do until i > nv;
       j = 1;
       do until j > nv;
          if j >= i;
              qsig[i,j] = xk[k1,1];
              k1 = k1 + 1;
          else;
              qsig[i,j] = qsig[j,i];
          endif;
      j = j+1;
      endo;
   i = i+1;
   endo;
   qth[nvsig1:(nvsig2 - 1),1] =
         reshape(qsig,(nv*nv),1);

   i = 1;
   do until i > nv;
       j = 1;
       do until j > nv;
          if j >= i;
              qsig[i,j] = xk[k1,1];
              k1 = k1 + 1;
          else;
              qsig[i,j] = qsig[j,i];
          endif;
      j = j+1;
      endo;
   i = i+1;
   endo;
   qth[nvsig2:nvend,1] =
         reshape(qsig,(nv*nv),1);

 retp(ofn(qth));
 endp;

@ ------------------------------------------------------------------------ @

"";

    "Variance-covariance matrix for parameters";
     if nv == 1;
         h = hessp(&ofn,th);
     else;
         h = hessp(&convert,xk);
     endif;
     bbb = eigrs(-h);
         if bbb[1,1] <= 0;
             "Hessian not positive definite";
         else;
              vc = invpd(-h);
              format /m3 /ro 14,6;
              vc;
              "";"Standard errors";;
              format /m0;
              stderr = sqrt(diag(vc));stderr;
          endif;
     "";kc = 2; call ofn(th);

@ Confirm that maximum has been found by calculating gradient @

    kc = 1;
    "";"numerical gradient:";
    if nv == 1;
       gradp(&ofn,th);
    else;
       gradp(&convert,xk);
    endif;

output off;