/*========= State-Space Model ===========

data : US GDP deflator inflation

<  Model  >
Y(t) = Z(t) + X(t)
Z(t)  = c + Z(t-1) + V(t), V(t) ~ iidN(0, sigz2))
X_t ~ AR(2)

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

new;
cls;
library pgraph optmum;
pqgwin many;

start = 1;
cut = 47;
k = 3;  @ dimension of Beta @
n = 1;  @ dimension of Y(t) @

t = 237;
load data3[t,2] = c:\data\inf\def_3TB.txt;
data = data3[cut+1:t,.];
inf = data[.,1];  @ GDP deflator inflation: 1959:Q1 - 2006:Q4 @
ymat = inf;
t = rows(ymat);

/*=============== Estimation =================*/
prmtr_in = 0.1|0.1|-0.8|-0.8|0.01|0.01; @ starting values @

{xout, fout,gout, cout} = optmum(&lik, prmtr_in);
xout_fnl = trans(xout);
cov = inv(hessp(&lik,xout));
grad = gradfd(&trans,xout);
cov_fnl = grad*cov*grad';
stde = sqrt(diag(cov_fnl));
t_val = xout_fnl./stde;


"======================================";
"    Estimates     S.D.     t value";
"--------------------------------------";
  xout_fnl~abs(stde)~abs(t_val);
"--------------------------------------";
"lnL" -fout;
"--------------------------------------";


/* ============ Filtering =============*/
{trend,cycle} = filter(xout);

seq = seqa(59.25,0.25,t);
graphset;
  _ptitlht=0.25;
  begwind;
  window(3,1,0);

    setwind(1);
    title(" Trend & actual ");
    xy(seq,trend[.,2]~ymat);

   setwind(2);
    title("Cycle");
   xy(seq,cycle~zeros(t,1)); 

   setwind(3);
   title(" Trend ");
   xy(seq,trend[.,1:3]);


endwind;


end;



/*================= Likelihood Function ===================*/
proc  lik(prmtr1);
local  lnf, Ht,prmtr,sigz2,sigx2,Q,R,F,beta_ll,p_ll, 
       lnL,itr, P_tl, eta_tl, f_tl, Kt, beta_tt, P_tt,beta_tl,
       phi1,phi2,c,mu,cor,cov;

prmtr = trans(prmtr1);
phi1 = prmtr[1]; @ AR(1) @
phi2 = prmtr[2]; @ AR(2) @
sigz2 = prmtr[3]; @ Var(permanent shock) @
sigx2 = prmtr[4]; @ Var(transitory shock) @
cor = prmtr[5];  @ correlation @
c = prmtr[6];    @ drift @

cov = cor*sqrt(sigz2*sigx2); @ covariance @

Q = (sigx2~cov~0 )|
    (cov~sigz2~0)|
    (0~0~0 );    @ k by k @     

R=zeros(n,n);   @ n by n @ 

F= 1~0~0|               
   0~phi1~phi2|  
   0~1~0;      @ k by k @

beta_ll = zeros(k,1); @ initial values for the latent variables, wild guess, k by 1 @

P_ll = eye(k)*1000; @ Penalty to the Wild guess for initial values @

Ht = 1~1~0;   @ n by k @
mu = c|0|0;

lnL = 0;
itr = 1;
do until itr>t;

  beta_tl = mu + 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 @

  if itr le start;
   lnf = 0;
  else; 
   lnf = lnpdfmvn(eta_tl',f_tl);
  endif;
  
  lnL = lnL + lnf;

  Kt = P_tl*Ht'*inv(f_tl);  @ Kalman gain, k by k @
  beta_tt =  beta_tl + Kt*eta_tl;  @ k by 1 @
  P_tt = P_tl - P_tl*Ht'*inv(f_tl)*Ht*P_tl; @ k by k @

  beta_ll = beta_tt;
  P_ll = P_tt;

  itr = itr + 1;
endo;

retp(-lnL);
endp;



/*================= Filtering ===================*/
proc(2) = filter(prmtr1);
local  lnf,Ht,prmtr,sigz2,sigx2,Q,R,F,beta_ll,p_ll, 
       lnL,itr, P_tl, eta_tl, f_tl, Kt, beta_tt, P_tt,beta_tl,
       phi1,phi2,c,mu,cor,cov,tend,cycle;

prmtr = trans(prmtr1);
phi1 = prmtr[1]; @ AR(1) @
phi2 = prmtr[2]; @ AR(2) @
sigz2 = prmtr[3]; @ Var(permanent shock) @
sigx2 = prmtr[4]; @ Var(transitory shock) @
cor = prmtr[5];  @ correlation @
c = prmtr[6];    @ drift @

cov = cor*sqrt(sigz2*sigx2);

Q = sigx2~cov~0|
    cov~sigz2~0|
    0~0~0;    @ k by k @     

R = zeros(n,n);   @ n by n @ 

F =  1~0~0|               
   0~phi1~phi2|  
   0~1~0;      @ k by k @

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

P_ll = eye(k)*1000; @ Penalty to the Wild guess for initial values @

Ht = (1~1~0);   @ n by k @
mu = c|0|0;

trend = zeros(t,3);
cycle = zeros(t,1);
itr = 1;
do until itr>t;

  beta_tl = mu + 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 @

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

  trend[itr,1]=beta_tt[1]-1.96*sqrt(P_tt[1,1]);
  trend[itr,2]=beta_tt[1];
  trend[itr,3]=beta_tt[1]+1.96*sqrt(P_tt[1,1]);

  cycle[itr,1]=beta_tt[2];

  beta_ll = beta_tt;
  P_ll = P_tt;

  itr = itr + 1;
endo;

retp(trend,cycle);
endp;


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

  lam1 = prmtr1[1]/(abs(prmtr1[1])+1); @ eigen value 1 @
  lam2 = prmtr1[2]/(abs(prmtr1[2])+1); @ eigen value 2 @

  prmtr[1] = lam1+lam2;  @ phi1 @
  prmtr[2] = -lam1*lam2; @ phi2 @
   
  prmtr[3:4] = exp(prmtr1[3:4]); @ sigz and sigx @
  prmtr[5] = ((prmtr1[5])^5)/(abs(prmtr1[5])^5+1); @ correlation @
 
retp(prmtr);
endp;