@ 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 spar,eps,alpha0,phi,sig0,a0,a1,g1,pm,b1,l1,cm,nu,mex,
    eta,eta2,heta,hsig,hq,lambada,hprob,xxx,serr,aerr,
    vyx,vyxx,yxx,ek,sk,nt,vex,msw,rho,vrho,ky,kq,q,lam,lamz,fl,fz,ip,vf,
    z00,a00,varf,iz,fm,chsi,it,f,fit,fx,hw,fj,ij,ap,const,pq,hk,ihk;


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

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

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

  if pphi > 0;
     phi = th[spar:spar+pphi-1,1];
  else;
     phi = 0;
  endif;
  spar = spar + pphi;
  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;
  a0 = th[spar,1];
  else;
  a0 = abs(th[spar,1]);
  endif;
  if karch >0;
  if earch ==1;
  a1 = th[spar+1:spar+karch,1];
    else;
  a1 = abs(th[spar+1:spar+karch,1]);
    endif;
  else;
  a1=0;
    endif;
  spar = spar+1+karch;

  pm = matpm(th[spar:spar+(ns*(ns-1))-1-rest,1]);
  spar = spar+(ns*(ns-1))-rest;
  g1 = 1|(abs(th[spar:spar+ns-2,1]));
@  g1 = 1|(th[spar,1]|th[spar+1,1]*100);@
  spar = spar+ns-1;
  if larch == 1;
    if earch == 1;
     l1 = th[spar,1];
    else;
     l1 = abs(th[spar,1]);
    endif;
     spar = spar+1;
  else;
  l1=0;
  endif;

  if vexog > 0;
     vex = abs(th[spar:spar+vexog-1,1]);
  else;
     vex = 0;
  endif;
  spar = spar + vexog;

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

  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 ------------------ @
  resid = y[pphi+1:capt,1] - alpha0;

    if mexog==1;
       resid = resid - xx[pphi+1:capt,.]*mex;
    endif;

    iz = 1;
    do until iz > pphi;
       resid = resid - phi[iz,1]*y[pphi+1-iz:capt-iz,1];
    iz = iz+1;
    endo;
    if pphi > 0;
       resid = (sig0*ones(pphi,1)) | resid;
    endif;

@------------- Collect karch+1 most recent squared residuals in
               the captst by karch + 1 matrix heta ------@
    eta2 = resid^2;
    if vexog >0;
    heta = a0*ones(captst,1) + vx[nk-wlag:capt-wlag,vexog]*vex;
    else;
    heta = a0*ones(captst,1);
    endif;

/*    heta = a0*ones(captst,1); */
   if karch == 0;
   heta = a0*ones(captst,1)~zeros(captst,1);
   endif;

      iz = 1;
      do until iz > karch;
        if earch == 1;
         heta = heta ~ abs(resid[nk-iz:capt-iz,1]);
        else;
         heta = heta ~ eta2[nk-iz:capt-iz,1];
        endif;
      iz = iz + 1;
      endo;

     if larch == 1;  @ Correct heta for leverage effects if desired @
        if earch == 1;
          heta[.,2] = heta[.,2]
                + (l1/a1[1,1])*resid[nk-1:nk-1+captst-1,1];
        else;
          heta[.,2] = heta[.,2]
                + (l1/a1[1,1])*
