### Aptech Systems, Inc. Worldwide Headquarters

Address:

Aptech Systems, Inc.

2350 East Germann Road, Suite #21

Chandler, AZ 85286Phone: 360.886.7100

FAX: 360.886.8922**Ready to Get Started?**### Request Quote & Product Information

### Industry Solutions

### Products

### Resources

### Support

### Training & Events

Want more guidance while learning about the full functionality of GAUSS and its capabilities? Get in touch for in-person training or browse additional references below.

### Tutorials

Step-by-step, informative lessons for those who want to dive into GAUSS and achieve their goals, fast.

### What’s New in GAUSS 17

### Want to find out more?

### Have a Specific Question?

### Q&A: Register and Login

### Support Plans

Premier Support and Platinum Premier Support are annually renewable membership programs that provide you with important benefits including technical support, product maintenance, and substantial cost-saving features for your GAUSS System or the GAUSS Engine.

### User Forums

Join our community to see why our users are considered some of the most active and helpful in the industry!

### Where to Buy

Available across the globe, you can have access to GAUSS no matter where you are.

### Recent Tags

applications character vectors CML CMLMT Constrained Optimization datasets dates dlibrary dllcall econometrics Editor error error codes error handling errors Excel file i/o floating network GAUSS Engine GAUSS Light graphics GUI hotkeys installation license licensing linux loading data loop loops matrix manipulation Maximum Likelihood Maxlik MaxLikMT Memory multidimensional array optimization Optmum output PQG graphics procs random numbers strings structures threading### Recent Questions

