/*===================================================================
  AN APPLICATION OF   "DYNAMIC LINEAR MODELS WITH MARKOV-SWITCHING"
  (KIM, CHANG-JIN (1994), JOURNAL OF ECONOMETRICS, 60, 1-22.)
  by  Myung-Jig Kim and Ji-Sung Yoo (1995, Journal of Monetary Economics)
  and Marcelle Chauvet (1998, International Economic Review)
  : Stock and Watson Model with Markov-Switching.


  WRITTEN BY CHANG-JIN KIM,
  DEPT. OF ECONOMICS
  KOREA UNIVERSITY,
  SEOUL, 136-701, KOREA
  E-MAIL: cjkim@korea.ac.kr


  MODEL: (based on demeaned log differences)

          yt = gamma_i*(Delta c_t) + z_{it},  i=1,2,3; t=1,...,T
          y_t =gamma_{40}(Delta c_t)+gamma_{41}(Delta c_{t-1})+
           +gamma_{42}(Delta c_{t-2})+gamma_{43}(Delta c_{t-3})+z_{4t},

         (Delta c_t)= mu_{s_t}+ phi1* (Delta c_{t-1})+(Delta c_{t-2}) + e_t
         z_{it}=psi_{i1}*z_{i,t-1} +psi_{i2}*z_{i,t-2} + v_{it}
         e_t -- iid N(0,1)
         v_{it}--iid N(0,sig_i^2)
         mu_{s_t} =(1-S_t)*mu_0 + mu_1*S_t

         Pr[S_t=1|S_{t-1}=1]=p,  Pr[S_t=0|S_{t-1}=0]=q,

         PARAMETERS TO BE ESTIMATED:
           phi1,phi2, mu_0,mu_1,
           psi_{i1},psi_{i2}, sig_i^2, i=1,2,3,4
       gamma_i, i=1,2,3; gamma_{40},gamma_{41},gamma_{42},gamma_{43},
       p,q
=======================================================================*/


new;
library optmum,PGRAPH;

format /m1 /rd 8,4;
load data[433,6]=sw_ms.prn;  @1959.1 -- 1995.1@
                @ month, ip, gmyxpq, mtq, lpnag, DOC index@

month=data[2:433,1];

y1=(ln(data[2:433,2]) - ln(data[1:432,2]))*100;
y2=(ln(data[2:433,3]) - ln(data[1:432,3]))*100;
y3=(ln(data[2:433,4]) - ln(data[1:432,4]))*100;
y4=(ln(data[2:433,5]) - ln(data[1:432,5]))*100;

doc_in= (data[2:433,6]);
doc_din= ((data[2:433,6])-(data[1:432,6]));
doc_mn=meanc(doc_din);
doc_sd=stdc(doc_din);

yy0=(y1~y2~y3~y4);

/*
@================================Option 1: Based on====================
                             Standardized and Demeaned data
                            [Easier to get convergence!!]
=======================================================================@
yy_s=stdc(yy0);
yy_org= (yy0'./yy_s);
yy_m=meanc(yy_org');
yy=(yy_org - yy_m);  @Standardized y1, y2,y3,y4:  @
                          @NOTE:  DIMENSION= 4xT  @

prmtr_in={0.5  0.1  0.5 0.1 0.5 0.1 0.5 0.1 0.5  0.1
          0.5 0.5  0.5 0.5
          0.5 0.5 0.5 0.5 0 0 0   -1.8 0.5 2 3              };
@For optimal initial prameters, estimate S&W model without switchig first,
 then............@
@=======================================================================@
*/

@================================Option 2: Based on======================
                     Demeaned data (without standardizing)
             [Difficult to get convergence if starting values are
             not good enough!!!]
========================================================================@

yy_m=meanc(yy0);
yy=yy0'-yy_m;   @4xT@

prmtr_in={
  0.4089   0.2092  -0.0085  -0.0085  -0.1904  -0.1904  -0.2204  -0.2204
  0.8664  -1.3182   0.5082   0.5535   0.8029   0.1360   0.5585   0.2151
  0.4468   0.1223   0.1136   0.1096   0.4063  -1.4985   0.2620   1.6654
  3.5777
 };


