/* 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;