/****************
This program estimates infrequent monetary policy shocks from a a bivariate model of 
SP500 and federal funds rate shocks (based on Wu-Xia 2016). The VARs are 5 variable, 5 lags (AIC).
The data is monthly. 738 datapoints, 1955m07 to 2016m12.

It decomposes the reduced form VAR residuals into three orthogonal parts. All three
parts are serially uncorrelated. 

   e1(t) = u1(t) + a12*S(t)*u3(t),
   e2(t) = u2(t) + S(t)*u3(t),

S(t) = 0, 1.

P(S(t)=0|S(t-1)=0) = p
P(S(t)=1|S(t-1)=1) = q. 


u3(t) = monetary policy shock. We have 6 parameters.

Only u3(t) is estimated as a state variable.

Computes the filtered (one-sided) shocks.

The program will work for Table 2 full sample, our core results. After restrictions, it can be used for 
Table 2 1955:07 - 2008:12 results. Also this cane be used, after appropriate sample or parameters restrictions, 
for Tables 5 and 6.

We used Gauss 9. Optmum optimizer with BFGS.

                                ************/


New;
library optmum;
#include optmum.ext
#include gradient.ext;
optset;


@ These are the reduced form VAR residuals for stock returns and Federal funds rate @

load data[738,2]=C:\arabinda\research\ds2014\rcdata04.txt; 

yy =data[1:738,.];


t=rows(yy);
format /m1 /rd 8,4;

output file=C:\arabinda\research\xxx\xxx.out reset;
output off;
@================= Initialize Global Variables============@

           START=1;  
           
           PRMTR_NO=6;   @ Number of parameters to be estimated @
           dmnsion=1;     @ dimension of transition equation @

pr_1=-2.8;
pr_2=-1.7;pr_3=-0.0;pr_4=-1.7;

pr_5=-2.1;

pr_6=-2.2;


@initial values of parameters@

PRMTR_IN=pr_1|pr_2|pr_3|pr_4|pr_5|pr_6;



@ Maximum Likelihood Estimation @

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

@ prmtr estimates, -log lik value, Gradient, return code @

  hout=inv(hessp(&lik_fcn,xout));


PRM_FNL=trans(xout);   @ Estimated  coefficients, constrained@

grdn_fnl=gradfd(&TRANS,xout);
Hsn_fnl=grdn_fnl*hout*grdn_fnl';
SD_fnl =sqrt(diag(Hsn_fnl)); @Std errors est. coeffs@


output on;
"==FINAL OUTPUT=============================================";
"The code is ";;-cout;
"The Likelihood value is ";; -fout;

"The Estimated parameters are:";
prm_fnl';
"Estimated parameters without transformation are:";
xout';
"Standard errors of parameters are:"; sd_fnl';


"===============================================================";
output off;

{eq_r,Prob_T0} =FILTER(XOUT);

output file=C:\arabinda\research\ds2014\prg021a.out reset;
Eq_R~PROB_T0;


@+++++++++++++++++++++++++++++++++++++++++++++++++++++++@
PROC LIK_FCN(PRMTR1);
Local  prmtr, PR_VAL,H, YT,ZT,LIKV, LIK, J_ITER,lik_m,

b11,b21,
B1,B2,

PR_TR,pr_tr1,
PR_TRF0,pr_trf1,
P, Q, F,qq,vecp,RR,

a12,sig_u1,sig_u2,sig_u3,

H0,H1,
Prob__0, Prob__1,
P_ll0, P_ll1,
b_ll0, b_ll1,

B_TL00, B_TL01,
B_TL10, B_TL11,

p_TL00, p_TL01,
p_TL10, p_TL11,

B_Tt00, B_Tt01,
B_Tt10, B_Tt11,

p_Tt00, p_Tt01,
p_Tt10, p_Tt11,

ss00,ss01,
ss10,ss11,

F_cast00,F_cast01,
F_cast10,F_cast11,

PRO_00, PRO_01,
PRO_10, PRO_11,

PR_VL00, PR_VL01,
PR_VL10, PR_VL11;

EXTERNAL PROC TRANS, V_PROB;
PRMTR=TRANS(PRMTR1);

        p=PRMTR[1,1];   @Pr[St=0|St-1=0}@
        SIG_u1=PRMTR[2,1];
        SIG_u2=PRMTR[3,1];
        SIG_u3=prmtr[4,1];

        a12=prmtr[5,1];
	
	q=PRMTR[6,1];   @Pr[St=1|St-1=1}@


@A matrix of Tr. prob@
PR_TR=  (p~1-q)|(1-p~q);

    @ Initial Probabilities @
           PROB__0=(1-q)/(2-p-q);
           PROB__1=1-prob__0;

           PR_TRF0=VEC(PR_TR);
           
        H0 = (0)|(0);
	    H1 = (a12)|(1);

        F = (0);

	RR= (sig_u1^2~0)|(0~sig_u2^2);

        QQ = (sig_u3^2);
        
        
    B_LL0 = (0);     @initial guess for state vector@
    B_LL1 = (0);


    vecp = inv(eye(1))*QQ;

    
    P_LL0 = vecp;
    P_ll1 = p_ll0;

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

