new;
library  optmum, pgraph;

format /m1 /rd 9,6;
load yy[167,1]=a:\lpgnp.txt; 
yy=100*yy;

t=rows(yy);

output file=a:\uc.out reset;
output off;

@ Maximum Likelihood Estimation @

START=2;        @ startup values for VEVD of likelihood @

        @Initial values of the parameters @

PRMTR_IN={   0.5	0.5	0	0	0 };

PRMTR_IN=PRMTR_IN';

@ maximize likelihood function @

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

@ prmtr estimates, -log lik value, Grandient, code@

prm_fnl=trans(xout);


"Calculating Hessian..... Please be patient!!!!";

                  hessn0=hessp(&lik_fcn,xout);

                  cov0=inv(hessn0);



grdn_fnl=gradfd(&TRANS,xout);

cov=grdn_fnl*cov0*grdn_fnl';

SD =sqrt(diag(cov)); @Standard errors of the estimated coefficients@


output on;

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

"initial values of prmtr is";

                 trans(prmtr_in)';

"==============================================================";

"likelihood value is ";; -fout;

"code";;cout;

"Estimated parameters are:";

prm_fnl';

XOUT';

"Standard errors of parameters are:"; sd';

OUTPUT OFF;

f=prm_fnl[4]~prm_fnl[5]|1~0;
ef=eig(f);
output on; ef; output off;


DATA = filter(xout);



output file=a:\uccycle.dta reset;

DATA[.,2]/100;

output off;


/*
XY(SEQA(1,1,ROWS(DATA)),YY[START:T,1]~DATA[.,1]);

XY(SEQA(1,1,ROWS(DATA)),DATA[.,2]); */


@ END OF MAIN PROGRAM @

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

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


PROC LIK_FCN(PRMTR1);

     LOCAL SN,SE,SNE,F,F1,H,Q,Q1,R,BETA_LL,P_LL,BETA_TL,P_TL,BETA_TT,P_TT,
     P_LL1, vt,ft, VAL,LIK_MAT,J_ITER, phi1,phi2,vecpll,prmtr,muvec,mu;

          @ transform hyperparameters to impose contraints @


     PRMTR=TRANS(PRMTR1);

@     PRMTR = PRMTR1; @

LOCATE 20,1;  PRMTR';

     SN=PRMTR[1];  @s.e. of the random walk component@
     SE=PRMTR[2];  @s.e. of the AR component@
     mu=PRMTR[3];
     phi1=PRMTR[4];
     phi2=PRMTR[5];

     muvec = mu|0|0;            @ drift vector @
     F=(1~0~0)|                @ transition matrix @
       (0~phi1~phi2)|
       (0~1~0);
     F1 = submat(F,2~3,2~3);   @ pick out I(0) part @
     H=1~1~0;                  @ measurement equation @


     Q= ZEROS(3,3);             @ E[V(t)V(t)'] @

     Q[1,1]=SN^2;
     Q[2,2]=SE^2;


     Q1 = submat(Q,2~3,2~3);    @ cov matrix of I(0) part @
     R= 0;                      @ cov matrix of error in measurement eqn @

     BETA_LL=yy[1]|0|0;             @ starting values @



     @ variance matrix of initial state vector @


          P_LL = zeros(3,3);
          P_LL1 = inv( (eye(4)-F1.*.F1) )*vec(Q1);
          P_LL1 = reshape(P_LL1',2,2);
          P_LL[1,1]=100000000;
          P_LL[2:3,2:3]=P_LL1;

     LIK_MAT=ZEROS(T,1);

J_ITER = 1;

DO UNTIL J_ITER>T;

      BETA_TL = muvec + F*BETA_LL ;
      P_TL = F*P_LL*F' + Q;

      vt=yy[j_iter,1] - H * BETA_TL ; @prediction error@

      ft= H * P_TL * H' + R;       @variance of forecast error@

      BETA_TT= BETA_TL + P_TL * H' * inv(ft) * vt;

      P_TT= P_TL - P_TL * H' * inv(ft) * H * P_TL;

      LIK_MAT[J_ITER,1] =  0.5*(LN(2*pi*ft) + vt^2/ft);

      BETA_LL=BETA_TT;

      P_LL=P_TT;

J_ITER = J_ITER+1;

ENDO;

VAL = SUMC(LIK_MAT[START:T]);

LOCATE 2,20;  VAL;

RETP(VAL);

ENDP;


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



PROC FILTER(PRMTR1);



LOCAL SN,SE,SNE,F,F1,H,Q,Q1,R,BETA_LL,P_LL,BETA_TL,P_TL,BETA_TT,P_TT,
P_LL1, vt,ft, VAL,LIK_MAT,J_ITER, phi1,phi2,vecpll,prmtr,muvec,mu,
BETA_MAT;


     BETA_MAT=ZEROS(T,2);



     LIK_MAT=ZEROS(T,1);



     PRMTR=TRANS(PRMTR1);



LOCATE 20,1;  PRMTR';

     SN=PRMTR[1];  @s.e. of the random walk component@
     SE=PRMTR[2];  @s.e. of the AR component@
     mu=PRMTR[3];
     phi1=PRMTR[4];
     phi2=PRMTR[5];

     muvec = mu|0|0;            @ drift vector @
     F=(1~0~0)|                @ transition matrix @
       (0~phi1~phi2)|
       (0~1~0);
     F1 = submat(F,2~3,2~3);   @ pick out I(0) part @
     H=1~1~0;                  @ measurement equation @


     Q= ZEROS(3,3);             @ E[V(t)V(t)'] @

     Q[1,1]=SN^2;
     Q[2,2]=SE^2;

     Q1 = submat(Q,2~3,2~3);    @ cov matrix of I(0) part @
     R= 0;                      @ cov matrix of error in measurement eqn @

     BETA_LL=yy[1]|0|0;             @ starting values @



     @ variance matrix of initial state vector @


          P_LL = zeros(3,3);
          P_LL1 = inv( (eye(4)-F1.*.F1) )*vec(Q1);
          P_LL1 = reshape(P_LL1',2,2);
          P_LL[1,1]=100000000;
          P_LL[2:3,2:3]=P_LL1;


J_ITER = 1;

DO UNTIL J_ITER>T;


      BETA_TL = muvec + F * BETA_LL ;

      P_TL = F * P_LL * F' + Q;

      vt=yy[j_iter,1] - H * BETA_TL ;

                  @prediction error@


      ft= H * P_TL * H' + R;       @variance of forecast error@

      BETA_TT= BETA_TL + P_TL * H' * inv(ft) * vt;

      P_TT= P_TL - P_TL * H' * inv(ft) * H * P_TL;


BETA_MAT[J_ITER,.]=BETA_TT[1,1]~BETA_TT[2,1];

      BETA_LL=BETA_TT;

      P_LL=P_TT;


J_ITER = J_ITER+1;

ENDO;


RETP(BETA_MAT[START:T,.]);

ENDP;



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

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

local c1,aaa,ccc;

      c1 = c0;


c1[1:2]=exp(-c0[1:2]);

        aaa=c0[4]./(1+abs(c0[4]));
     	ccc=(1-abs(aaa))*c0[5]./(1+abs(c0[5]))+abs(aaa)-aaa^2;

     c1[4]=2*aaa;
     c1[5]= -1* (aaa^2+ccc) ; 


retp(c1);

endp;