% ucsv_agg.m -- aggregate inflation
% June 5, 2015, MWW
%


clear all;
small = 1.0e-10;
big = 1.0e+6;                  
rng(663436);
 
  % -- File Directories  
  outdir = '/Users/mwatson/Dropbox/p_comp/ddisk/matlab/out/';
  figdir = '/Users/mwatson/Dropbox/p_comp/ddisk/matlab/fig/';
  matdir = '/Users/mwatson/Dropbox/p_comp/ddisk/matlab/mat/';
    
  % -- Read in Data --- 
  load_data = 1;  % 1 if reloading data from Excel, etc 
  pcomp_data_calendar_mtoq;

  % Dates
  first_date = [1960 1];
  last_date = [2015 1];
  ismpl = smpl(calvec,first_date,last_date,4);
  calvec_ismpl = calvec(ismpl==1);
  dnobs_ismpl = size(calvec_ismpl,1);
  str_tmp = [matdir 'calvec_ismpl_ucsv']; save(str_tmp,'calvec_ismpl');
  
  % Parameters for UCSV Draws
  n_burnin = 10000;   % Discarded Draws
  n_draws = 100000;   % Number of Draws after burnin
  
  % Parameters for scale mixture of epsilon component
  scl_eps_vec = [1;5;10];
  prob_scl_eps_vec = [0.975;(1/60);(1/120)];
  prob_scl_eps_vec(1) = 1-sum(prob_scl_eps_vec(2:end));
  if sum(prob_scl_eps_vec) ~= 1;
      error('Prob_scl_eps_vec does not sum to one');
  end;
  
  % -- Parameters for RW Innovation Variance -- 
  g_values = linspace(0,0.5,11)';
  p_g_values = ones(size(g_values,1),1)/size(g_values,1);
  g_eps_prior = [g_values p_g_values];
  g_dtau_prior = [g_values p_g_values];
   
% ulabel = ['_PCEagg'];
% y = dp_agg(ismpl==1,:);
ulabel = ['_PCEagg_xfe'];
y = dp_agg_xfe(ismpl==1,:);
% ulabel = ['_PCEagg_xe'];
% y = dp_agg_xe(ismpl==1,:);
 

  [tau_draws,tau_f_draws,sigma_dtau_draws,sigma_eps_draws,g_eps_draws,g_dtau_draws,scl_eps_draws,sigmatotal_eps_draws] = ucsv_outlier(y,n_burnin,n_draws,g_eps_prior,g_dtau_prior,scl_eps_vec,prob_scl_eps_vec);
 
  % IMA(1,1) Parameters
  lam0=(sigma_dtau_draws.^2)+2*(sigma_eps_draws.^2);
  rho1=-(sigma_eps_draws.^2)./lam0;
  theta_draws=-(ones(size(rho1))-sqrt(ones(size(rho1))-4*rho1.^2))./(2*rho1);
  var_a=lam0./(ones(size(rho1))+theta_draws.^2);
  sigma_a_draws=sqrt(var_a);
 
  % Compute Posterior for g's
  g_eps_prior_post = g_summary(g_eps_prior,g_eps_draws);
  g_dtau_prior_post = g_summary(g_dtau_prior,g_dtau_draws);
  str_tmp = [matdir 'g_eps_prior_post' ulabel]; save(str_tmp,'g_eps_prior_post');
  str_tmp = [matdir 'g_dtau_prior_post' ulabel]; save(str_tmp,'g_dtau_prior_post');
 
  % Compute posterior for variances
  pctvec = [0.05 1/6 0.50 5/6 0.95]';
  tmp = sigma_dtau_draws.^2;
  var_dtau_mean_pct = post_mean_pct(tmp',pctvec);
  tmp = sigma_eps_draws.^2;
  var_eps_mean_pct = post_mean_pct(tmp',pctvec);
  tmp = sigmatotal_eps_draws.^2;
  vartotal_eps_mean_pct = post_mean_pct(tmp',pctvec);
  str_tmp = [matdir 'var_eps_mean_pct' ulabel]; save(str_tmp,'var_eps_mean_pct');
  str_tmp = [matdir 'var_dtau_mean_pct' ulabel]; save(str_tmp,'var_dtau_mean_pct');
  str_tmp = [matdir 'vartotal_eps_mean_pct' ulabel]; save(str_tmp,'vartotal_eps_mean_pct');
 
  % Compute posteriors for trend values
  tau_mean_pct = post_mean_pct(tau_draws',pctvec);
  tau_f_mean_pct = post_mean_pct(tau_f_draws',pctvec);
  str_tmp = [matdir 'tau_mean_pct' ulabel]; save(str_tmp,'tau_mean_pct');
  str_tmp = [matdir 'tau_f_mean_pct' ulabel]; save(str_tmp,'tau_f_mean_pct');
  tau_f_mean = mean(tau_f_draws,2);
  tmp = mean(tau_f_draws.^2,2);
  tau_f_var = tmp-tau_f_mean.^2;
  tau_f_sd = sqrt(tau_f_var);
  tau_f_mean_sd = [tau_f_mean tau_f_sd];
  str_tmp = [matdir 'tau_f_mean_sd' ulabel]; save(str_tmp,'tau_f_mean_sd');
  
  tau_mean = mean(tau_draws,2);
  tmp = mean(tau_draws.^2,2);
  tau_var = tmp-tau_mean.^2;
  tau_sd = sqrt(tau_var);
  tau_mean_sd = [tau_mean tau_sd];
  str_tmp = [matdir 'tau_mean_sd' ulabel]; save(str_tmp,'tau_mean_sd');
 
  % Compute posteriors for MA parameter
  theta_mean_pct = post_mean_pct(theta_draws',pctvec);
  sigma_a_mean_pct = post_mean_pct(sigma_a_draws',pctvec);
  str_tmp = [matdir 'theta_mean_pct' ulabel]; save(str_tmp,'theta_mean_pct');
  str_tmp = [matdir 'sigma_a_mean_pct' ulabel]; save(str_tmp,'sigma_a_mean_pct');
 
  % Compute posteriors for scale_eps_values
  pctvec = [1/6 0.50 5/6]';
  scl_eps_pct = post_mean_pct(scl_eps_draws',pctvec);
  str_tmp = [matdir 'scl_eps_pct' ulabel]; save(str_tmp,'scl_eps_pct');