@ These GAUSS procedures evaluate Markov transition matrix and evaluate bivariate
   likelihood function @

/* ============================================================== */
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 spar,eps,alpha0,phi1,phi2,sig0,a0,a1,g1,pm,b1,l1,cm,nu,
    eta,eta2,heta,hsig,hq,lambada,hprob,a01,a02,a11,a12,g2,l2,alpha02,
    vyx,vyxx,yxx,ek,sk,nt,resid1,resid2,lama,heta1,heta2,kp,p,tj,tt,j,i,
    kl,iz,fm,chsi,it,f,fit,fx,hw,fj,ij,ap,const,pq,hk,ihk,corr,hsig1,hsig2,
    hq1,hq2,lhh11,lhh12,lhh21,lhh22,lhh1,lhh2,hhw11,hhw12,hhw21,hhw22,
    cl,mex1,mex2,yxy,f1,f2,ff,iota,t1,t2,t3,t4,t5,sr,wv,iv,varv;

  if ismiss(th);
        retp(error(0));
  endif;

@ -------------- Convert parameter vector to convenient form --------------@
  spar = 2;
  alpha0 = th[1,1];
  alpha02 = th[2,1];
  spar = spar + 1;

  if mexog1 > 0;
     mex1 = th[spar:spar+mexog1-1,1];
  else;
     mex1 = 0;
  endif;
  spar = spar + mexog1;

  if mexog2 > 0;
     mex2 = th[spar:spar+mexog2-1,1];
  else;
     mex2 = 0;
  endif;
  spar = spar + mexog2;

  if pphi1 > 0;
     phi1 = th[spar:spar+pphi1-1,1];
  else;
     phi1 = 0;
  endif;
  spar = spar + pphi1;

  if pphi2 > 0;
     phi2 = th[spar:spar+pphi2-1,1];
  else;
     phi2 = 0;
  endif;
  spar = spar + pphi2;

  if barch == 0;
      sig0 = abs(th[spar,1]);
      spar = spar + 1;
  else;
      sig0 = 10000;   @ this should not be used with this option.  The absurd
                          large value should flag if used inadvertently @
  endif;

  if earch ==1;
  a01 = th[spar,1];
    else;
  a01 = abs(th[spar,1]);
    endif;
  spar = spar+1;

  if karch1 >0;
    if earch ==1;
      a11 = th[spar:spar+karch1-1,1];
    else;
      a11 = abs(th[spar:spar+karch1-1,1]);
    endif;
   else;
     a11=0;
  endif;
  spar = spar+karch1;

  if earch ==1;
  a02 = th[spar,1];
    else;
  a02 = abs(th[spar,1]);
    endif;
  spar = spar+1;

  if karch2 >0;
    if earch ==1;
     a12 = th[spar:spar+karch2-1,1];
    else;
     a12 = abs(th[spar:spar+karch2-1,1]);
    endif;
  else;
    a12=0;
  endif;
  spar = spar+karch2;

  pm = matpm(th[spar:spar+(ns*(ns-1))-1-rest,1]);
  spar = spar+(ns*(ns-1))-rest;
  g1 = 1|abs(th[spar:spar+ns1-2,1]);
  spar = spar+(ns1-1);
  g2 = 1|abs(th[spar:spar+ns2-2,1]);
  spar = spar+(ns2-1);

  if larch1 == 1;
    if earch == 1;
     l1 = th[spar,1];
    else;
     l1 = abs(th[spar,1]);
    endif;
     spar = spar+1;
  endif;

  if larch2 == 1;
    if earch == 1;
     l2 = th[spar,1];
    else;
     l2 = abs(th[spar,1]);
    endif;
     spar = spar+1;
  endif;

   if cov > 0;
@    corr = abs(th[spar:spar+cov-1,1]);@
    corr = 1/(1+th[spar:spar+cov-1,1]^2);
    spar = spar + cov;
   else;
    corr = 0;
   endif;

  if tarch == 1;
       nu = 2 + abs(th[spar,1]);
       lambada = (1/sqrt(((nu - 2)*pi)))*gamma((nu+1)/2)/gamma(nu/2);
       spar = spar + 1;
  endif;