- Matrices are not conformable (ss=gf*omega*gf')
- G0121 matrix not positive definite with my data
- G0025 : Undefined symbol: 't'
- some question about installing library in student lite version in Mac OS
- Why not load the library in student lite version?
- Group Column Sum
- A C program
- Mixed-effects models
- envelope
- "do while" command problems

### Features

### Time Series 2.0 MT

### Industry Solutions

### Find out more now

### Time Series MT 2.1

### Find out more now

### Find out more now

# Resources

# Error G0008

While running the code "MAXSEEK" (downloaded from Hamilton's website), I get the following error. Hamilton said he did not face the following errors. Can someone please advise? Is there any option of attaching a file?

Line 153 in C:\gauss10\mrs\procs

Syntax error G0008 : 'endp'

Line 182 in C:\gauss10\mrs\MAXSEEK

Nested procedure definition G0155

Line 182 in C:\gauss10\mrs\MAXSEEK

Syntax error G0008 : 'proc matpm(xth)'

## 2 Answers

**MAXSEEK and PROCS codes are given below **

MAXSEEK

@ This GAUSS file reads in data, sets options, and calls numerical

optimization routine for numerical estimation of switching model @

output file=junk reset;

/* ======================================================================= */

@ Replace this first section with a section to read in your data

Example 1 reproduces the GNP results from Hamilton, "A New

Approach to the Economic Analysis of Nonstationary Time Series

and the Business Cycle," Econometrica, March 1989

Example 2 reproduces the interest rate results from Hamilton,

"Rational Expectations Econometric Analysis of Changes in Regime:

An Investigation of the Term Structure of Interest Rates," Journal

of Economic Dynamics and Control, June 1988

Example 3 reproduces the results for the French exchange rate 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).

Example 4 fits a three-state model to US real GNP data; estimation

algorithm bogs down due to bumping against corners with nonnegative

definite hessian; Example 5 reestimates with zeros imposed @

/*===============================================

EXAMPLE 1:*/

@ The following lines read in the raw data, and convert to

100 times the first difference of the log @

capt = 135; @ capt is the sample size @

load bbbb[capt+1,1] = C:\gauss10\swarch\gnp82.dat;

bbbb = bbbb[1:capt,1]~bbbb[2:capt+1,1];

y = 100*ln(bbbb[.,2]./bbbb[.,1]);

@ Adjust any of the following to control specification desired @

ns = 2; @ ns is the number of primitive states @

ps = 4; @ ps is the number of lagged states that matter

for y; use ps = 0 if only the current state

matters @

pphi = 4; @ pphi is the number of lags in autoregression

for y; pphi should be greater than or equal

to ps @

isig = 1; @ isig = 1 for constant variances,

isig =ns for changing variances @

ipm = 1; @ ipm specifies way in which transition probs

are parameterized

ipm = 1 implies p11 and p22 estimated

ipm = 2 implies pij for i=1,..,n j=1,..,n-1

ipm = 3 user input code @

a=0;b=0;c=0; @ These parameters control the Bayesian prior

as suggested in Hamilton, "A Quasi-Bayesian

Approach to Estimating Parameters for Mixtures

of Normal Distributions", Journal of Business

and Economic Statistics, January 1991. If

no Bayesian adjustment is desired, set all

parameters to zero @

/* Input initial values for parameters

The order in which variables are represented is as follows

first ns elements: means for states

next pphi elements: autoregressive coefficients

next isig elements: variances for states when izz =1

std. deviations for states when izz =2

next elements: when izz =1 these are the transition probs

when izz =2 these are params v(i,j) such that

p(i,j) = v(i,j)^2 / sum j v(i,j)^2

elements are ordered as p11, p22 when ns =2

ordered as p(1,1),p(2,1),...,p(ns,1),

p(2,1),...,p(ns-1,ns) when ns > 2 */

nth = 9; @ nth is the number of params to be estimated@

mu1 = 1; mu2 = 0; phi1 = .1; phi2 = 0; phi3 = 0; phi4 = 0;

sig = 1; p11 = 1.5; p22 = 1.5;

th = mu1|mu2|phi1|phi2|phi3|phi4|sig|p11|p22;

/* =============================================

EXAMPLE 2:

capt = 103;

load y[capt,1] = rradj.dat;

ns = 2;

ps = 4;

pphi = 4;

isig = 2;

ipm = 1;

a=0;b=0;c=0;

nth = 10;

mu1 = 2; mu2 = 1; phi1 = .9; phi2 = 0.; phi3 = 0.;

phi4 = 0; sig1 = .8; sig2 = .2; p11 = 1; p22 = 1;

th = mu1|mu2|phi1|phi2|phi3|phi4|sig1|sig2|p11|p22; */

/*==============================================

EXAMPLE 3

capt = 58;

load y[capt,3] = exdata.txt; @ this loads exchange rate data

for France, U.K. and Germany @

y = y[.,1]; @ this selects French data @

ns = 2;

ps = 0;

pphi = 0;

isig = 2;

ipm = 1;

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

nth = 6;

mu1 = 3; mu2 = -3; sig1 = 3; sig2 = 6; p11 = 2; p22 = 2;

th = mu1|mu2|sig1|sig2|p11|p22; */

/* ============================================================

EXAMPLE 4

capt = 135; @ capt is the sample size @

load bbbb[capt+1,1] = gnp82.dat;

bbbb = bbbb[1:capt,1]~bbbb[2:capt+1,1];

y = 100*ln(bbbb[.,2]./bbbb[.,1]);

a=.1;b=.1;c=.1;

ns = 3; ps = 1; pphi = 1; isig = 1; ipm = 2;

nth = 11;

mu1 = -1; mu2 = .2; mu3 = 1.2; phi1 = .1; sig = .8;

p11 = 1.1; p21 = 1.8; p31 = .1; p12 = .15; p22 = 3; p23 = .33;

th = mu1|mu2|mu3|phi1|sig|p11|p21|p31|p12|p22|p23; */

/* ============================================================

EXAMPLE 5

capt = 135; @ capt is the sample size @

load bbbb[capt+1,1] = gnp82.dat;

bbbb = bbbb[1:capt,1]~bbbb[2:capt+1,1];

y = 100*ln(bbbb[.,2]./bbbb[.,1]);

a=.1;b=.1;c=.1;

ns = 3; ps = 1; pphi = 1; isig = 1; ipm = 3;

nth = 8;

mu1 = -1; mu2 = .2; mu3 = 1.2; phi1 = .1; sig = .8;

p11 = 1; p21 = 1; p23 = 1;

th = mu1|mu2|mu3|phi1|sig|p11|p21|p23; */

/* ======================================================================= */

/* ======================================================================= */

@ In general no parts of this section should be changed @

proc startval; @ This defines starting value for iteration to be th @

retp(th); endp;

nk = pphi+1; @ nk is the first observation for which the

likelihood will be evaluated @

izz = 1; @ izz = 1 when params read in so as to assure inequality

constraints; izz = 2 for final reporting of results @

n = ns^(ps+1); @ n is the dimension of the state vector @

kc = 1; @ kc = 2 to echo parameter values @

ks = 1; @ ks = 2 if smoothed probs are to be calculated @

captst = capt - nk +1; @ captst is the effective sample size @

skif = zeros(captst,n); @ skif is the matrix of filtered probs @

skis = zeros(captst,n); @ skis is the matrix of smoothed probs @

id = eye(ns); @ used in certain calculations below @

"Bayesian prior used";

"a=";;a;;"b=";;b;;"c=";;c;

proc pattern1; @ This proc returns a (ps+1)*ns x n matrix. The ith

column contains a one in row j if st = j, contains a

one in row ns+j if st-1 = j, and so on @

local i1,ix,iq,na;

na = n/ns;

ix = eye(ns).*.ones(1,na);

i1 = 1;

do until i1 > ps;

na = na/ns;

iq = ones(1,ns^i1).*.(eye(ns).*.ones(1,na));

ix = iq|ix;

i1 = i1+1;

endo;

retp(ix);

endp;

hp = pattern1;

#include procs;

/* ================================================================= */

/* ======================================================================= */

proc matpm(xth); @This proc defines the user's conventions for reading

elements of Markov transition probabilities from

parameter vector @

local pm,ixth;

ixth = rows(xth);

pm = zeros(ns,ns);

if ipm == 1; @ for ns =2 this option has parameters as p11 and p22 @

if izz == 1;

pm[1,1] = xth[1,1]^2/(1 +xth[1,1]^2);

pm[2,2] = xth[2,1]^2/(1 + xth[2,1]^2);

else;

pm[1,1] = xth[1,1];

pm[2,2] = xth[2,1];

endif;

pm[2,1] = 1 - pm[1,1];

pm[1,2] = 1 - pm[2,2];

elseif ipm == 2; @ general case has parameters pij for i = 1,...,n and

j = 1,...,n-1 @

pm[1:ns-1,.] = reshape(xth[1:ixth,1],ns-1,ns);

if izz == 1;

pm[ns,.] = ones(1,ns);

pm = pm^2;

pm = pm./(sumc(pm)');

else;

pm[ns,.] = (1 - sumc(pm))';

endif;

elseif ipm == 3; @ This section can be rewritten by user to impose zeros

and ones where desired @

if izz == 1;

pm[1,1] = xth[1,1]^2/(1 + xth[1,1]^2);

pm[1,2] = xth[2,1]^2/(1 + xth[2,1]^2);

pm[2,3] = xth[3,1]^2/(1 + xth[3,1]^2);

elseif izz == 2;

pm[1,1] = xth[1,1];

pm[1,2] = xth[2,1];

pm[2,3] = xth[3,1];

endif;

pm[3,1] = 1 - pm[1,1];

pm[2,2] = 1 - pm[1,2];

pm[3,3] = 1 - pm[2,3];

endif;

retp(pm);

endp;

/* ======================================================================= */

/* ==================================================================== */

@ Set parameters to use Gauss numerical optimizer @

library optmum;

#include optmum.ext;

__btol = 1.e-06; @ This controls convergence criterion for coefficients@

__gtol = 1.e-06; @ This controls convergence criterion for gradient @

__algr = 1; @ This chooses BFGS optimization @

_opalgr = 2; @ This chooses BFGS optimization for newer versions of

GAUSS @

__miter = 150; @ This controls the maximum number of iterations @

__output = 1; @ This causes extra output to be displayed @

__covp = 0; @ This speeds up return from OPTMUM; note that the

program makes a reparameterization to calculate

std. errors @

@ Next call the GAUSS numerical optimizer @

output off;

{x,f,g,h} =optmum(&ofn,startval);

output file=junk on;

"";"";"MLE as parameterized for numerical optimization ";

"Coefficients:";x';

"";"Value of log likelihood:";;-f;

"";"Gradient vector:";g';

/* ======================================================================= */

/*======================================================================== */

@ In general no parts of this section need be changed @

@ Reparameterize for reporting final results @

izz = 2;

"";"Vector is reparameterized to report final results as follows";

"Means for each state:";x[1:ns,1]';

If pphi > 0;

"Autoregressive coefficients:";x[ns+1:ns+pphi,1]';

endif;

ncount = ns + pphi ;

x[ncount+1:ncount+isig,1] = x[ncount+1:ncount+isig,1]^2;

"Variances:";x[ncount+1:ncount+isig,1];

ncount = ncount + isig;

/* ======================================================================= */

/* ======================================================================= */

proc pmth; @ This proc converts the last elements of parameter vector

(which relates to transition probabilities) from the

th[i,j]^2/{sum j th[i,j]^2} form that is used for

numerical estimation into the p[i,j] form that is used

to calculate standard errors @

local pm;

if ipm == 1;

x[ncount+1,1] = x[ncount+1,1]^2/(1 + x[ncount+1,1]^2);

x[ncount+2,1] = x[ncount+2,1]^2/(1 + x[ncount+2,1]^2);

elseif ipm == 2;

pm = zeros(ns,ns);

pm[1:ns-1,.] = reshape(x[ncount+1:nth],ns-1,ns);

pm[ns,.] = ones(1,ns);

pm = pm^2;

pm = pm./(sumc(pm)');

x[ncount+1:nth,1] = reshape(pm[1:ns-1,.],ns*(ns-1),1);

elseif ipm == 3; @ User may want to alter these next lines @

x[ncount+1:nth,1]

= (x[ncount+1:nth,1]^2)./(1 + x[ncount+1:nth,1]^2);

endif;

retp(x);

endp;

/* ======================================================================= */

/* ======================================================================= */

@ In general no changes are necessary from here out @

call pmth;

h = (hessp(&ofn,x));

va = eigrs(h);

kc = 2;

ks = 2;

call ofn(x);

if minc(eigrs(h)) <= 0;

"Negative of Hessian is not positive definite";

"Either you have not found local maximum, or else estimates are up "

"against boundary condition. In latter case, impose the restricted "

"params rather than estimate them to calculate standard errors";

else;

h = invpd(h);

std = diag(h)^.5;

"For vector of coefficients parameterized as follows,";x';

"the standard errors are";std';

endif;

"";"-------------------------------";"";

"Probabilities for primitive states";

"filtered probabilities";format /rd 1,0;

"Obs ";;

t = 0;

do until t > ps;

i = 1;

do until i == ns;

"P(st-";;t;;"=";;i;;") ";;

i = i+1;

endo;

t = t+1;

endo;"";

format /rd 6,4;

skif = (skif*hp')*(eye(ps+1).*.id[.,1:ns-1]);

skif = seqa(nk,1,captst)~skif;skif;

"";"smoothed probabilities";

format /rd 1,0;

"Obs ";;

i = 1;

do until i > ns;

"P(st = ";;i;;") ";;

i = i+1;

endo;

format /rd 6,4;

skis = skis*hp';

skis = seqa(nk,1,captst)~skis[.,1:ns];skis;

/*========================================================================*/

PROCS

@ These GAUSS procedures evaluate Markov transition matrix, evaluate

likelihood function, and evaluate filter and smoothed probababilities @

/* ============================================================== */

proc matf(pm); @This proc returns the (n x n) matrix F of Markov

transition probabilities for state vector @

local iz,iw,ib,na,nb,nc,fz,fm;

@ set initial values for use with iteration @

na = 1;

nb = ns;

nc = ns*ns;

fm = pm;

iz = 1;

do until iz > ps;

fz = fm;

fm = zeros(nc,nc);

iw = 1;

do until iw > ns;

fm[((iw-1)*nb+1):(iw*nb),((iw-1)*na+1):(iw*na)]

= fz[1:nb,((iw-1)*na+1):iw*na];

iw = iw+1;

endo;

ib = 2;

do until ib > ns;

fm[1:nc,((ib-1)*nb+1):ib*nb] = fm[1:nc,1:nb];

ib = ib+1;

endo;

na = na*ns;

nb = nb*ns;

nc = nc*ns;

iz = iz+1;

endo;

retp(fm);

endp;

/* ==================================================================== */

/* ================================================================= */

proc ofn(th); @ this proc evaluates filter probs and likelihood @

local ncount,mu,phi,sig,pm,eta,iz,fm,chsi,it,f,fit,fx,hw,fj,

ij,ap,const,pq,hk,ihk;

@ Convert parameter vector to convenient form @

mu = th[1:ns,1];

ncount = ns + 1; @ ncount is the number of params read @

if pphi == 0;

phi = 1;

else;

phi = 1 | -th[ncount:ncount+pphi-1];

endif;

ncount = ncount + pphi;

if isig == 1;

sig = th[ncount,1]*ones(n,1);

ncount = ncount +1;

elseif isig == ns;

sig = th[ncount:ncount+ns-1,1].*hp[1:ns,.];

sig = sumc(sig);

ncount = ncount + ns;

endif;

if izz == 1;

sig = sig^2;

endif;

pm = matpm(th[ncount:nth,1]);

@ Construct constant term to be subtracted for each observation @

const = phi .*. mu;

const = const'*hp;

@ Convert data to AR resids @

eta = y[nk:capt,1];

iz = 1;

do until iz > pphi;

eta = eta ~y[nk-iz:capt-iz,1];

iz = iz+1;

endo;

eta = ((eta*phi - const)^2)./sig';

@ Adjust data to avoid numeric overflows and underflows

eta = reshape(eta,1,n*(captst));

eta = minc(eta|800*ones(1,n*captst));

eta = reshape(eta',captst,n); @

eta = (1./sqrt(sig')).*exp(-eta/2);

@ Calculate ergodic probabilities @

fm = matf(pm);

ap = (eye(n)-fm)|ones(1,n);

chsi = sumc((invpd(ap'*ap))');

chsi = maxc(chsi'|zeros(1,n)); @ This line eliminates roundoff error @

if kc > 1;

"";"Matrix of Markov transition probabilities:";pm;

"";"Ergodic probs for full state vector:";chsi';

"";"Ergodic probs for primitive states:";

pq = hp*chsi;pq[1:ns,1]';

endif;

@ Filter iteration @

f = 0;

it = 1;

do until it > captst;

fx = chsi .* eta[it,.]';

fit = sumc(fx);

skif[it,.] = fx'/fit;

f = f + ln(fit);

chsi = fm*fx/fit;

it = it+1;

endo;

@ Make Bayesian adjustment to likelihood if desired @

if a > 0;

fj = 0;

ij = 1;

do until ij > ns;

fj = fj + a*ln(sig[ij,1]) + b/sig[ij,1]

+ c*(mu[ij,1]^2/sig[ij,1]);

ij = ij +1;

endo;

f = f - fj/2;

endif;

@ Calculate smoothed probs if desired @

if ks == 2;

skis[captst,.] = skif[captst,.];

it = 1;

do until it == captst;

if minc(skif[captst-it,.]') > 1.e-150;

skis[captst-it,.] = skif[captst-it,.].*

((skis[captst-it+1,.]./(skif[captst-it,.]*fm'))*fm);

else; @ adjust code so as not to divide by zero @

hk = skif[captst-it,.]*fm';

ihk = 1;

do until ihk > n;

if hk[1,ihk] > 1.e-150;

hk[1,ihk] = skis[captst-it+1,ihk]/hk[1,ihk];

else;

hk[1,ihk] = 0;

endif;

ihk = ihk + 1;

endo;

skis[captst-it,.] = skif[captst-it,.].*(hk*fm);

endif;

it = it+1;

endo;

endif;

@ Print out value of log likelihood if desired @

if kc == 2;

"";"Log likelihood:";f;

endif;

retp(-f);

endp;

/*=========================================================================*/

I copied and pasted your code from above into two separate files, `markov_mod.gss` and `procs`. I had to do some find and replace, because the forum website changed some of the quotation marks to two single quotes, subtraction signs to dashes, etc. I also had to close the comment at the very last line of your post. As a note for others, I also did not copy the line that says MAXSEEK or the line that says PROCS into my files. I used those as demarcations between the two files.

After this, I downloaded the `gnp82.dat` dataset from the link that says "Programs for estimation of Markov switching models by numerical optimization" on this page. I placed the `gnp82.dat` file in the same directory as the other two files.

I was then able to run the code successfully with GAUSS version 13.1 on Linux. I have attached a zip file with these contents.

Final note: this code requires the GAUSS Optmum application module.

## Your Answer

## 2 Answers

**MAXSEEK and PROCS codes are given below **

MAXSEEK

@ This GAUSS file reads in data, sets options, and calls numerical

optimization routine for numerical estimation of switching model @

output file=junk reset;

/* ======================================================================= */

@ Replace this first section with a section to read in your data

Example 1 reproduces the GNP results from Hamilton, "A New

Approach to the Economic Analysis of Nonstationary Time Series

and the Business Cycle," Econometrica, March 1989

Example 2 reproduces the interest rate results from Hamilton,

"Rational Expectations Econometric Analysis of Changes in Regime:

An Investigation of the Term Structure of Interest Rates," Journal

of Economic Dynamics and Control, June 1988

Example 3 reproduces the results for the French exchange rate 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).

Example 4 fits a three-state model to US real GNP data; estimation

algorithm bogs down due to bumping against corners with nonnegative

definite hessian; Example 5 reestimates with zeros imposed @

/*===============================================

EXAMPLE 1:*/

@ The following lines read in the raw data, and convert to

100 times the first difference of the log @

capt = 135; @ capt is the sample size @

load bbbb[capt+1,1] = C:\gauss10\swarch\gnp82.dat;

bbbb = bbbb[1:capt,1]~bbbb[2:capt+1,1];

y = 100*ln(bbbb[.,2]./bbbb[.,1]);

@ Adjust any of the following to control specification desired @

ns = 2; @ ns is the number of primitive states @

ps = 4; @ ps is the number of lagged states that matter

for y; use ps = 0 if only the current state

matters @

pphi = 4; @ pphi is the number of lags in autoregression

for y; pphi should be greater than or equal

to ps @

isig = 1; @ isig = 1 for constant variances,

isig =ns for changing variances @

ipm = 1; @ ipm specifies way in which transition probs

are parameterized

ipm = 1 implies p11 and p22 estimated

ipm = 2 implies pij for i=1,..,n j=1,..,n-1

ipm = 3 user input code @

a=0;b=0;c=0; @ These parameters control the Bayesian prior

as suggested in Hamilton, "A Quasi-Bayesian

Approach to Estimating Parameters for Mixtures

of Normal Distributions", Journal of Business

and Economic Statistics, January 1991. If

no Bayesian adjustment is desired, set all

parameters to zero @

/* Input initial values for parameters

The order in which variables are represented is as follows

first ns elements: means for states

next pphi elements: autoregressive coefficients

next isig elements: variances for states when izz =1

std. deviations for states when izz =2

next elements: when izz =1 these are the transition probs

when izz =2 these are params v(i,j) such that

p(i,j) = v(i,j)^2 / sum j v(i,j)^2

elements are ordered as p11, p22 when ns =2

ordered as p(1,1),p(2,1),...,p(ns,1),

p(2,1),...,p(ns-1,ns) when ns > 2 */

nth = 9; @ nth is the number of params to be estimated@

mu1 = 1; mu2 = 0; phi1 = .1; phi2 = 0; phi3 = 0; phi4 = 0;

sig = 1; p11 = 1.5; p22 = 1.5;

th = mu1|mu2|phi1|phi2|phi3|phi4|sig|p11|p22;

/* =============================================

EXAMPLE 2:

capt = 103;

load y[capt,1] = rradj.dat;

ns = 2;

ps = 4;

pphi = 4;

isig = 2;

ipm = 1;

a=0;b=0;c=0;

nth = 10;

mu1 = 2; mu2 = 1; phi1 = .9; phi2 = 0.; phi3 = 0.;

phi4 = 0; sig1 = .8; sig2 = .2; p11 = 1; p22 = 1;

th = mu1|mu2|phi1|phi2|phi3|phi4|sig1|sig2|p11|p22; */

/*==============================================

EXAMPLE 3

capt = 58;

load y[capt,3] = exdata.txt; @ this loads exchange rate data

for France, U.K. and Germany @

y = y[.,1]; @ this selects French data @

ns = 2;

ps = 0;

pphi = 0;

isig = 2;

ipm = 1;

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

nth = 6;

mu1 = 3; mu2 = -3; sig1 = 3; sig2 = 6; p11 = 2; p22 = 2;

th = mu1|mu2|sig1|sig2|p11|p22; */

/* ============================================================

EXAMPLE 4

capt = 135; @ capt is the sample size @

load bbbb[capt+1,1] = gnp82.dat;

bbbb = bbbb[1:capt,1]~bbbb[2:capt+1,1];

y = 100*ln(bbbb[.,2]./bbbb[.,1]);

a=.1;b=.1;c=.1;

ns = 3; ps = 1; pphi = 1; isig = 1; ipm = 2;

nth = 11;

mu1 = -1; mu2 = .2; mu3 = 1.2; phi1 = .1; sig = .8;

p11 = 1.1; p21 = 1.8; p31 = .1; p12 = .15; p22 = 3; p23 = .33;

th = mu1|mu2|mu3|phi1|sig|p11|p21|p31|p12|p22|p23; */

/* ============================================================

EXAMPLE 5

load bbbb[capt+1,1] = gnp82.dat;

bbbb = bbbb[1:capt,1]~bbbb[2:capt+1,1];

y = 100*ln(bbbb[.,2]./bbbb[.,1]);

a=.1;b=.1;c=.1;

ns = 3; ps = 1; pphi = 1; isig = 1; ipm = 3;

nth = 8;

mu1 = -1; mu2 = .2; mu3 = 1.2; phi1 = .1; sig = .8;

p11 = 1; p21 = 1; p23 = 1;

th = mu1|mu2|mu3|phi1|sig|p11|p21|p23; */

/* ======================================================================= */

/* ======================================================================= */

@ In general no parts of this section should be changed @

proc startval; @ This defines starting value for iteration to be th @

retp(th); endp;

nk = pphi+1; @ nk is the first observation for which the

likelihood will be evaluated @

izz = 1; @ izz = 1 when params read in so as to assure inequality

constraints; izz = 2 for final reporting of results @

n = ns^(ps+1); @ n is the dimension of the state vector @

kc = 1; @ kc = 2 to echo parameter values @

ks = 1; @ ks = 2 if smoothed probs are to be calculated @

captst = capt - nk +1; @ captst is the effective sample size @

skif = zeros(captst,n); @ skif is the matrix of filtered probs @

skis = zeros(captst,n); @ skis is the matrix of smoothed probs @

id = eye(ns); @ used in certain calculations below @

"Bayesian prior used";

"a=";;a;;"b=";;b;;"c=";;c;

proc pattern1; @ This proc returns a (ps+1)*ns x n matrix. The ith

column contains a one in row j if st = j, contains a

one in row ns+j if st-1 = j, and so on @

local i1,ix,iq,na;

na = n/ns;

ix = eye(ns).*.ones(1,na);

i1 = 1;

do until i1 > ps;

na = na/ns;

iq = ones(1,ns^i1).*.(eye(ns).*.ones(1,na));

ix = iq|ix;

i1 = i1+1;

endo;

retp(ix);

endp;

hp = pattern1;

#include procs;

/* ================================================================= */

/* ======================================================================= */

proc matpm(xth); @This proc defines the user's conventions for reading

elements of Markov transition probabilities from

parameter vector @

local pm,ixth;

ixth = rows(xth);

pm = zeros(ns,ns);

if ipm == 1; @ for ns =2 this option has parameters as p11 and p22 @

if izz == 1;

pm[1,1] = xth[1,1]^2/(1 +xth[1,1]^2);

pm[2,2] = xth[2,1]^2/(1 + xth[2,1]^2);

else;

pm[1,1] = xth[1,1];

pm[2,2] = xth[2,1];

endif;

pm[2,1] = 1 - pm[1,1];

pm[1,2] = 1 - pm[2,2];

elseif ipm == 2; @ general case has parameters pij for i = 1,...,n and

j = 1,...,n-1 @

pm[1:ns-1,.] = reshape(xth[1:ixth,1],ns-1,ns);

if izz == 1;

pm[ns,.] = ones(1,ns);

pm = pm^2;

pm = pm./(sumc(pm)');

else;

pm[ns,.] = (1 - sumc(pm))';

endif;

elseif ipm == 3; @ This section can be rewritten by user to impose zeros

and ones where desired @

if izz == 1;

pm[1,1] = xth[1,1]^2/(1 + xth[1,1]^2);

pm[1,2] = xth[2,1]^2/(1 + xth[2,1]^2);

pm[2,3] = xth[3,1]^2/(1 + xth[3,1]^2);

elseif izz == 2;

pm[1,1] = xth[1,1];

pm[1,2] = xth[2,1];

pm[2,3] = xth[3,1];

endif;

pm[3,1] = 1 - pm[1,1];

pm[2,2] = 1 - pm[1,2];

pm[3,3] = 1 - pm[2,3];

endif;

retp(pm);

endp;

/* ======================================================================= */

/* ==================================================================== */

@ Set parameters to use Gauss numerical optimizer @

library optmum;

#include optmum.ext;

__btol = 1.e-06; @ This controls convergence criterion for coefficients@

__gtol = 1.e-06; @ This controls convergence criterion for gradient @

__algr = 1; @ This chooses BFGS optimization @

_opalgr = 2; @ This chooses BFGS optimization for newer versions of

GAUSS @

__miter = 150; @ This controls the maximum number of iterations @

__output = 1; @ This causes extra output to be displayed @

__covp = 0; @ This speeds up return from OPTMUM; note that the

program makes a reparameterization to calculate

std. errors @

@ Next call the GAUSS numerical optimizer @

output off;

{x,f,g,h} =optmum(&ofn,startval);

output file=junk on;

"";"";"MLE as parameterized for numerical optimization ";

"Coefficients:";x';

"";"Value of log likelihood:";;-f;

"";"Gradient vector:";g';

/* ======================================================================= */

/*======================================================================== */

@ In general no parts of this section need be changed @

@ Reparameterize for reporting final results @

izz = 2;

"";"Vector is reparameterized to report final results as follows";

"Means for each state:";x[1:ns,1]';

If pphi > 0;

"Autoregressive coefficients:";x[ns+1:ns+pphi,1]';

endif;

ncount = ns + pphi ;

x[ncount+1:ncount+isig,1] = x[ncount+1:ncount+isig,1]^2;

"Variances:";x[ncount+1:ncount+isig,1];

ncount = ncount + isig;

/* ======================================================================= */

/* ======================================================================= */

proc pmth; @ This proc converts the last elements of parameter vector

(which relates to transition probabilities) from the

th[i,j]^2/{sum j th[i,j]^2} form that is used for

numerical estimation into the p[i,j] form that is used

to calculate standard errors @

local pm;

if ipm == 1;

x[ncount+1,1] = x[ncount+1,1]^2/(1 + x[ncount+1,1]^2);

x[ncount+2,1] = x[ncount+2,1]^2/(1 + x[ncount+2,1]^2);

elseif ipm == 2;

pm = zeros(ns,ns);

pm[1:ns-1,.] = reshape(x[ncount+1:nth],ns-1,ns);

pm[ns,.] = ones(1,ns);

pm = pm^2;

pm = pm./(sumc(pm)');

x[ncount+1:nth,1] = reshape(pm[1:ns-1,.],ns*(ns-1),1);

elseif ipm == 3; @ User may want to alter these next lines @

x[ncount+1:nth,1]

= (x[ncount+1:nth,1]^2)./(1 + x[ncount+1:nth,1]^2);

endif;

retp(x);

endp;

/* ======================================================================= */

/* ======================================================================= */

@ In general no changes are necessary from here out @

call pmth;

h = (hessp(&ofn,x));

va = eigrs(h);

kc = 2;

ks = 2;

call ofn(x);

if minc(eigrs(h)) <= 0;

"Negative of Hessian is not positive definite";

"Either you have not found local maximum, or else estimates are up "

"against boundary condition. In latter case, impose the restricted "

"params rather than estimate them to calculate standard errors";

else;

h = invpd(h);

std = diag(h)^.5;

"For vector of coefficients parameterized as follows,";x';

"the standard errors are";std';

endif;

"";"-------------------------------";"";

"Probabilities for primitive states";

"filtered probabilities";format /rd 1,0;

"Obs ";;

t = 0;

do until t > ps;

i = 1;

do until i == ns;

"P(st-";;t;;"=";;i;;") ";;

i = i+1;

endo;

t = t+1;

endo;"";

format /rd 6,4;

skif = (skif*hp')*(eye(ps+1).*.id[.,1:ns-1]);

skif = seqa(nk,1,captst)~skif;skif;

"";"smoothed probabilities";

format /rd 1,0;

"Obs ";;

i = 1;

do until i > ns;

"P(st = ";;i;;") ";;

i = i+1;

endo;

format /rd 6,4;

skis = skis*hp';

skis = seqa(nk,1,captst)~skis[.,1:ns];skis;

/*========================================================================*/

PROCS

@ These GAUSS procedures evaluate Markov transition matrix, evaluate

likelihood function, and evaluate filter and smoothed probababilities @

/* ============================================================== */

proc matf(pm); @This proc returns the (n x n) matrix F of Markov

transition probabilities for state vector @

local iz,iw,ib,na,nb,nc,fz,fm;

@ set initial values for use with iteration @

na = 1;

nb = ns;

nc = ns*ns;

fm = pm;

iz = 1;

do until iz > ps;

fz = fm;

fm = zeros(nc,nc);

iw = 1;

do until iw > ns;

fm[((iw-1)*nb+1):(iw*nb),((iw-1)*na+1):(iw*na)]

= fz[1:nb,((iw-1)*na+1):iw*na];

iw = iw+1;

endo;

ib = 2;

do until ib > ns;

fm[1:nc,((ib-1)*nb+1):ib*nb] = fm[1:nc,1:nb];

ib = ib+1;

endo;

na = na*ns;

nb = nb*ns;

nc = nc*ns;

iz = iz+1;

endo;

retp(fm);

endp;

/* ==================================================================== */

/* ================================================================= */

proc ofn(th); @ this proc evaluates filter probs and likelihood @

local ncount,mu,phi,sig,pm,eta,iz,fm,chsi,it,f,fit,fx,hw,fj,

ij,ap,const,pq,hk,ihk;

@ Convert parameter vector to convenient form @

mu = th[1:ns,1];

ncount = ns + 1; @ ncount is the number of params read @

if pphi == 0;

phi = 1;

else;

phi = 1 | -th[ncount:ncount+pphi-1];

endif;

ncount = ncount + pphi;

if isig == 1;

sig = th[ncount,1]*ones(n,1);

ncount = ncount +1;

elseif isig == ns;

sig = th[ncount:ncount+ns-1,1].*hp[1:ns,.];

sig = sumc(sig);

ncount = ncount + ns;

endif;

if izz == 1;

sig = sig^2;

endif;

pm = matpm(th[ncount:nth,1]);

@ Construct constant term to be subtracted for each observation @

const = phi .*. mu;

const = const'*hp;

@ Convert data to AR resids @

eta = y[nk:capt,1];

iz = 1;

do until iz > pphi;

eta = eta ~y[nk-iz:capt-iz,1];

iz = iz+1;

endo;

eta = ((eta*phi - const)^2)./sig';

@ Adjust data to avoid numeric overflows and underflows

eta = reshape(eta,1,n*(captst));

eta = minc(eta|800*ones(1,n*captst));

eta = reshape(eta',captst,n); @

eta = (1./sqrt(sig')).*exp(-eta/2);

@ Calculate ergodic probabilities @

fm = matf(pm);

ap = (eye(n)-fm)|ones(1,n);

chsi = sumc((invpd(ap'*ap))');

chsi = maxc(chsi'|zeros(1,n)); @ This line eliminates roundoff error @

if kc > 1;

"";"Matrix of Markov transition probabilities:";pm;

"";"Ergodic probs for full state vector:";chsi';

"";"Ergodic probs for primitive states:";

pq = hp*chsi;pq[1:ns,1]';

endif;

@ Filter iteration @

f = 0;

it = 1;

do until it > captst;

fx = chsi .* eta[it,.]';

fit = sumc(fx);

skif[it,.] = fx'/fit;

f = f + ln(fit);

chsi = fm*fx/fit;

it = it+1;

endo;

@ Make Bayesian adjustment to likelihood if desired @

if a > 0;

fj = 0;

ij = 1;

do until ij > ns;

fj = fj + a*ln(sig[ij,1]) + b/sig[ij,1]

+ c*(mu[ij,1]^2/sig[ij,1]);

ij = ij +1;

endo;

f = f - fj/2;

endif;

@ Calculate smoothed probs if desired @

if ks == 2;

skis[captst,.] = skif[captst,.];

it = 1;

do until it == captst;

if minc(skif[captst-it,.]') > 1.e-150;

skis[captst-it,.] = skif[captst-it,.].*

((skis[captst-it+1,.]./(skif[captst-it,.]*fm'))*fm);

else; @ adjust code so as not to divide by zero @

hk = skif[captst-it,.]*fm';

ihk = 1;

do until ihk > n;

if hk[1,ihk] > 1.e-150;

hk[1,ihk] = skis[captst-it+1,ihk]/hk[1,ihk];

else;

hk[1,ihk] = 0;

endif;

ihk = ihk + 1;

endo;

skis[captst-it,.] = skif[captst-it,.].*(hk*fm);

endif;

it = it+1;

endo;

endif;

@ Print out value of log likelihood if desired @

if kc == 2;

"";"Log likelihood:";f;

endif;

retp(-f);

endp;

/*=========================================================================*/

I copied and pasted your code from above into two separate files, `markov_mod.gss` and `procs`. I had to do some find and replace, because the forum website changed some of the quotation marks to two single quotes, subtraction signs to dashes, etc. I also had to close the comment at the very last line of your post. As a note for others, I also did not copy the line that says MAXSEEK or the line that says PROCS into my files. I used those as demarcations between the two files.

After this, I downloaded the `gnp82.dat` dataset from the link that says "Programs for estimation of Markov switching models by numerical optimization" on this page. I placed the `gnp82.dat` file in the same directory as the other two files.

I was then able to run the code successfully with GAUSS version 13.1 on Linux. I have attached a zip file with these contents.

Final note: this code requires the GAUSS Optmum application module.