/*================================================
         Dynamic factor Model and Kalman Filter 

Y1= r1*X1  + Z1

Y2= r2*X1  + Z2

X1, Z1, Z2 and Z3 ~ AR(1)

=================================================*/

new;
cls;
library pgraph;

/*=============== DGP ============*/
t=200;  @ Sample size @
n=2;    @ # of variables @
k=1;    @ Dimension of beta @

gam=2;
phi=0.5;
sig21=1;
sig22=2;
sig2v=0.8;
tru_para=gam|
         phi|
         sig2v|
         sig21|
         sig22;

beta=zeros(t,1);
b_L=0;
b_var=sig2v/(1-phi^2);
itr=1;
do until itr>t;
   beta[itr]=phi*b_L+sqrt(sig2v)*rndn(1,1);
   b_L=beta[itr];
   itr=itr+1;
endo;

y1=beta+sqrt(sig21)*rndn(t,1);
y2=gam*beta+sqrt(sig22)*rndn(t,1);

 ymat=y1~y2;


/*=============== MLE ============*/
_qn_PrintIters=1; @ print iteration information @

prmtr_in=1|0|ln(1)|ln(1)|ln(1);

{xout,fout,gout,cout}=qnewton(&lik_f, prmtr_in);
xout_fnl=trans(xout);

"-----------------------------------";
"        true values   estimates";
"-----------------------------------";
tru_para~xout_fnl;
"-----------------------------------";

/*=============== filtering  ===============*/
xmat=filter(xout);

@=======Filter=======@
title("Common Factor");
xy(t , xmat~beta~zeros(t,1));


/*================= Lik_f ===================*/
proc lik_f(prmtr1);
local  lnf, Ht,prmtr,Q,R,F,beta_ll,p_ll,lnL,itr, pphi,
       P_tl, eta_tl, f_tl, Kt, beta_tt, P_tt,beta_tl,
       ggam,ssig21,ssig22,ssig2v;

prmtr=trans(prmtr1);
ggam=prmtr[1];
pphi=prmtr[2];
ssig2v=prmtr[3];
ssig21=prmtr[4];
ssig22=prmtr[5];

Q=ssig2v;

R=ssig21~0|0~ssig22;                           @ 1 by 1 @ 

F=pphi;        @ k by k @

beta_ll=zeros(k,1);                @ k by 1 @

P_ll=ssig2v/(1-pphi^2);

Ht=1|ggam;  @ n by k @

lnL=0;
itr=1;
do until itr>t;
  beta_tl=F*beta_ll; @ k by 1 @
  P_tl=F*P_ll*F' +Q;  @ k by k @

 eta_tl=ymat[itr,.]'-Ht*beta_tl; @ n by 1 @
 f_tl=Ht*P_tl*Ht' +R;  @ n by n @
 lnL=lnL - 0.5*ln(2*pi*det(f_tl))-0.5*eta_tl'*inv(f_tl)*eta_tl;

 Kt=P_tl*Ht'*inv(f_tl);
 beta_tt= beta_tl +Kt*eta_tl;
 P_tt=P_tl - Kt*Ht*P_tl;

 beta_ll=beta_tt;
 P_ll=P_tt;

 itr=itr+1;
endo;

retp(-lnL);
endp;


/*================= Filter ===================*/
proc Filter(prmtr1);
local  xmat,lnf, Ht,prmtr,Q,R,F,beta_ll,p_ll,lnL,itr,pphi, 
       P_tl, eta_tl, f_tl, Kt, beta_tt, P_tt,beta_tl,
       ggam,ssig21,ssig22,ssig2v;

prmtr=trans(prmtr1);
ggam=prmtr[1];
pphi=prmtr[2];
ssig2v=prmtr[3];
ssig21=prmtr[4];
ssig22=prmtr[5];

Q=ssig2v;

R=ssig21~0|0~ssig22;                           @ 1 by 1 @ 

F=pphi;        @ k by k @

beta_ll=zeros(k,1);                @ k by 1 @

P_ll=ssig2v/(1-pphi^2);

Ht=1|ggam;  @ n by k @

xmat=zeros(t,1);

lnL=0;
itr=1;
do until itr>t;
  beta_tl=F*beta_ll; @ k by 1 @
  P_tl=F*P_ll*F' +Q;  @ k by k @

 eta_tl=ymat[itr,.]'-Ht*beta_tl; @ n by 1 @
 f_tl=Ht*P_tl*Ht' +R;  @ n by n @
 lnL=lnL - 0.5*ln(2*pi*det(f_tl))-0.5*eta_tl'*inv(f_tl)*eta_tl;

 Kt=P_tl*Ht'*inv(f_tl);
 beta_tt= beta_tl +Kt*eta_tl;
 P_tt=P_tl - Kt*Ht*P_tl;

 xmat[itr]=beta_tt;

 beta_ll=beta_tt;
 P_ll=P_tt;

 itr=itr+1;
endo;

retp(xmat);
endp;



/*================  Trans  ====================*/
proc trans(prmtr1);
local prmtr;
     prmtr=prmtr1;

     prmtr[2]=prmtr1[2]/(abs(prmtr1[2])+1);
     prmtr[3:5]=exp(prmtr1[3:5]);

retp(prmtr);
endp;