/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

arima212.opt

This program generates Table 2 in Morley, Nelson, and Zivot, 
"Why Are Beveridge-Nelson and Unobserved-Component Decompositions 
of GDP So Different?" (forthcoming, RESTAT).

If you use this program, please cite the paper.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

new;
library  optmum, pgraph;

format /m1 /rd 9,6;

            @Real GDP 1947.1 - 1998.2; From CITIBASE@

load data[206,1]=lngdpq.txt;  
        @Real GDP 1947.1 - 1998.2; From CITIBASE@



yy= 100*data;
dyy=yy[2:206]-yy[1:205];
t=rows(dyy);

output file=arima212.out reset;
output off;

@ Maximum Likelihood Estimation @

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

        @Initial values of the parameters @


PRMTR_IN={  0.012160  1.5 -0.6  0.809842    0   0 };
PRMTR_IN={  0   0   0   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;

DATA = filter(prm_fnl);



output file=arima212.dta reset;

dyy~DATA;

output off;







@ 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,G,
     P_LL1, vt,ft, VAL,LIK_MAT,J_ITER, phi1,phi2,vecpll,prmtr,muvec,mu,theta1,
theta2;

          @ transform hyperparameters to impose contraints @


     PRMTR=TRANS(PRMTR1);

@     PRMTR = PRMTR1; @

LOCATE 20,1;  PRMTR';

     SE=PRMTR[1];  @s.e. of the random walk component@
     phi1=PRMTR[2];  @s.e. of the AR component@
     phi2=PRMTR[3]; @ COV(n,e) @
     mu=PRMTR[4];
      theta1=PRMTR[5];
     theta2=PRMTR[6];


     F=(phi1~phi2~theta1~theta2)|                @ transition matrix @
       (1~0~0~0)|
       (0~0~0~0)|
       (0~0~1~0);

     H=1~0~0~0;                  @ measurement equation @


     Q=SE^2;             @ E[V(t)V(t)'] @
    G=1|0|1|0;


     R= 0;                      @ cov matrix of error in measurement eqn @

    BETA_LL=0|0|0|0;        @ starting values @



     @ variance matrix of initial state vector @

          P_LL = inv( (eye(16)-F.*.F) )*vec(G*Q*G');
          P_LL=  reshape(P_LL,4,4);



     LIK_MAT=ZEROS(T,1);

J_ITER = 1;

DO UNTIL J_ITER>T;

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

      vt=dyy[j_iter,1] - mu - 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,G,
P_LL1, vt,ft, VAL,LIK_MAT,J_ITER, phi1,phi2,vecpll,prmtr,muvec,mu,theta1,theta2,
BETA_MAT;


     BETA_MAT=ZEROS(T,1);



     LIK_MAT=ZEROS(T,1);



     PRMTR=PRMTR1;



LOCATE 20,1;  PRMTR';

     SE=PRMTR[1];  @s.e. of the random walk component@
     phi1=PRMTR[2];  @s.e. of the AR component@
     phi2=PRMTR[3]; @ COV(n,e) @
     mu=PRMTR[4];
     theta1=PRMTR[5];
     theta2=PRMTR[6];


     F=(phi1~phi2~theta1~theta2)|                @ transition matrix @
       (1~0~0~0)|
       (0~0~0~0)|
       (0~0~1~0);

     H=1~0~0~0;                  @ measurement equation @


     Q=SE^2;             @ E[V(t)V(t)'] @
    G=1|0|1|0;


     R= 0;                      @ cov matrix of error in measurement eqn @

    BETA_LL=0|0|0|0;         @ starting values @



     @ variance matrix of initial state vector @

          P_LL = inv( (eye(16)-F.*.F) )*vec(G*Q*G');
          P_LL=  reshape(P_LL,4,4);




J_ITER = 1;

DO UNTIL J_ITER>T;

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

      vt=dyy[j_iter,1] - mu - 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_MAT[J_ITER,.]=BETA_TT[3,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,p,varcov,aaa,ccc,e1,e2;

      c1 = c0;


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




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

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


retp(c1);

endp;