@ ---------- Convert data to AR resids ----------@
if mexog1 >0;
  resid1 = y1[pphi1+1:capt,1] - alpha0 - xx[pphi1+1:capt,1]*mex1;
else;
  resid1 = y1[pphi1+1:capt,1] - alpha0;
endif;

if mexog2 >0;
  resid2 = y2[pphi2+1:capt,1] - alpha02 - xx[pphi2+1:capt,.]*mex2;
else;
  resid2 = y2[pphi2+1:capt,1] - alpha02;
endif;

    iz = 1;
    do until iz > pphi1;
       resid1 = resid1 - phi1[iz,1]*y1[pphi1+1-iz:capt-iz,1];
       resid2 = resid2 - phi2[iz,1]*y2[pphi2+1-iz:capt-iz,1];
    iz = iz+1;
    endo;
    if pphi1 > 0;
       resid1 = (sig0*ones(pphi1,1)) | resid1;
       resid2 = (sig0*ones(pphi2,1)) | resid2;
    endif;

@------------- Collect karch+1 most recent squared residuals in
               the captst by karch + 1 matrix heta ------@
    resid=resid1~resid2;
    eta2 = resid^2;
    heta1 = a01*ones(captst,1);
    heta2 = a02*ones(captst,1);
      iz = 1;
    if karch1 >0;
      do until iz > karch1;
        if earch == 1;
         heta1 = heta1 ~ abs(resid1[nk-iz:capt-iz,1]);
        else;
         heta1 = heta1 ~ eta2[nk-iz:capt-iz,1];
        endif;
      iz = iz + 1;
      endo;
    else;
         heta1 = heta1 ~ zeros(captst,1);
         a11=1;
    endif;

      iz=1;
    if karch2 >0;
      do until iz > karch2;
        if earch == 1;
         heta2 = heta2 ~ abs(resid2[nk-iz:capt-iz,1]);
        else;
         heta2 = heta2 ~ eta2[nk-iz:capt-iz,2];
        endif;
      iz = iz + 1;
      endo;
    else;
         heta2 = heta2 ~ zeros(captst,1);
         a12=1;
    endif;



     if larch1 == 1;  @ Correct heta for leverage effects if desired @
        if earch == 1;
          heta1[.,2] = heta1[.,2]
                + (l1/a11[1,1])*resid1[nk-1:nk-1+captst-1,1];
        else;
          heta1[.,2] = heta1[.,2]
                + (l1/a11[1,1])*
                 ((resid1[nk-1:nk-1+captst-1,1] .< 0).*
                   eta2[nk-1:nk-1+captst-1,1]);
        endif;
      endif;

     if larch2 == 1;  @ Correct heta for leverage effects if desired @
        if earch == 1;
          heta2[.,2] = heta2[.,2]
                + (l2/a12[1,1])*resid2[nk-1:nk-1+captst-1,1];
        else;
          heta2[.,2] = heta2[.,2]
                + (l2/a12[1,1])*
                 ((resid2[nk-1:nk-1+captst-1,1] .< 0).*
                   eta2[nk-1:nk-1+captst-1,2]);
        endif;
    endif;

@---------------- Construct karch + 1 by n matrix hsig that captures
           correct weights on lagged squared residuals for each of the
           n possible configurations for most recent states ---------@
hsig1 = zeros(ps+1,(ps+1)*ns1);
hsig2 = zeros(ps+1,(ps+1)*ns2);
hsig1[1,1:ns1] = ones(1,ns1);
hsig2[1,1:ns2] = ones(1,ns2);
iz = 2;
do until iz > ps+1;
   hsig1[iz,(iz-1)*ns1+1:iz*ns1] = a11[iz-1,1]./g1';
   hsig2[iz,(iz-1)*ns2+1:iz*ns2] = a12[iz-1,1]./g2';
iz = iz+1;
endo;

hsig1 = hsig1*hp1;
hsig2 = hsig2*hp2;

@----------- Postmultiply heta by hsig to obtain the captst by n matrix of
        variances that would apply for each of the n possible configurations
        for recent states ----------@

heta1 = heta1*hsig1;

heta2 = heta2*hsig2;

    if earch == 1;
     heta1=exp(heta1);
     heta2=exp(heta2);
    endif;

 hq1 = ones(1,ns1^kzarch1).*.g1';
 hq2 = ones(1,ns2^kzarch2).*.g2';

