/*
    uc_ngfilt.gss, mww, 9/20/2005
    Compute UC trend+irregular model for inflation Series
    Non Gaussian Filter for variances
    
    Use Starting Date = 1953 for GDP Deflator and CPIU
                      = 1960 for PCE_All and PCE_Core
    
    Corrected for error in draw_eps_eta procedure, 4/25.2007, mww
    Added, lower bound on variances 4/26/2007
    
*/
new;
/*
library pgraph;
graphset;
*/
outwidth 200;
@fmtdir="c:\\qepd\\ddisk_finalpaper_corrected\\gss\\fmt\\";@            @ FMT directior @

rndseed 878325;  

@ Parameters @
p_agg=0;         @ Inflation temporal aggregator 0 = average; 3 = last month of quarter @
                   
small=1.0e-06;
big=1.0e+06;

/*
#include fcst.prc; 
#include calendar.gss;          @ Calendars @
#include infl_in.src;           @ Read in Inflation Data @
*/

str_inf =
"cpiu" ;

load data[180,1]=cpiu_q.csv;
p = data[.,1];


meanc(p)~stdc(p);
/*
str_inf =
"gdp" $|
"cpiu" ;
*/
s0 = "cpiu";
vague = 1000;

burnin=100;
ndraw=5100;

@ -- Parameters for log-chi-squared errors -- @
r_p=.086;
r_m1=-7.472;
r_m2=-0.698;
r_sig2=1.411;
r_sig=sqrt(r_sig2);

@ -- Parameters for RW Innovation Variance -- @
tau1=.20;    @ SD regime 1 @
tau2=.20;    @ SD regime 2 @
q_p=1.0;     @ Prob of Regime 1 @
q1=tau1^2;
q2=tau2^2;


firstdate = 1960.0;
lastdate = 2004.99;

n = rows(data);
dnobs = n-1;


@ Carry Out Analysis for Each P-Measure @
format /rd 6,2;
  p=ln(p);
  dp=miss(zeros(dnobs,1),0);
  dp[2:dnobs]=p[2:dnobs]-p[1:dnobs-1];
  dp = trimr(dp,1,0);
@   y=selif(dp,ismpl);@
  y=400*dp;  @ PAAR @
n = rows(y);

@ Construct LT Matrix of ones @
ltone=lowmat(ones(n,n));
bsum=lowmat(ones(n+1,n+1));

  
  @ Parameters for Initial Conditions, bounds and so forth @
    @ Set tau0 = average of first few observations @
    tau0=meanc(y[1:4]);
    
    @ Compute variance of dy for use in constructing bounds and initial conditions @
    dy=y[2:rows(y)]-y[1:rows(y)-1];
    var_dy = (stdc(dy))^2;
    
    /* Lower Bounds on variances -- 
      -- theser are needed to keep algorith away from boundary
      -- the values are problem specific: these values work well for US Inflation */
    var_eta_min = 0.015*var_dy;
    var_eps_min = 0.005*var_dy;
    
    @ Initial Values: these values work well for US Inflation @
    var_eps_initial = var_dy/3;
    var_eta_initial = var_dy/3;
    
    "";
    "Carry Out Analysis for series ";;s0;
    "Estimate of tau0 ";;tau0;
    "Lower bound on SD eta = ";;sqrt(var_eta_min);
    "Lower bound on SD eps = ";;sqrt(var_eps_min);  
    "Initial guess of sd_eta = ";;sqrt(var_eta_initial);
    "Initial guess of sd_eps = ";;sqrt(var_eps_initial);     
    
   y=y-tau0; @ Eliminate intial value of tau from analysis @
   
   @ -- step 0 -- @
    @ -- initial values of indicator prob vectors -- @
     r_pt_eps = r_p*ones(n,1);
     q_pt_eps = q_p*ones(n,1);
     r_pt_eta = r_p*ones(n,1);
     q_pt_eta = q_p*ones(n,1);
    @ -- initial value of var_eps_n  and var_eta_n @
     var_eps_n = var_eps_initial*ones(n,1);
     var_eta_n = var_eta_initial*ones(n,1);
    
    sd_eps_save=zeros(n,ndraw-burnin);
    sd_eta_save=zeros(n,ndraw-burnin);
    tau_save=zeros(n,ndraw-burnin);
    itmp=0;
    idraw=1; do while idraw <= ndraw;
     {eps,eta,tau}=draw_eps_eta(y,var_eps_n,var_eta_n);
     {var_eps_n,r_pt_eps,q_pt_eps}=draw_var(eps,r_pt_eps,q_pt_eps,var_eps_min);
     {var_eta_n,r_pt_eta,q_pt_eta}=draw_var(eta,r_pt_eta,q_pt_eta,var_eta_min);  
     if idraw .> burnin;
      sd_eps_n=sqrt(var_eps_n);
      sd_eta_n=sqrt(var_eta_n);
      sd_eps_save[.,idraw-burnin]=sd_eps_n;
      sd_eta_save[.,idraw-burnin]=sd_eta_n;     
      tau_save[.,idraw-burnin]=tau;
     endif;
     itmp=itmp+1;
     if itmp .== 100;
      itmp=0;
      print /flush;;idraw;;timestr(time);
     endif;
   idraw=idraw+1; endo;  

   m_eta = meanc(sd_eta_save'); s_eta = stdc(sd_eta_save');
   m_eps = meanc(sd_eps_save'); s_eps = stdc(sd_eps_save');
   m_tau = meanc(tau_save');    s_tau = stdc(tau_save');

   m_eta~s_eta~m_eps~s_eps~m_tau~s_tau;
   
