@ 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;
    if fit > 0;
     f = f + ln(fit);
    else;
     f = f - 10;
    endif;
     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;
/*=========================================================================*/