prmtr_in={
  0.2091   0.2091  -0.0085  -0.0085  -0.1904  -0.1904  -0.2204  -0.2204
  0.8667  -1.3184   0.5082   0.5535   0.8029   0.1360   0.5584   0.2151
  0.4468   0.1223   0.1136   0.1096   0.4063  -1.4986   0.2619   1.6649
  3.5774 };

@For initial prameters, estimate S&W model without switchig first,
 then............@

@=====================================================================@

t=cols(yy);

/*
output file=sw_ms.out reset;
output off;

@ Maximum Likelihood Estimation @
@==================================================@


START=1;
prmtr_in=prmtr_in';

{xout,fout,gout,cout}=optmum(&lik_fcn,PRMTR_IN);

PRM_FNL=TRANS(xout);             @ Estimated  coefficients, constrained@
output on;

"==FINAL OUTPUT========================================================";

"likelihood value is ";; -1*fout;
"Estimated parameters are:";
prm_fnl';
output off;

"Calculating Hessian..... Please be patient!!!!";
                  hout0=hessp(&lik_fcn,xout);
                  hout=inv(hout0);

grdn_fnl=gradfd(&TRANS,xout);
Hsn_fnl=grdn_fnl*hout*grdn_fnl';
SD_fnl =sqrt(diag(Hsn_fnl)); @Standard errors of the estimated coefficients@

output on;

"--------------------------------------------------";
"log likelihood value is ";; -fout;
"-------------------------------------------------- ";

"Estimated parameters and standard errors are:";
" ";
prm_fnl~sd_fnl;
" ";
"---------------------------------------------------";
"Xout =";
xout';
" ";
"For Normal Convergence, The following number is 0:";
cout;
"---------------------------------------------------";
output off;
*/ xout=prmtr_in';

start=1;


{DLT_CTT,PRTT0,PRTL0} =FILTER(XOUT);
smooth0=SMOOTH(prtt0,prtl0); @Pr[S_t=0|Y_T], Smoothed probabilities@

xy(seqa(1,1,rows(prtt0)),prtt0);
xy(seqa(1,1,rows(smooth0)),smooth0);


sd_ctt=stdc(dlt_ctt);  @standard dev. of (Delta c_{t|t})@

new_ctt=dlt_ctt.*doc_sd./sd_ctt;

        new_ind=zeros(t,1);
        new_ind[1,1]=doc_in[1,1];
        ii=2;
        do until ii>t;
           new_ind[ii,1]=new_ctt[ii,1] + new_ind[ii-1,1] +doc_mn;
        ii=ii+1;
        endo;

new_ind=new_ind -new_ind[120,1]+doc_in[120,1];

output file=sw_ms.dta reset;
month~new_ind~doc_in~prtt0~smooth0;
output off;

xy(seqa(1,1,rows(new_ind[START:T])),new_ind[START:T]~doc_in[START:T]);
xy(seqa(1,1,rows(new_ind)),new_ind);


END;

@ END OF MAIN PROGRAM @
@========================================================================@
@========================================================================@


PROC LIK_FCN(PRMTR1);

LOCAL PRMTR, PPR,QPR, PHI1,PHI2, VAR_W, MU0,MU1,R,Q,H,F,B_LL0,B_LL1,
      VECP_LL,P_LL0,P_LL1,PROB__1,PROB__0,LIKV,J_ITER,B_TL00,B_TL01,
      B_TL10,B_TL11,P_TL00,P_TL01,P_TL10,P_TL11,F_CAST00,F_CAST01,
      F_CAST10,F_CAST11,SS00,SS01,SS10,SS11,B_TT00,B_TT01,B_TT10,
      B_TT11,P_TT00,P_TT01,P_TT10,P_TT11,PR_VL00,PR_VL01,PR_VL10,PR_VL11,
      PR_VAL,PRO_00,PRO_01,PRO_10,PRO_11,MU_M0,MU_M1,
      PSI11,PSI12,PSI21,PSI22,PSI31,PSI32,PSI41,PSI42,
      SIG1,SIG2,SIG3,SIG4,GAMMA1,GAMMA2,GAMMA3,GAMMA4,SIG_E,gamma41,
      gamma42,gamma43;