@                 ((resid[nk-1:nk-1+captst-1,1] .< 0).*   GJR @
                 ((resid[nk-1:nk-1+captst-1,1] .> 0).*
                   eta2[nk-1:nk-1+captst-1,1]);
        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 ---------@
hsig = zeros(ps+1,(ps+1)*ns);
hsig[1,1:ns] = ones(1,ns);
iz = 2;
do until iz > ps+1;
   hsig[iz,(iz-1)*ns+1:iz*ns] = a1[iz-1,1]./g1';
iz = iz+1;
endo;
hsig = hsig*hp;

@----------- 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 ----------@
heta = heta*hsig;
/*
"karch1:";;karch1;
"ns:";;ns;
"rows(ones(1,ns^karch1)):";;rows(ones(1,ns^karch1));
"cols(ones(1,ns^karch1)):";;cols(ones(1,ns^karch1));
"g1:";;g1';

"rheta:";;rows(heta);
"cheta:";;cols(heta);
*/
 hq = ones(1,ns^karch1).*.g1';
    if earch == 1;
     heta=exp(heta);
    endif;
/*
"rhq:";;rows(hq);
"chq:";;cols(hq);
*/

heta = heta.*hq;
/*
    if vexog >0;
    heta = heta + vx[nk-wlag:capt-wlag,vexog]*vex;
    endif;
*/
@---------- 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;

if mswarch==1;
    eta2[nk:nk+captst-1,1] = (resid[nk:nk+captst-1,1]-heta*chsi*msw)^2;
endif;

@------------ 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 ---@
    if tarch == 0;    @ normal density here @
        eta = eta2[nk:nk+captst-1,1]./heta;
        eta = 0.39894228 * (1./sqrt(heta)).* exp(-eta/2);
    elseif tarch == 1; @ t density here @
         eta = ones(captst,n) + (eta2[nk:nk+captst-1,1]./(heta*(nu-2)));
         eta = eta^(-(nu+1)/2)*lambada./sqrt(heta);
    endif;

@ ------------The basic filter iteration is contained here ----------- @
   f = 0;
   it = 1;
   do until it > captst;

       if ks == 2;  @calculate predicted variances and outlier probs if
                     desired @
          varfor[it+nk-1,1] = heta[it,.]*chsi;
          if tarch == 0;
             hprob = cdfnc(resid[it+nk-1,1]./sqrt(heta[it,.]));
          elseif tarch == 1;
             hprob = cdftc(resid[it+nk-1,1]./
             sqrt(heta[it,.]*(nu-2)/nu),nu);
          endif;
          outprob[it+nk-1,1] = hprob*chsi;
        endif;

     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;

@ --------- 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 ----------@

if kc==2;
  "log likelihood= ";; f;
  "Variance= ";;(resid[nk:capt,1]'*resid[nk:capt,1])/(captst-1);
  "Residuals Sum of Squares= ";;(resid[nk:capt,1]'*resid[nk:capt,1]);
  "Mean of sqrt(h)= ";;meanc(sqrt(varfor[nk:capt,1]));
  "Mean of h= ";;meanc(varfor[nk:capt,1]);
  "";
  "Analisis of residuals";

yxx=resid[nk:capt,1]./sqrt(varfor[nk:capt,1]);
"mean = ";;meanc(yxx);
vyx = ((yxx-meanc(yxx))'*(yxx-meanc(yxx)));
vyxx = ((yxx-meanc(yxx))'*(yxx-meanc(yxx)))/(captst-1);
"st. deviation = ";;sqrt(vyxx);
sk = (ones(captst,1)'*(yxx-meanc(yxx))^3)/captst;
sk=sk/((vyx/captst)^1.5);
"SK= ";;sk;
ek = (ones(captst,1)'*(yxx-meanc(yxx))^4)/captst;
ek=(ek/((vyx/captst)^2))-3;
"EK= ";;ek;
nt=(captst/6)*(sk^2+.25*ek^2);
"Normality Test -Jarque-Bera (1980):";;nt;
"";
"Sample autocorrelations for levels";
"kau:";;kau;
rho=zeros(kau,1);

kx=0;
vrho=(yxx[1:captst-kx,1]-meanc(yxx))'*(yxx[kx+1:captst,1]-meanc(yxx))/captst;

kx=1;
do until kx > kau;
rho[kx,1]=(yxx[1:captst-kx,1]-meanc(yxx))'*(yxx[kx+1:captst,1]-meanc(yxx));
rho[kx,1]=rho[kx,1]/(captst*vrho);
kx=kx+1;
endo;
"rho:";
rho';
"";
"Standard Error:";;sqrt(1/captst);
q=0;
kq=1;
do until kq > (kau/2);
q = q + rho[kq,1]^2/(captst-kq);
kq=kq+1;
endo;
"Ljung-Box-(kau/2):";;q*captst*(captst+2);
kq=1;
q=0;
do until kq > kau;
q = q + rho[kq,1]^2/(captst-kq);
kq=kq+1;
endo;
"Ljung-Box-kau:";;q*captst*(captst+2);
"";
"Sample autocorrelations for squared residuals";
yxx=yxx^2;
rho=zeros(kau,1);

kx=0;
vrho=(yxx[1:captst-kx,1]-meanc(yxx))'*(yxx[kx+1:captst,1]-meanc(yxx))/captst;

kx=1;
do until kx > kau;
rho[kx,1]=(yxx[1:captst-kx,1]-meanc(yxx))'*(yxx[kx+1:captst,1]-meanc(yxx));
rho[kx,1]=rho[kx,1]/(captst*vrho);
kx=kx+1;
endo;
"rho:";
rho';
"";
"Standard Error:";;sqrt(1/captst);
q=0;
kq=1;
do until kq > (kau/2);
q = q + rho[kq,1]^2/(captst-kq);
kq=kq+1;
endo;
"Ljung-Box-(kau/2):";;q*captst*(captst+2);
kq=1;
q=0;
do until kq > kau;
q = q + rho[kq,1]^2/(captst-kq);
kq=kq+1;
endo;
"Ljung-Box-kau:";;q*captst*(captst+2);
"";

if mk>0;
"Nstep Forecast of Variance";
"";
z00 = zeros(mk,1);

  z00 = y[capt+1:capt+mk,1] - alpha0 ;

    if mexog==1;
       z00 = z00 - xx[capt+1:capt+mk,1]*mex;
    endif;

    iz = 1;
    do until iz > pphi;
       z00 = z00 - phi[iz,1]*y[capt:capt+mk-1,1];
    iz = iz+1;
    endo;

@----------- Calculate fl (the mk-period-ahead transition (n x n) matrix)
      and lam (the 2 x 1 mk-period ahead decay factor for ARCH component)----@
    lam = zeros(karch1,karch1);
    lam[1,.] = a1';
    if larch == 1;
       lam[1,1] = lam[1,1] + .5*l1;
    endif;
    a00 = a0;
    if karch > 1;
      lam[2:karch1,1:karch1-1] = eye(karch1-1);
      a00 = a0|zeros(karch1-1,1);
    endif;
    "lam is";lam;
"";
"Last forecasted Variance:";;varfor[capt,1];
"";
    fl = eye(n);
    it = 1;
    do until it >= mk;
      fl = fl*fm;
    it = it + 1;
    endo;
    hp = pattern1;
    lamz = zeros((ps+1)*ns,karch1);
    ip = 1;
    do until ip > karch1;
       lamz[ip*ns+1:(ip+1)*ns,ip] = g1;
    ip = ip+1;
    endo;
    lamz = lamz'*hp;

@----------------- calculate predicted variances ---------------------@
          varf=zeros(mk,1);
          fz = (eta2[captst+nk-1,1]*ones(1,n)) ./ lamz[1,.];
               if karch > 1;
                  ip = 2;
                  do until ip > karch;
                   fz = fz | ((eta2[captst+nk-ip,1]*ones(1,n)) ./ lamz[ip,.]);
                  ip = ip+1;
                  endo;
               endif;
          vf = a0 + a1'*fz;
          if resid[captst+nk-1,1] < 0;
             vf = vf + l1*((eta2[captst+nk-1,1]*ones(1,n)) ./ lamz[1,.]) ;
          endif;

               if karch > 1;
                  fz[2:karch,.] = fz[1:karch-1,.];
               endif;
               fz[1,.] = vf;

               ip=1;
               do until ip > mk;
                fl = fl*fm;
                  if ip==1;
                    vf = hq .* vf;
                  else;
                    fz = a00 + lam*fz;
                    vf = (hq*fl) .* fz[1,.];
                  endif;
                varf[ip,1] = vf*chsi;
   if captx>capt;
     serr=(z00[ip,1]^2-varf[ip,1])^2;
     aerr=abs(z00[ip,1]^2-varf[ip,1]);
     ip~varf[ip,1]~serr~aerr;
   else;
    ip~varf[ip,1];
   endif;
               ip = ip+1;
               endo;
endif;

endif;

retp(-f);
endp;