@ PROCS @
@ ----------------------------------- @
proc(3) = draw_eps_eta(y,var_eps_n,var_eta_n);
 @ Construct draws of eps and eta @
 local n, cov_eps, cov_tau, diag_y, cov_y, tau, kappa, mutau_y, covtau_y,chol_covtau_y,eta,eps;
 n=rows(y);
 cov_eps=diagrv(eye(n),var_eps_n);
 cov_tau=ltone*cov_eps*ltone'; 
 diag_y=diag(cov_tau)+var_eta_n;
 cov_y=diagrv(cov_tau,diag_y);
 kappa=cov_tau*invpd(cov_y);
 mutau_y=kappa*y;
 covtau_y=cov_tau-kappa*cov_tau';
 chol_covtau_y=chol(covtau_y);
 tau=mutau_y+chol_covtau_y'rndn(n,1);
 eta=y-tau;
 eps=tau[1]|(tau[2:n]-tau[1:n-1]);
 retp(eps,eta,tau);
endp;
@ ----------------------------------------------------@
proc(3) = draw_var(x,r_pt,q_pt,var_min);

@ Construct draw of variance @
local n,lnres2,tmp,ir,iq,mut,qt,vd,valpha,vy,cy,diagvy,vyi,kgain,ye;
local ahat0,ahat1,vhat0,cvhat0,adraw0,adraw1,edraw,udraw;
local f1,f2,fe,fu,vardraw;
    
n=rows(x);
 @ Construct estimate of var_eps @
 lnres2=ln(x.^2);
     
 @ -- Step 1 -- initial draws of Indicator Variables -- @
 tmp=rndu(n,1);
 ir=tmp .<= r_pt;
 tmp=rndu(n,1);
 iq=tmp .<= q_pt;


 @ -- Step 2; compute system parameters given indicators -- @
 mut = (ir*r_m1) + ((1-ir)*r_m2);
 qt = (iq*q1) + ((1-iq)*q2);

 @ -- Compute Covariance Matrix  -- @
 vd=diagrv(eye(n+1),(vague|qt));
 valpha=bsum*vd*bsum';
 vy=valpha[2:n+1,2:n+1];
 cy=valpha[1:n+1,2:n+1];
 diagvy=diag(vy)+r_sig2;
 vy=diagrv(vy,diagvy);
 vyi=invpd(vy);
 kgain=cy*vyi;

 @ -- Compute draws of state and shocks -- @
 ye=lnres2-mut;
 ahat0=kgain*ye;
 ahat1=ahat0[2:n+1];
 vhat0=valpha-kgain*cy';
 cvhat0=chol(vhat0);
 adraw0=ahat0+cvhat0'rndn(n+1,1);
 adraw1=adraw0[2:n+1];
 edraw=lnres2-adraw1;
 udraw=adraw0[2:n+1]-adraw0[1:n];

 @ -- Compute Mixture Probabilities -- @
  @ -- e shocks -- (Note sigma is the same) -- @
  f1=exp(   (-0.5)* (((edraw-r_m1)./r_sig).^2)  );
  f2=exp(   (-0.5)* (((edraw-r_m2)./r_sig).^2)  );
  fe= r_p*f1 + (1-r_p)*f2;
  r_pt=(r_p*f1)./fe;
  
  @ -- u shocks -- Means are both zero @
  f1=(1/tau1)*exp(   (-0.5)* ((udraw./tau1).^2)  );
  f2=(1/tau1)*exp(   (-0.5)* ((udraw./tau2).^2)  );
  fu= q_p*f1 + (1-q_p)*f2;
  q_pt=(q_p*f1)./fu;
    
 @ -- Compute Variance of AR Shocks and Variance of Annual Difs -- @
  vardraw = exp(adraw1); @ Variance Draw @;
  
 @ -- Impose minimum value --- @
 tmp=var_min*ones(rows(vardraw),1);
 vardraw=maxc((vardraw~tmp)');
  
 retp(vardraw,r_pt,q_pt);
endp;   