EXTERNAL PROC TRANS, V_PROB;

           PRMTR=TRANS(PRMTR1);

           LOCATE 20,1;  PRMTR';

    phi1=prmtr[1,1];
    phi2=prmtr[2,1];
    psi11=prmtr[3,1];
    psi12=prmtr[4,1];
    psi21=prmtr[5,1];
    psi22=prmtr[6,1];
    psi31=prmtr[7,1];
    psi32=prmtr[8,1];
    psi41=prmtr[9,1];
    psi42=prmtr[10,1];

    sig1=prmtr[11,1];
    sig2=prmtr[12,1];
    sig3=prmtr[13,1];
    sig4=prmtr[14,1];

    gamma1=prmtr[15,1];
    gamma2=prmtr[16,1];
    gamma3=prmtr[17,1];
    gamma4=prmtr[18,1];

    gamma41=prmtr[19,1];
    gamma42=prmtr[20,1];
    gamma43=prmtr[21,1];

    MU0=PRMTR[22,1];
    MU1=PRMTR[23,1];

    QPR=PRMTR[24,1];    @Pr[St=0/St-1=0]@
    PPR=PRMTR[25,1];    @Pr[St=1/St-1=1]@

@=====================================================================@

 F=(phi1~phi2~  0~  0~    0~   0~   0~   0~  0~  0~  0~  0)|
   (   1~   0~  0~  0~    0~   0~   0~   0~  0~  0~  0~  0)|
   (   0~   1~  0~  0~    0~   0~   0~   0~  0~  0~  0~  0)|
   (   0~   0~  1~  0~    0~   0~   0~   0~  0~  0~  0~  0)|
   (   0~   0~  0~  0~psi11~psi12~  0~   0~  0~  0~  0~  0)|
   (   0~   0~  0~  0~    1~   0~   0~   0~  0~  0~  0~  0)|
   (   0~   0~  0~  0~    0~   0~psi21~psi22~  0~  0~  0~  0)|
   (   0~   0~  0~  0~    0~   0~   1~   0~  0~  0~  0~  0)|
   (   0~   0~  0~  0~    0~   0~   0~   0~psi31~psi32~  0~  0)|
   (   0~   0~  0~  0~    0~   0~   0~   0~  1~  0~  0~  0)|
   (   0~   0~  0~  0~    0~   0~   0~   0~  0~  0~psi41~psi42)|
   (   0~   0~  0~  0~    0~   0~   0~   0~  0~  0~  1~  0);

   H=(gamma1~0~0~0~                  1~0~0~0~0~0~0~0)|
     (gamma2~0~0~0~                  0~0~1~0~0~0~0~0)|
     (gamma3~0~0~0~                  0~0~0~0~1~0~0~0)|
     (gamma4~gamma41~gamma42~gamma43~0~0~0~0~0~0~1~0);

    sig_e=1;

    Q=(sig_e^2~0~0~0~    0~0~    0~0~    0~0~    0~0)|
      (   0~   0~0~0~    0~0~    0~0~    0~0~    0~0)|
      (   0~   0~0~0~    0~0~    0~0~    0~0~    0~0)|
      (   0~   0~0~0~    0~0~    0~0~    0~0~    0~0)|
      (   0~0~0~0~  sig1^2~0~    0~0~    0~0~    0~0)|
      (   0~0~0~0~       0~0~    0~0~    0~0~    0~0)|
      (   0~0~0~0~       0~0~sig2^2~0~   0~0~    0~0)|
      (   0~0~0~0~       0~0~    0~0~    0~0~    0~0)|
      (   0~0~0~0~       0~0~    0~0~sig3^2~0~   0~0)|
      (   0~0~0~0~       0~0~    0~0~    0~0~    0~0)|
      (   0~0~0~0~       0~0~    0~0~    0~0~sig4^2~0)|
      (   0~0~0~0~       0~0~    0~0~    0~0~    0~0);

    R=0;

     MU_M0=MU0|ZEROS(11,1);
     MU_M1=MU1|ZEROS(11,1);