LIK_m = zeros(t,1);
J_ITER = 1;
DO UNTIL J_ITER>T;

YT=YY[J_ITER,.];


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

             B_TL00 =  F * B_LL0;  @WHEN S_{t-1}=0, S_{t}=0@
             B_TL01 =  F * B_LL0;  @WHEN S_{t-1}=0, S_{t}=1@

             B_TL10 =  F * B_LL1;  @WHEN S_{t-1}=1, S_{t}=0@
             B_TL11 =  F * B_LL1;  @WHEN S_{t-1}=1, S_{t}=1@

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


             P_TL00 = F * P_LL0 * F' + QQ;
             P_TL01 = F * P_LL0 * F' + QQ;

             P_TL10 = F * P_LL1 * F' + QQ;
             P_TL11 = F * P_LL1 * F' + QQ;


         @======= Prediction Errors ==@

             F_CAST00= YT' - H0*B_Tl00;
             F_CAST01= YT' - H1*B_Tl01;

             F_CAST10= YT' - H0*B_Tl10;
             F_CAST11= YT' - H1*B_Tl11;

    @++++++++++++++++@

             SS00= H0 * P_TL00 * H0' + RR ; @VARIANCE OF FORECAST ERROR@
             SS01= H1 * P_TL01 * H1' + RR;


             SS10= H0 * P_TL10 * H0' + RR; @VARIANCE OF FORECAST ERROR@
             SS11= H1 * P_TL11 * H1' + RR; 

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


             B_TT00 = B_TL00 + (P_TL00 * H0' *INV(SS00)) * F_CAST00;
             B_TT01 = B_TL01 + (P_TL01 * H1' *INV(SS01)) * F_CAST01;

             B_TT10 = B_TL10 + (P_TL10 * H0' *INV(SS10)) * F_CAST10;
             B_TT11 = B_TL11 + (P_TL11 * H1' *INV(SS11)) * F_CAST11;

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

             P_TT00 = (EYE(1) - P_TL00 * H0'*INV(SS00) * H0 ) * P_TL00;
             P_TT01 = (EYE(1) - P_TL01 * H1'*INV(SS01) * H1 ) * P_TL01;

             P_TT10 = (EYE(1) - P_TL10 * H0'*INV(SS10) * H0 ) * P_TL10;
             P_TT11 = (EYE(1) - P_TL11 * H1'*INV(SS11) * H1 ) * P_TL11;

    @+++++++++++++++++++@
                /* HAMILTON FILTER */
       /*==================================================*/
          @Pr[St,Yt/Yt-1]@

          PR_VL00=V_PROB(F_CAST00,SS00)*PR_TRF0[1]*PROB__0;
          PR_VL01=V_PROB(F_CAST01,SS01)*PR_TRF0[2]*PROB__0;

          PR_VL10=V_PROB(F_CAST10,SS10)*PR_TRF0[3]*PROB__1;
          PR_VL11=V_PROB(F_CAST11,SS11)*PR_TRF0[4]*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})@

lik = -1*ln(pr_val);


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

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

             PRO_10=PR_VL10/PR_VAL;       @Pr[St,St-1/Yt]@
             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]@

             @Pr[St=i/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;


LIK_M[j_iter,.] = LIK; 

J_ITER = J_ITER+1;
ENDO;



likv=sumc(lik_m);
              
RETP(LIKV);
ENDP;

@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
PROC TRANS(c0); @ constraining values of reg. coeff.@
     local cc;

cc = c0;

    cc[1,.]= exp(-1*c0[1,1])./ (1 + exp(-1*c0[1,1]));

    cc[2:4,.] = exp(-1*c0[2:4,.]);

    cc[6,.]= exp(-1*c0[6,1])./ (1 + exp(-1*c0[6,1]));

    
   retp(cc);
endp;


@++++++++++++++++++++++++++++++++@

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

VAL=(1/SQRT(((2*pi)^2)*DET(HE)))*EXP(-0.5*EV'*inv(HE)*EV);


RETP(VAL);
ENDP;

@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@

proc(2)= filter(prmtr1);
Local  prmtr, PR_VAL,H, YT,ZT,LIKV, LIK, J_ITER,lik_m,
r_eq,pr_t0,

b11,b21,
B1,B2,

PR_TR,pr_tr1,
PR_TRF0,pr_trf1,
P, Q, F,qq,vecp,RR,

a12,sig_u1,sig_u2,sig_u3,

H0,H1,
Prob__0, Prob__1,
P_ll0, P_ll1,
b_ll0, b_ll1,

B_TL00, B_TL01,
B_TL10, B_TL11,

p_TL00, p_TL01,
p_TL10, p_TL11,

B_Tt00, B_Tt01,
B_Tt10, B_Tt11,

p_Tt00, p_Tt01,
p_Tt10, p_Tt11,

ss00,ss01,
ss10,ss11,

F_cast00,F_cast01,
F_cast10,F_cast11,

PRO_00, PRO_01,
PRO_10, PRO_11,

PR_VL00, PR_VL01,
PR_VL10, PR_VL11;

EXTERNAL PROC TRANS, V_PROB;
PRMTR=TRANS(PRMTR1);



        p=PRMTR[1,1];   @Pr[St=0|St-1=0}@
        SIG_u1=PRMTR[2,1];
        SIG_u2=PRMTR[3,1];
        SIG_u3=prmtr[4,1];

        a12=prmtr[5,1];
	
	q=PRMTR[6,1];   @Pr[St=1|St-1=1}@


@A matrix of Tr. prob@
PR_TR=  (p~1-q)|(1-p~q);

    @ Initial Probabilities @
           PROB__0=(1-q)/(2-p-q);
           PROB__1=1-prob__0;

           PR_TRF0=VEC(PR_TR);
           
        H0 = (0)|(0);
	    H1 = (a12)|(1);

        F = (0);

	RR= (sig_u1^2~0)|(0~sig_u2^2);

        QQ = (sig_u3^2);
        
        
    B_LL0 = (0);     @initial guess for state vector@
    B_LL1 = (0);


    vecp = inv(eye(1))*QQ;

    
    P_LL0 = vecp;
    P_ll1 = p_ll0;


r_eq=zeros(t,1);

PR_T0=zeros(t,1);

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

LIK_m = zeros(t,1);
J_ITER = 1;
DO UNTIL J_ITER>T;

YT=YY[J_ITER,.];

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

             B_TL00 =  F * B_LL0;  @WHEN S_{t-1}=0, S_{t}=0@
             B_TL01 =  F * B_LL0;  @WHEN S_{t-1}=0, S_{t}=1@

             B_TL10 =  F * B_LL1;  @WHEN S_{t-1}=1, S_{t}=0@
             B_TL11 =  F * B_LL1;  @WHEN S_{t-1}=1, S_{t}=1@

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


             P_TL00 = F * P_LL0 * F' + QQ;
             P_TL01 = F * P_LL0 * F' + QQ;

             P_TL10 = F * P_LL1 * F' + QQ;
             P_TL11 = F * P_LL1 * F' + QQ;


         @======= Prediction Errors ==@

             F_CAST00= YT' - H0*B_Tl00;
             F_CAST01= YT' - H1*B_Tl01;

             F_CAST10= YT' - H0*B_Tl10;
             F_CAST11= YT' - H1*B_Tl11;

    @++++++++++++++++@

             SS00= H0 * P_TL00 * H0' + RR ; @VARIANCE OF FORECAST ERROR@
             SS01= H1 * P_TL01 * H1' + RR;


             SS10= H0 * P_TL10 * H0' + RR; @VARIANCE OF FORECAST ERROR@
             SS11= H1 * P_TL11 * H1' + RR; 

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


             B_TT00 = B_TL00 + (P_TL00 * H0' *INV(SS00)) * F_CAST00;
             B_TT01 = B_TL01 + (P_TL01 * H1' *INV(SS01)) * F_CAST01;

             B_TT10 = B_TL10 + (P_TL10 * H0' *INV(SS10)) * F_CAST10;
             B_TT11 = B_TL11 + (P_TL11 * H1' *INV(SS11)) * F_CAST11;

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

             P_TT00 = (EYE(1) - P_TL00 * H0'*INV(SS00) * H0 ) * P_TL00;
             P_TT01 = (EYE(1) - P_TL01 * H1'*INV(SS01) * H1 ) * P_TL01;

             P_TT10 = (EYE(1) - P_TL10 * H0'*INV(SS10) * H0 ) * P_TL10;
             P_TT11 = (EYE(1) - P_TL11 * H1'*INV(SS11) * H1 ) * P_TL11;


    @+++++++++++++++++++@
                 /* HAMILTON FILTER */
       /*==================================================*/
          @Pr[St,Yt/Yt-1]@

          PR_VL00=V_PROB(F_CAST00,SS00)*PR_TRF0[1]*PROB__0;
          PR_VL01=V_PROB(F_CAST01,SS01)*PR_TRF0[2]*PROB__0;

          PR_VL10=V_PROB(F_CAST10,SS10)*PR_TRF0[3]*PROB__1;
          PR_VL11=V_PROB(F_CAST11,SS11)*PR_TRF0[4]*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})@

lik = -1*ln(pr_val);


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

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

             PRO_10=PR_VL10/PR_VAL;       @Pr[St,St-1/Yt]@
             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]@

             @Pr[St=i/Yt]@


r_eq[j_iter,.]=PRO_00*B_TT00[1,1] + PRO_10*B_TT10[1,1]
        +PRO_01*B_TT01[1,1] + PRO_11*B_TT11[1,1];

PR_T0[j_iter,.]=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;



LIK_M[j_iter,.] = LIK;
J_ITER = J_ITER+1;
ENDO;



retp(r_eq,pr_t0);
endp;


@~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~@

end;