heta1 = heta1.*hq1;
heta2 = heta2.*hq2;

@---------- Calculate ergodic probabilities -----------@
     fm = matf(pm);
     ap = (eye(n)-fm)|ones(1,n);
ap; 
    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;
@ ------------The basic filter iteration is contained here ----------- @

@------------ Calculate a captst by n matrix of conditional densities;
   row t of this matrix gives the n possible values the likelihood would
   take at date t for each of the n possible configurations of states ---@
lama=zeros(kk,kk);
eta=zeros(captst,16);
   f = 0;
   it = 1;
 do until it > captst;

hhw11=heta1[it,.]*hw1;
hhw12=heta1[it,.]*hww1[1:4,.];

hhw21=heta2[it,.]*hw2;
hhw22=heta2[it,.]*hww2[1:4,.];

lhh11=hhw11'.*hh[.,1:4];
lhh12=hhw12'.*hh[.,1:4];
lhh1=lhh11~lhh12;

lhh21=hhw21'.*hh[.,1:2];
lhh22=hhw22'.*hh[.,3:4];
lhh2=lhh21~lhh22~lhh21~lhh22;
lama=lhh1+lhh2;

j=1;
p=1;
    varv=zeros(nm,nm);
if cov==0;
  do until j >kk;
   i=1;
     do until i >kk;
eta[it,p]= exp(-.5*resid[nk+it-1,.]*inv(lama[i:i+1,j:j+1])*resid[nk+it-1,.]');
 eta[it,p]=eta[it,p]/sqrt(lama[i,j]*lama[i+1,j+1]);
  varv=varv+lama[i:i+1,j:j+1]*chsi[p,1];
      i=i+2;
      p=p+1;
    endo;
   j=j+2;
  endo;
else;
  do until j >kk;
   i=1;
   cl=1;
     do until i >kk;
      lama[i+1,j]=corr[cl,1]*sqrt(lama[i,j]*lama[i+1,j+1]);
      lama[i,j+1]=lama[i+1,j];
  varv=varv+lama[i:i+1,j:j+1]*chsi[p,1];
eta[it,p]= exp(-.5*resid[nk+it-1,.]*inv(lama[i:i+1,j:j+1])*resid[nk+it-1,.]');
  eta[it,p]=eta[it,p]/sqrt(det(lama[i:i+1,j:j+1]));
      i=i+2;
      p=p+1;
        if i==kk/2+1;
        cl=cov;
        endif;
    endo;
   j=j+2;
  endo;

  if kvar==1;
       output file=junkw.var on;
    f1 = alpha0 + phi1[1,1]*y1[nk+it-2,1];
    f2 = alpha02 + phi2[1,1]*y2[nk+it-2,1];
    ff=f1|f2;
    iota=ones(nm,1);
    iv=inv(varv);
    t1=(iota'*iv*iota);
    t3=iv*iota/(iota'*iv*iota);
    t2=(iota'*iv*ff);

    wv=((1/lam)*iv*ff)+((1- t2/lam)/t1)*iv*iota;

    t4=wv'*ff;
    t5=wv'*varv*wv;
    sr=t4/sqrt(t5);

    it~wv'~sr;
       output file=junkw.var off;
  endif;

  if kvar==2;
       output file=junkw.var on;
    f1 = alpha0 + phi1[1,1]*y1[nk+it-2,1];
    f2 = alpha02 + phi2[1,1]*y2[nk+it-2,1];
    ff=f1|f2;
    iota=ones(nm,1);
    iv=inv(varv);
    t1=(iota'*iv*iota);
    t3=iv*iota/(iota'*iv*iota);
    wv=t3;
    t4=wv'*ff;
    t5=wv'*varv*wv;
    sr=t4/sqrt(t5);

    it~wv'~sr;
       output file=junkw.var off;
  endif;

endif;
     fx = chsi .* eta[it,.]';
     fit = sumc(fx);
     skif[it,.] = fx'/fit;
     f = f + ln(fit);
     chsi = fm*fx/fit;

   it = it+1;
  endo;

@ --------- Calculate smoothed  probabilities 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 ----------@

f= f - captst*nm*0.918938534;

if kc==2;
  "log likelihood= "; f;
endif;
retp(-f);
endp;