/*  NOTE:  KIM'S FILTER = HAMILTON FILTER + KALMAN FILTER               */
/*  INITIALIZING THE FILTER: FOR BOTH KALMAN FILTER AND HAMILTON FILTER */

         B_LL0= INV(EYE(12)-F)* MU_M0;
         B_LL1= INV(EYE(12)-F)* MU_M1;

         VECP_LL=INV(EYE(144) - F.*.F)*vec(Q);

    P_LL0=vecP_LL[1:12,1]~vecP_LL[13:24,1]~vecP_LL[25:36,1]~vecP_LL[37:48,1]
        ~vecP_LL[49:60,1]~vecP_LL[61:72,1]~vecP_LL[73:84,1]~vecP_LL[85:96,1]
        ~vecP_LL[97:108,1]~vecP_LL[109:120,1]~vecP_LL[121:132,1]~
         vecP_LL[133:144,1];

    P_LL1=P_LL0;

         PROB__1=(1-QPR)/(2-PPR-QPR);   @Pr[S_0=1/Y_0], STEADY STATE PROB.@
         PROB__0=1-PROB__1  ;           @Pr[S_0=0/Y_0], STEADY STATE PROB @

/*  START ITERATION */
/*===================================================================== */

       LIKV=0.0;
       J_ITER = 1;
       DO UNTIL J_ITER>T;

       /* KALMAN FILTER */
       /*==============================================================*/

             B_TL00 = MU_M0+ F * B_LL0;  @WHEN S_{t-1}=0, S_{t}=0@
             B_TL01 = MU_M1+ F * B_LL0;  @WHEN S_{t-1}=0, S_{t}=1@
             B_TL10 = MU_M0+ F * B_LL1;  @WHEN S_{t-1}=1, S_{t}=0@
             B_TL11 = MU_M1+ F * B_LL1;  @WHEN S_{t-1}=1, S_{t}=1@

             P_TL00 = F * P_LL0 * F' + Q;
             P_TL01 = F * P_LL0 * F' + Q;
             P_TL10 = F * P_LL1 * F' + Q;
             P_TL11 = F * P_LL1 * F' + Q;

             F_CAST00= yy[.,J_ITER]- H * B_TL00;
             F_CAST01= yy[.,J_ITER]- H * B_TL01;
             F_CAST10= yy[.,J_ITER]- H * B_TL10;
             F_CAST11= yy[.,J_ITER]- H * B_TL11;

             SS00= H * P_TL00 * H' +R; @VARIANCE OF FORECAST ERROR@
             SS01= H*  P_TL01 * H' +R;
             SS10= H * P_TL10 * H' +R;
             SS11= H * P_TL11 * H' +R;

             B_TT00 = B_TL00 + (P_TL00 * H' *INV(SS00)) * F_CAST00;
             B_TT01 = B_TL01 + (P_TL01 * H' *INV(SS01)) * F_CAST01;
             B_TT10 = B_TL10 + (P_TL10 * H' *INV(SS10)) * F_CAST10;
             B_TT11 = B_TL11 + (P_TL11 * H' *INV(SS11)) * F_CAST11;

             P_TT00 = (EYE(12) - (P_TL00 * H'*INV(SS00)) * H ) * P_TL00;
             P_TT01 = (EYE(12) - (P_TL01 * H'*INV(SS01)) * H ) * P_TL01;
             P_TT10 = (EYE(12) - (P_TL10 * H'*INV(SS10)) * H ) * P_TL10;
             P_TT11 = (EYE(12) - (P_TL11 * H'*INV(SS11)) * H ) * P_TL11;

       /* HAMILTON FILTER */
       /*==============================================================*/

             PR_VL00=V_PROB(F_CAST00,SS00)*  QPR*PROB__0;    @Pr[St,Yt/Yt-1]@
             PR_VL01=V_PROB(F_CAST01,SS01)*(1-QPR)*PROB__0;
             PR_VL10=V_PROB(F_CAST10,SS10)*(1-PPR)*PROB__1;
             PR_VL11=V_PROB(F_CAST11,SS11)*    PPR*PROB__1;
                                          @CONDITIONAL DENSITY TIMES WEIGHT@

             PR_VAL=PR_VL00+PR_VL01+PR_VL10+PR_VL11;
                                  @WEIGHTED AVERAGE OF CONDITIONAL DENSITIES:
                                   f(y_t|Y_{t-1})@

             PRO_00=PR_VL00/PR_VAL;       @Pr[St,St-1/Yt]@
             PRO_01=PR_VL01/PR_VAL;
             PRO_10=PR_VL10/PR_VAL;
             PRO_11=PR_VL11/PR_VAL;

             PROB__0=PRO_00+PRO_10;       @Pr[St=0/Yt]@
             PROB__1=PRO_01+PRO_11;       @Pr[St=1/Yt]@


       /* COLLAPSING TERMS */
       /*==============================================================*/


              B_LL0=(PRO_00*B_TT00 + PRO_10*B_TT10)/PROB__0;
              B_LL1=(PRO_01*B_TT01 + PRO_11*B_TT11)/PROB__1;


              P_LL0=(PRO_00*(P_TT00+(B_LL0-B_TT00)*(B_LL0-B_TT00)')+
                     PRO_10*(P_TT10+(B_LL0-B_TT10)*(B_LL0-B_TT10)'))/PROB__0;

              P_LL1=(PRO_01*(P_TT01+(B_LL1-B_TT01)*(B_LL1-B_TT01)')+
                     PRO_11*(P_TT11+(B_LL1-B_TT11)*(B_LL1-B_TT11)'))/PROB__1;

              IF J_ITER < START; goto skip; endif;
              LIKV = LIKV -LN(PR_VAL);

                skip:
       J_ITER = J_ITER+1;
       ENDO;

/*  END ITERATION    */
/*===================*/

RETP(LIKV);
ENDP;


@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
proc trans(c0);                 @..Constraining Values of Reg. Coff..@

  local c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,cc;

    cc=c0;

    c1=c0[1,1]/(1+abs(c0[1,1]));
    c2=c0[2,1]/(1+abs(c0[2,1]));
    cc[1,1]=c1+c2;
    cc[2,1]=-1*c1*c2;

    c3=c0[3,1]/(1+abs(c0[3,1]));
    c4=c0[4,1]/(1+abs(c0[4,1]));
    cc[3,1]=c3+c4;
    cc[4,1]=-1*c3*c4;

    c5=c0[5,1]/(1+abs(c0[5,1]));
    c6=c0[6,1]/(1+abs(c0[6,1]));
    cc[5,1]=c5+c6;
    cc[6,1]=-1*c5*c6;

    c7=c0[7,1]/(1+abs(c0[7,1]));
    c8=c0[8,1]/(1+abs(c0[8,1]));
    cc[7,1]=c7+c8;
    cc[8,1]=-1*c7*c8;

    c9=c0[9,1]/(1+abs(c0[9,1]));
    c10=c0[10,1]/(1+abs(c0[10,1]));
    cc[9,1]=c9+c10;
    cc[10,1]=-1*c9*c10;

    cc[19:21]=c0[19:21]/10;
    cc[24:25,1]=EXP(C0[24:25,1]) ./ (1 + EXP(c0[24:25,1]));

retp(cc);
endp;


@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
PROC V_PROB(EV, HE);      @ CALCULATES    Pr[Yt/St,Yt-1]@

         LOCAL VAL;
         VAL=(1/SQRT(DET(HE)))*EXP(-0.5*EV'*INV(HE)*EV);

RETP(VAL);
ENDP;

@========================================================================@
@========================================================================@
PROC (3) = FILTER(PRMTR1);

LOCAL PRMTR, PPR,QPR, PHI1,PHI2, VAR_W, MU0,MU1,R,Q,H,F,B_LL0,B_LL1,
      VECP_LL,P_LL0,P_LL1,PROB__1,PROB__0,LIKV,J_ITER,B_TL00,B_TL01,
      B_TL10,B_TL11,P_TL00,P_TL01,P_TL10,P_TL11,F_CAST00,F_CAST01,
      F_CAST10,F_CAST11,SS00,SS01,SS10,SS11,B_TT00,B_TT01,B_TT10,
      B_TT11,P_TT00,P_TT01,P_TT10,P_TT11,PR_VL00,PR_VL01,PR_VL10,PR_VL11,
      PR_VAL,PRO_00,PRO_01,PRO_10,PRO_11,MU_M0,MU_M1,DCCI_M,
      PSI11,PSI12,PSI21,PSI22,PSI31,PSI32,PSI41,PSI42,PR_TT0M,PR_TL0M,
      SIG1,SIG2,SIG3,SIG4,GAMMA1,GAMMA2,GAMMA3,GAMMA4,SIG_E,B_TT,
      gamma41,gamma42,gamma43,K_GN00,K_GN01,K_GN10,K_GN11,SS_K_GN;

           PRMTR=TRANS(PRMTR1);

           LOCATE 20,1;  PRMTR';

    phi1=prmtr[1,1];
    phi2=prmtr[2,1];
    psi11=prmtr[3,1];
    psi12=prmtr[4,1];
    psi21=prmtr[5,1];
    psi22=prmtr[6,1];
    psi31=prmtr[7,1];
    psi32=prmtr[8,1];
    psi41=prmtr[9,1];
    psi42=prmtr[10,1];

    sig1=prmtr[11,1];
    sig2=prmtr[12,1];
    sig3=prmtr[13,1];
    sig4=prmtr[14,1];

    gamma1=prmtr[15,1];
    gamma2=prmtr[16,1];
    gamma3=prmtr[17,1];
    gamma4=prmtr[18,1];

    gamma41=prmtr[19,1];
    gamma42=prmtr[20,1];
    gamma43=prmtr[21,1];

    MU0=PRMTR[22,1];
    MU1=PRMTR[23,1];

    QPR=PRMTR[24,1];    @Pr[St=0/St-1=0]@
    PPR=PRMTR[25,1];    @Pr[St=1/St-1=1]@

@=====================================================================@

 F=(phi1~phi2~  0~  0~    0~   0~   0~   0~  0~  0~  0~  0)|
   (   1~   0~  0~  0~    0~   0~   0~   0~  0~  0~  0~  0)|
   (   0~   1~  0~  0~    0~   0~   0~   0~  0~  0~  0~  0)|
   (   0~   0~  1~  0~    0~   0~   0~   0~  0~  0~  0~  0)|
   (   0~   0~  0~  0~psi11~psi12~  0~   0~  0~  0~  0~  0)|
   (   0~   0~  0~  0~    1~   0~   0~   0~  0~  0~  0~  0)|
   (   0~   0~  0~  0~    0~   0~psi21~psi22~  0~  0~  0~  0)|
   (   0~   0~  0~  0~    0~   0~   1~   0~  0~  0~  0~  0)|
   (   0~   0~  0~  0~    0~   0~   0~   0~psi31~psi32~  0~  0)|
   (   0~   0~  0~  0~    0~   0~   0~   0~  1~  0~  0~  0)|
   (   0~   0~  0~  0~    0~   0~   0~   0~  0~  0~psi41~psi42)|
   (   0~   0~  0~  0~    0~   0~   0~   0~  0~  0~  1~  0);

   H=(gamma1~0~0~0~                  1~0~0~0~0~0~0~0)|
     (gamma2~0~0~0~                  0~0~1~0~0~0~0~0)|
     (gamma3~0~0~0~                  0~0~0~0~1~0~0~0)|
     (gamma4~gamma41~gamma42~gamma43~0~0~0~0~0~0~1~0);

    sig_e=1;

    Q=(sig_e^2~0~0~0~    0~0~    0~0~    0~0~    0~0)|
      (   0~   0~0~0~    0~0~    0~0~    0~0~    0~0)|
      (   0~   0~0~0~    0~0~    0~0~    0~0~    0~0)|
      (   0~   0~0~0~    0~0~    0~0~    0~0~    0~0)|
      (   0~0~0~0~  sig1^2~0~    0~0~    0~0~    0~0)|
      (   0~0~0~0~       0~0~    0~0~    0~0~    0~0)|
      (   0~0~0~0~       0~0~sig2^2~0~   0~0~    0~0)|
      (   0~0~0~0~       0~0~    0~0~    0~0~    0~0)|
      (   0~0~0~0~       0~0~    0~0~sig3^2~0~   0~0)|
      (   0~0~0~0~       0~0~    0~0~    0~0~    0~0)|
      (   0~0~0~0~       0~0~    0~0~    0~0~sig4^2~0)|
      (   0~0~0~0~       0~0~    0~0~    0~0~    0~0);

    R=0;

     MU_M0=MU0|ZEROS(11,1);
     MU_M1=MU1|ZEROS(11,1);


DCCI_M=ZEROS(T,1);  @will store FIRST DIFF. OF DEMEANED NEW INDEX:
                     (1,1) element of B_TT@
PR_TT0M=ZEROS(T,1);  @will store Pr[S_t=0|Y_t]@
PR_TL0M=ZEROS(T,1);  @will store Pr[S_t=0|Y_{t-1}]@

/*  NOTE:  KIM'S FILTER = HAMILTON FILTER + KALMAN FILTER               */
/*  INITIALIZING THE FILTER: FOR BOTH KALMAN FILTER AND HAMILTON FILTER */

         B_LL0= INV(EYE(12)-F)* MU_M0;
         B_LL1= INV(EYE(12)-F)* MU_M1;

         VECP_LL=INV(EYE(144) - F.*.F)*vec(Q);

    P_LL0=vecP_LL[1:12,1]~vecP_LL[13:24,1]~vecP_LL[25:36,1]~vecP_LL[37:48,1]
        ~vecP_LL[49:60,1]~vecP_LL[61:72,1]~vecP_LL[73:84,1]~vecP_LL[85:96,1]
        ~vecP_LL[97:108,1]~vecP_LL[109:120,1]~vecP_LL[121:132,1]~
         vecP_LL[133:144,1];

    P_LL1=P_LL0;

         PROB__1=(1-QPR)/(2-PPR-QPR);   @Pr[S_0=1/Y_0], STEADY STATE PROB.@
         PROB__0=1-PROB__1  ;           @Pr[S_0=0/Y_0], STEADY STATE PROB @

/*  START ITERATION */
/*===================================================================== */

       LIKV=0.0;
       J_ITER = 1;
       DO UNTIL J_ITER>T;

PR_TL0M[J_ITER,1]= QPR*PROB__0 + (1-PPR)*PROB__1;     @Pr[St=0/Yt-1]@

       /* KALMAN FILTER */
       /*==============================================================*/

             B_TL00 = MU_M0+ F * B_LL0;  @WHEN S_{t-1}=0, S_{t}=0@
             B_TL01 = MU_M1+ F * B_LL0;  @WHEN S_{t-1}=0, S_{t}=1@
             B_TL10 = MU_M0+ F * B_LL1;  @WHEN S_{t-1}=1, S_{t}=0@
             B_TL11 = MU_M1+ F * B_LL1;  @WHEN S_{t-1}=1, S_{t}=1@

             P_TL00 = F * P_LL0 * F' + Q;
             P_TL01 = F * P_LL0 * F' + Q;
             P_TL10 = F * P_LL1 * F' + Q;
             P_TL11 = F * P_LL1 * F' + Q;

             F_CAST00= yy[.,J_ITER]- H * B_TL00;
             F_CAST01= yy[.,J_ITER]- H * B_TL01;
             F_CAST10= yy[.,J_ITER]- H * B_TL10;
             F_CAST11= yy[.,J_ITER]- H * B_TL11;

             SS00= H * P_TL00 * H' +R; @VARIANCE OF FORECAST ERROR@
             SS01= H*  P_TL01 * H' +R;
             SS10= H * P_TL10 * H' +R;
             SS11= H * P_TL11 * H' +R;

             B_TT00 = B_TL00 + (P_TL00 * H' *INV(SS00)) * F_CAST00;
             B_TT01 = B_TL01 + (P_TL01 * H' *INV(SS01)) * F_CAST01;
             B_TT10 = B_TL10 + (P_TL10 * H' *INV(SS10)) * F_CAST10;
             B_TT11 = B_TL11 + (P_TL11 * H' *INV(SS11)) * F_CAST11;

             P_TT00 = (EYE(12) - (P_TL00 * H'*INV(SS00)) * H ) * P_TL00;
             P_TT01 = (EYE(12) - (P_TL01 * H'*INV(SS01)) * H ) * P_TL01;
             P_TT10 = (EYE(12) - (P_TL10 * H'*INV(SS10)) * H ) * P_TL10;
             P_TT11 = (EYE(12) - (P_TL11 * H'*INV(SS11)) * H ) * P_TL11;

       /* HAMILTON FILTER */
       /*==============================================================*/


             PR_VL00=V_PROB(F_CAST00,SS00)*  QPR*PROB__0;    @Pr[St,Yt/Yt-1]@
             PR_VL01=V_PROB(F_CAST01,SS01)*(1-QPR)*PROB__0;
             PR_VL10=V_PROB(F_CAST10,SS10)*(1-PPR)*PROB__1;
             PR_VL11=V_PROB(F_CAST11,SS11)*    PPR*PROB__1;
                                          @CONDITIONAL DENSITY TIMES WEIGHT@

             PR_VAL=PR_VL00+PR_VL01+PR_VL10+PR_VL11;
                                  @WEIGHTED AVERAGE OF CONDITIONAL DENSITIES:
                                   f(y_t|Y_{t-1})@

             PRO_00=PR_VL00/PR_VAL;       @Pr[St,St-1/Yt]@
             PRO_01=PR_VL01/PR_VAL;
             PRO_10=PR_VL10/PR_VAL;
             PRO_11=PR_VL11/PR_VAL;

B_TT=B_TT00*PRO_00+B_TT01*PRO_01+B_TT10*PRO_10+B_TT11*PRO_11;
DCCI_M[J_ITER,1]=B_TT[1,1];

             PROB__0=PRO_00+PRO_10;       @Pr[St=0/Yt]@
             PROB__1=PRO_01+PRO_11;       @Pr[St=1/Yt]@

PR_TT0M[J_ITER,1]=PROB__0;

       /* COLLAPSING TERMS */
       /*==============================================================*/


              B_LL0=(PRO_00*B_TT00 + PRO_10*B_TT10)/PROB__0;
              B_LL1=(PRO_01*B_TT01 + PRO_11*B_TT11)/PROB__1;


              P_LL0=(PRO_00*(P_TT00+(B_LL0-B_TT00)*(B_LL0-B_TT00)')+
                     PRO_10*(P_TT10+(B_LL0-B_TT10)*(B_LL0-B_TT10)'))/PROB__0;

              P_LL1=(PRO_01*(P_TT01+(B_LL1-B_TT01)*(B_LL1-B_TT01)')+
                     PRO_11*(P_TT11+(B_LL1-B_TT11)*(B_LL1-B_TT11)'))/PROB__1;


              IF J_ITER < START; goto skip; endif;
              LIKV = LIKV -LN(PR_VAL);

                skip:
       J_ITER = J_ITER+1;
       ENDO;

/*  END ITERATION    */
/*===================*/



RETP(DCCI_M,PR_TT0M,PR_TL0M);
ENDP;

@==========================================================@

PROC SMOOTH(pr_tt0,pr_tl0);

         @pr_TT0 contains Pr[S_t|Y_t]@
         @pr_TL0 contains Pr[S_t|Y_{t-1}]@

local ppr, qpr, pr_sm0,pr_sm1, j_iter,pr_sm00,pr_sm01,pr_sm10,pr_sm11,
      pr_tt1,pr_tl1,TT,prmtr;

          TT=ROWS(PR_TT0);
          prmtr=trans(xout);

           PPR=PRMtr[25,1];    @Pr[St=1/St-1=1]@
           QPR=PRMtr[24,1];    @Pr[St=0/St-1=0]@

           pr_tt1=1-pr_tt0;
           pr_tl1=1-pr_tl0;

           pr_sm0=pr_tt0;     @ pr_sm0 will contain Pr[S_t|Y_T]@
           pr_sm1=pr_tt1;

j_iter=TT-1;
do until j_iter < 1;

    @The followings are P[S_t, S_t+1|Y_T] @

   pr_sm00=pr_sm0[j_iter+1,1]*qpr*    pr_tt0[j_iter,1]/ pr_tl0[j_iter+1,1];

   pr_sm01=pr_sm1[j_iter+1,1]*(1-qpr)*pr_tt0[j_iter,1]/ pr_tl1[j_iter+1,1];

   pr_sm10=pr_sm0[j_iter+1,1]*(1-ppr)*pr_tt1[j_iter,1]/ pr_tl0[j_iter+1,1];

   pr_sm11=pr_sm1[j_iter+1,1]*ppr*    pr_tt1[j_iter,1]/ pr_tl1[j_iter+1,1];

   pr_sm0[j_iter,1]=pr_sm00+pr_sm01;
   pr_sm1[j_iter,1]=pr_sm10+pr_sm11;

j_iter=j_iter -1;
endo;

RETP(pr_sm0); @This proc returns Pr[S_t=0|Y_T]@
endp;
@++++++++++++++++++++++++++++++++++++++++++++++++++++++++@