% Summary Calculations 
% May 17, 2015
%


clear all;
small = 1.0e-10;
big = 1.0e+6;                  
rng(63761);
 
  % -- File Directories  
  outdir = '/Users/mwatson/Dropbox/p_comp/Matlab/out/';
  figdir = '/Users/mwatson/Dropbox/p_comp//Matlab/fig/';
  matdir = '/Users/mwatson/Dropbox/p_comp/Matlab/mat/';
  mat_save_dir = '/Users/mwatson/Dropbox/p_co/Matlab/mat_save/';  % Location for Saving Various Files
  
  % -- Read in Data --- 
  load_data = 0;  % 1 if reloading data from Excel, etc 
  pcomp_data_calendar_mtoq;

  % Dates
  first_date = [1960 1];
  last_date = [2015 1];
  ismpl_60_end = smpl(calvec,first_date,last_date,4);
  cal_60_end = calvec(ismpl_60_end==1);
  dnobs_60_end = size(cal_60_end,1);
  
  dp_agg_60_end = dp_agg(ismpl_60_end==1,:);
  dp_agg_xe_60_end = dp_agg_xe(ismpl_60_end==1,:);
  dp_agg_xfe_60_end = dp_agg_xfe(ismpl_60_end==1,:);
  dp_disagg_60_end = dp_disagg(ismpl_60_end==1,:);
  share_avg_60_end = share_avg(ismpl_60_end==1,:);
  
  tmp = dp_agg;
  tmp1 = NaN*zeros(dnobs,1);
  tmp1(4:end) = (tmp(4:end)+tmp(3:end-1)+tmp(2:end-2)+tmp(1:end-3))/4;
  dp4_agg_60_end = tmp1(ismpl_60_end==1);
  
  tmp = dp_agg_xe;
  tmp1 = NaN*zeros(dnobs,1);
  tmp1(4:end) = (tmp(4:end)+tmp(3:end-1)+tmp(2:end-2)+tmp(1:end-3))/4;
  dp4_agg_xe_60_end = tmp1(ismpl_60_end==1);
  
  tmp = dp_agg_xfe;
  tmp1 = NaN*zeros(dnobs,1);
  tmp1(4:end) = (tmp(4:end)+tmp(3:end-1)+tmp(2:end-2)+tmp(1:end-3))/4;
  dp4_agg_xfe_60_end = tmp1(ismpl_60_end==1);
  
  tmp = dp_disagg;
  tmp1 = NaN*zeros(dnobs,17);
  tmp1(4:end,:) = (tmp(4:end,:)+tmp(3:end-1,:)+tmp(2:end-2,:)+tmp(1:end-3,:))/4;
  dp4_disagg_60_end = tmp1(ismpl_60_end==1,:);
  
  % -- Univariate Filtered and Smoothed Values
  ulabel = ['_PCEagg'];
  str_tmp = [matdir 'tau_mean_pct' ulabel]; load(str_tmp);
  str_tmp = [matdir 'tau_f_mean_pct' ulabel]; load(str_tmp);
  tau_f_h = tau_f_mean_pct(:,1);
  tau_s_h = tau_mean_pct(:,1);
  tau_f_mean_pct_h = tau_f_mean_pct;
  tau_s_mean_pct_h = tau_mean_pct;
  
  ulabel = ['_PCEagg_xe'];
  str_tmp = [matdir 'tau_mean_pct' ulabel]; load(str_tmp);
  str_tmp = [matdir 'tau_f_mean_pct' ulabel]; load(str_tmp);
  tau_f_xe = tau_f_mean_pct(:,1);
  tau_s_xe = tau_mean_pct(:,1);
  tau_f_mean_pct_xe = tau_f_mean_pct;
  tau_s_mean_pct_xe = tau_mean_pct;
  
  ulabel = ['_PCEagg_xfe'];
  str_tmp = [matdir 'tau_mean_pct' ulabel]; load(str_tmp);
  str_tmp = [matdir 'tau_f_mean_pct' ulabel]; load(str_tmp);
  tau_f_xfe = tau_f_mean_pct(:,1);
  tau_s_xfe = tau_mean_pct(:,1);tau_f_mean_pct_xfe = tau_f_mean_pct;
  tau_s_mean_pct_xfe = tau_mean_pct;
  tau_f_mean_pct_xfe = tau_f_mean_pct;
  tau_s_mean_pct_xfe = tau_mean_pct;
  
  % -- Some results for the 3-component model
  mlabel = '_mucsv_fl_tvp_3comp';
  str_tmp = [matdir 'calvec_ismpl' mlabel]; load(str_tmp);
  % check dates
  tmp = calvec_ismpl-cal_60_end;
  if sum(tmp) ~= 0;
    error('Dates not synchronized');
  end;  
  str_tmp = [matdir 'agg_f_tau_mean_sd' mlabel]; load(str_tmp);
  str_tmp = [matdir 'agg_tau_mean_sd' mlabel]; load(str_tmp);
  tau_f_3c = agg_f_tau_mean_sd(:,1);
  tau_s_3c = agg_tau_mean_sd(:,1);
  str_tmp = [matdir 'agg_f_tau_mean_pct' mlabel]; load(str_tmp);
  str_tmp = [matdir 'agg_tau_mean_pct' mlabel]; load(str_tmp);
  tau_f_mean_pct_3c = agg_f_tau_mean_pct;
  tau_s_mean_pct_3c = agg_tau_mean_pct;
  str_tmp = [matdir 'disagg_f_tau_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'disagg_tau_mean' mlabel]; load(str_tmp);
  dis_tau_f_3c = disagg_f_tau_mean;
  dis_tau_s_3c = disagg_tau_mean;
  % -- Mean Values of Parameters
  str_tmp = [matdir 'alpha_eps_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'alpha_tau_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'sigma_eps_common_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'sigma_eps_common_total_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'sigma_eps_unique_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'sigma_eps_unique_total_mean' mlabel]; load(str_tmp); 
  str_tmp = [matdir 'sigma_dtau_common_mean' mlabel]; load(str_tmp); 
  str_tmp = [matdir 'sigma_dtau_unique_mean' mlabel]; load(str_tmp);
  alpha_eps_3c = alpha_eps_mean;
  alpha_tau_3c = alpha_tau_mean;
  sigma_eps_common_3c = sigma_eps_common_mean;
  sigma_eps_common_total_3c = sigma_eps_common_total_mean;
  sigma_eps_unique_3c = sigma_eps_unique_mean;
  sigma_eps_unique_total_3c = sigma_eps_unique_mean;
  sigma_dtau_common_3c = sigma_dtau_common_mean;
  sigma_dtau_unique_3c = sigma_dtau_unique_mean;
  
   
  % -- Some results for the 17-component model
  mlabel = '_mucsv_fl_tvp';
  str_tmp = [matdir 'calvec_ismpl' mlabel]; load(str_tmp);
  % check dates
  tmp = calvec_ismpl-cal_60_end;
  if sum(tmp) ~= 0;
    error('Dates not synchronized');
  end;  
  str_tmp = [matdir 'agg_f_tau_mean_sd' mlabel]; load(str_tmp);
  str_tmp = [matdir 'agg_tau_mean_sd' mlabel]; load(str_tmp);
  tau_f_17c = agg_f_tau_mean_sd(:,1);
  tau_s_17c = agg_tau_mean_sd(:,1);
  str_tmp = [matdir 'agg_f_tau_mean_pct' mlabel]; load(str_tmp);
  str_tmp = [matdir 'agg_tau_mean_pct' mlabel]; load(str_tmp);
  tau_f_mean_pct_17c = agg_f_tau_mean_pct;
  tau_s_mean_pct_17c = agg_tau_mean_pct;
  str_tmp = [matdir 'disagg_f_tau_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'disagg_tau_mean' mlabel]; load(str_tmp);
  dis_tau_f_17c = disagg_f_tau_mean;
  dis_tau_s_17c = disagg_tau_mean;
  % Compute 3-Components from 17
  s = share_avg(ismpl_60_end==1,:);
  ifood = 2;
  ienergy = [10 15];
  icore = [1 3:9 11:14 16:17];
  s_energy = s(:,ienergy);
  s_energy = s_energy./repmat(sum(s_energy,2),1,2);
  s_core = s(:,icore);
  s_core = s_core./repmat(sum(s_core,2),1,14);
  tau_f_17c_food = dis_tau_f_17c(:,ifood);
  tau_f_17c_energy = sum(dis_tau_f_17c(:,ienergy).*s_energy,2);
  tau_f_17c_core = sum(dis_tau_f_17c(:,icore).*s_core,2);
  tau_s_17c_food = dis_tau_s_17c(:,ifood);
  tau_s_17c_energy = sum(dis_tau_s_17c(:,ienergy).*s_energy,2);
  tau_s_17c_core = sum(dis_tau_s_17c(:,icore).*s_core,2); 
  % -- Mean Values of Parameters
  str_tmp = [matdir 'alpha_eps_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'alpha_tau_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'sigma_eps_common_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'sigma_eps_common_total_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'sigma_eps_unique_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'sigma_eps_unique_total_mean' mlabel]; load(str_tmp); 
  str_tmp = [matdir 'sigma_dtau_common_mean' mlabel]; load(str_tmp); 
  str_tmp = [matdir 'sigma_dtau_unique_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'var_eps_common_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'var_eps_common_total_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'var_eps_unique_mean' mlabel]; load(str_tmp);
  str_tmp = [matdir 'var_eps_unique_total_mean' mlabel]; load(str_tmp); 
  str_tmp = [matdir 'var_dtau_common_mean' mlabel]; load(str_tmp); 
  str_tmp = [matdir 'var_dtau_unique_mean' mlabel]; load(str_tmp);
 
  
  tmp = [cal_60_end dp_agg_60_end tau_f_mean_pct_17c(:,3:5) tau_f_mean_pct_h(:,3:5) tau_f_mean_pct_xe(:,3:5) tau_f_mean_pct_xfe(:,3:5)];
  
  prtmat_comma_screen(tmp,'%4.2f','\n');
  
  
%   alpha_eps_17c = alpha_eps_mean;
%   alpha_tau_17c = alpha_tau_mean;
%   sigma_eps_common_17c = sigma_eps_common_mean;
%   sigma_eps_common_total_17c = sigma_eps_common_total_mean;
%   sigma_eps_unique_17c = sigma_eps_unique_mean;
%   sigma_eps_unique_total_17c = sigma_eps_unique_mean;
%   sigma_dtau_common_17c = sigma_dtau_common_mean;
%   sigma_dtau_unique_17c = sigma_dtau_unique_mean;
%   
%   var_eps_common_17c = var_eps_common_mean;
%   var_eps_common_total_17c = var_eps_common_total_mean;
%   var_eps_unique_17c = var_eps_unique_mean;
%   var_eps_unique_total_17c = var_eps_unique_mean;
%   var_dtau_common_17c = var_dtau_common_mean;
%   var_dtau_unique_17c = var_dtau_unique_mean;
%  
%   
%   alpha_eps_17c_energy = sum(alpha_eps_17c(:,ienergy).*s_energy,2);
%   alpha_tau_17c_energy = sum(alpha_tau_17c(:,ienergy).*s_energy,2);
%   sigma_eps_unique_17c_energy = sqrt(sum((sigma_eps_unique_17c(:,ienergy).*s_energy).^2,2));
%   sigma_eps_unique_total_17c_energy = sqrt(sum((sigma_eps_unique_total_17c(:,ienergy).*s_energy).^2,2));
%   sigma_dtau_unique_17c_energy = sqrt(sum((sigma_dtau_unique_17c(:,ienergy).*s_energy).^2,2));
%   
% %   % Plot model parameters for energy -- 3c and 17c models
% %   figure;
% %    subplot(3,2,1);
% %    plot(cal_60_end,alpha_tau_3c(:,3),'- g','LineWidth',4);
% %    hold on;
% %      plot(cal_60_end,alpha_tau_17c_energy,'- m','LineWidth',4);
% %      plot(cal_60_end,alpha_tau_17c(:,ienergy(1)),'-- r','LineWidth',1);
% %      plot(cal_60_end,alpha_tau_17c(:,ienergy(2)),'-- b','LineWidth',1);
% %    hold off;
% %    title('Alpha tau');
% %    legend('3C','17C','17C Gas','17C Utilities');
% %    legend('Location','SouthWest');
% %    ax = gca;
% %    ax.FontSize = 30;
% %    ylim([0 1]);
% %    ax.XLim = [1960 2015];
% %    
% %    subplot(3,2,2);
% %    plot(cal_60_end,alpha_eps_3c(:,3),'- g','LineWidth',4);
% %    hold on;
% %      plot(cal_60_end,alpha_eps_17c_energy,'- m','LineWidth',4);
% %      plot(cal_60_end,alpha_eps_17c(:,ienergy(1)),'-- r','LineWidth',1);
% %      plot(cal_60_end,alpha_eps_17c(:,ienergy(2)),'-- b','LineWidth',1);
% %    hold off;
% %    title('Alpha eps');
% %    ax = gca;
% %    ax.FontSize = 30;
% %    ylim([0 1]);
% %    ax.XLim = [1960 2015];
% %    
% %    subplot(3,2,3);
% %    plot(cal_60_end,sigma_eps_common_total_3c,'- g','LineWidth',4);
% %    hold on;
% %      plot(cal_60_end,sigma_eps_common_total_17c,'- m','LineWidth',4);
% %    hold off;
% %    title('Sigma eps common');
% %    ax = gca;
% %    ax.FontSize = 30;
% %    ax.XLim = [1960 2015];
% %    
% %    subplot(3,2,4);
% %    plot(cal_60_end,sigma_dtau_common_3c,'- g','LineWidth',4);
% %    hold on;
% %      plot(cal_60_end,sigma_dtau_common_17c,'- m','LineWidth',4);
% %    hold off;
% %    title('Sigma dtau common');
% %    ax = gca;
% %    ax.FontSize = 30;
% %    ax.XLim = [1960 2015];
% %    
% %    subplot(3,2,5);
% %    plot(cal_60_end,sigma_eps_unique_total_3c(:,3),'- g','LineWidth',4);
% %    hold on;
% %      plot(cal_60_end,sigma_eps_unique_total_17c_energy,'- m','LineWidth',4);
% %      plot(cal_60_end,sigma_eps_unique_total_17c(:,ienergy(1)),'-- r','LineWidth',1);
% %      plot(cal_60_end,sigma_eps_unique_total_17c(:,ienergy(2)),'-- b','LineWidth',1);
% %    hold off;
% %    title('Sigma eps unique');
% %    ax = gca;
% %    ax.FontSize = 30;
% %    ax.XLim = [1960 2015];
% %    
% %    subplot(3,2,6);
% %    plot(cal_60_end,sigma_dtau_unique_3c(:,3),'- g','LineWidth',4);
% %    hold on;
% %      plot(cal_60_end,sigma_dtau_unique_17c_energy,'- m','LineWidth',4);
% %      plot(cal_60_end,sigma_dtau_unique_17c(:,ienergy(1)),'-- r','LineWidth',1);
% %      plot(cal_60_end,sigma_dtau_unique_17c(:,ienergy(2)),'-- b','LineWidth',1);
% %    hold off;
% %    title('Sigma dtau unique');
% %    ax = gca;
% %    ax.FontSize = 30;
% %    ax.XLim = [1960 2015];
% %    
% %    error('tmp');
%   
%   % Plot 3 Inflation Rates
% %   plot(cal_60_end,dp_agg_60_end,'- k','LineWidth',3);
% %   hold on;
% %    plot(cal_60_end,dp_agg_xe_60_end,'- b','LineWidth',3);
% %    plot(cal_60_end,dp_agg_xfe_60_end,'- r','LineWidth',3);
% %   hold off;
% %   legend('Headline','X-Energy','X-Food&Energy');
% %   legend('Location','NorthEast');
% %   ax = gca;
% %   ax.FontSize = 20;
% %   ax.XLim = [1960 2015];
%   
% %   plot(cal_60_end,dp4_agg_60_end,'- k','LineWidth',3);
% %   hold on;
% %    plot(cal_60_end,dp4_agg_xe_60_end,'- b','LineWidth',3);
% %    plot(cal_60_end,dp4_agg_xfe_60_end,'- r','LineWidth',3);
% %   hold off;
% %   legend('Headline','X-Energy','X-Food&Energy');
% %   legend('Location','NorthEast');
% %   ax = gca;
% %   ax.FontSize = 20;
% %   ax.XLim = [1960 2015];
%   
%   % Plot Trend Inflation Rates
% %   plot(cal_60_end,dp_agg_60_end,'- k','LineWidth',1);
% %   hold on;
% %    plot(cal_60_end,tau_s_h,'- k','LineWidth',3);
% %    plot(cal_60_end,tau_s_xe,'- b','LineWidth',3);
% %    plot(cal_60_end,tau_s_xfe,'- r','LineWidth',3);
% %   hold off;
% %   legend('Headline','Tau-Headline','Tau-X-Energy','TauX-Food&Energy');
% %   legend('Location','NorthEast');
% %   ax = gca;
% %   ax.FontSize = 20;
% %   ax.XLim = [1960 2015];
% 
% %   % Plot Trend Inflation Rates
% %   plot(cal_60_end,dp_agg_60_end,'- k','LineWidth',1);
% %   hold on;
% %    plot(cal_60_end,tau_s_h,'- k','LineWidth',3);
% %    plot(cal_60_end,tau_s_xe,'- b','LineWidth',3);
% %    plot(cal_60_end,tau_s_xfe,'- r','LineWidth',3);
% %    plot(cal_60_end,tau_s_3c,'- g','LineWidth',3);
% %   hold off;
% %   legend('Headline','Tau-Headline','Tau-X-Energy','TauX-Food&Energy','M3comp');
% %   legend('Location','NorthEast');
% %   ax = gca;
% %   ax.FontSize = 20;
% %   ax.XLim = [1960 2015];
% %   
% %   % Plot Trend Inflation Rates
% %   plot(cal_60_end,dp_agg_60_end,'- k','LineWidth',1);
% %   hold on;
% %    plot(cal_60_end,tau_f_h,'- k','LineWidth',3);
% %    plot(cal_60_end,tau_f_xe,'- b','LineWidth',3);
% %    plot(cal_60_end,tau_f_xfe,'- r','LineWidth',3);
% %    plot(cal_60_end,tau_f_3c,'- g','LineWidth',3);
% %   hold off;
% %   legend('Headline','Tau-Headline','Tau-X-Energy','TauX-Food&Energy','M3comp');
% %   legend('Location','NorthEast');
% %   ax = gca;
% %   ax.FontSize = 20;
% %   ax.XLim = [1960 2015];
% 
% %   % Plot Trend Inflation Rates
% %   plot(cal_60_end,dp_agg_60_end,'- k','LineWidth',1);
% %   hold on;
% %    plot(cal_60_end,tau_s_3c,'- g','LineWidth',3);
% %    plot(cal_60_end,tau_s_17c,'- m','LineWidth',3);
% %   hold off;
% %   legend('Headline','M3comp','M17comp');
% %   legend('Location','NorthEast');
% %   ax = gca;
% %   ax.FontSize = 20;
% %   ax.XLim = [1960 2015];
%   
% 
% %   % Plot Trend Inflation Rates
% %   plot(cal_60_end,dp_agg_60_end,'- k','LineWidth',1);
% %   hold on;
% %    plot(cal_60_end,tau_f_3c,'- g','LineWidth',3);
% %    plot(cal_60_end,tau_f_17c,'- m','LineWidth',3);
% %   hold off;
% %   legend('Headline','M3comp','M17comp');
% %   legend('Location','NorthEast');
% %   ax = gca;
% %   ax.FontSize = 20;
% %   ax.XLim = [1960 2015];
% 
%   
% %   plot(cal_60_end,dp_agg_60_end,'- k','LineWidth',1);
% %   hold on;
% %    plot(cal_60_end,tau_f_h,'- k','LineWidth',3);
% %    plot(cal_60_end,tau_f_xe,'- b','LineWidth',3);
% %    plot(cal_60_end,tau_f_xfe,'- r','LineWidth',3);
% %   hold off;
% %   legend('Headline','Tau-Headline','Tau-X-Energy','TauX-Food&Energy');
% %   legend('Location','NorthEast');
% %   ax = gca;
% %   ax.FontSize = 20;
% %   ax.XLim = [1960 2015];
% % 
% % % Plot Components of Trend Inflation -- 3 Component
% %   plot(cal_60_end,tau_s_3c,'- k','LineWidth',1);
% %   hold on;
% %       plot(cal_60_end,dis_tau_s_3c(:,1),'- b','LineWidth',1);
% %       plot(cal_60_end,dis_tau_s_3c(:,2),'- r','LineWidth',1);
% %       plot(cal_60_end,dis_tau_s_3c(:,3),'- g','LineWidth',1);
% %   hold off;
% %   legend('Total','Core','Food','Energy');
% %   legend('Location','NorthEast');
% %   ax = gca;
% %   ax.FontSize = 20;
% %   ax.XLim = [1960 2015];
% %   error('tmp');
%   
% % %  Plot Components of Trend Inflation -- 3 Component
% %   plot(cal_60_end,tau_f_3c,'- k','LineWidth',1);
% %   hold on;
% %       plot(cal_60_end,dis_tau_f_3c(:,1),'- b','LineWidth',1);
% %       plot(cal_60_end,dis_tau_f_3c(:,2),'- r','LineWidth',1);
% %       plot(cal_60_end,dis_tau_f_3c(:,3),'- g','LineWidth',1);
% %   hold off;
% %   legend('Total','Core','Food','Energy');
% %   legend('Location','NorthEast');
% %   ax = gca;
% %   ax.FontSize = 20;
% %   ax.XLim = [1960 2015];
% 
% % %  Plot Components of Trend Inflation -- 3 and 17 Components
% % figure;
% % subplot(2,2,1);
% %   plot(cal_60_end,dis_tau_s_3c(:,1),'- g','LineWidth',4);
% %   hold on;
% %       plot(cal_60_end,tau_s_17c_core,'- m','LineWidth',4);;
% %   hold off;
% %   legend('3C','17C');
% %   legend('Location','NorthWest');
% %   title('Core');
% %   ax = gca;
% %   ax.FontSize = 30;
% %   ax.XLim = [1960 2015];
% % 
% %  subplot(2,2,2);
% %   plot(cal_60_end,dis_tau_s_3c(:,2),'- g','LineWidth',4);
% %   hold on;
% %       plot(cal_60_end,tau_s_17c_food,'- m','LineWidth',4);;
% %   hold off;
% %   legend('3C','17C');
% %   legend('Location','NorthWest');
% %   title('Food');
% %   ax = gca;
% %   ax.FontSize = 30;
% %   ax.XLim = [1960 2015];
% %  
% % subplot(2,2,3);
% %   plot(cal_60_end,dis_tau_s_3c(:,3),'- g','LineWidth',4);
% %   hold on;
% %       plot(cal_60_end,tau_s_17c_energy,'- m','LineWidth',4);;
% %   hold off;
% %   legend('3C','17C');
% %   legend('Location','NorthWest');
% %   title('Energy');
% %   ax = gca;
% %   ax.FontSize = 30;
% %   ax.XLim = [1960 2015];
%   
% % figure;
% % subplot(2,2,1);
% %   plot(cal_60_end,dis_tau_f_3c(:,1),'- g','LineWidth',4);
% %   hold on;
% %       plot(cal_60_end,tau_f_17c_core,'- m','LineWidth',4);;
% %   hold off;
% %   legend('3C','17C');
% %   legend('Location','NorthWest');
% %   title('Core');
% %   ax = gca;
% %   ax.FontSize = 30;
% %   ax.XLim = [1960 2015];
% % 
% %  subplot(2,2,2);
% %   plot(cal_60_end,dis_tau_f_3c(:,2),'- g','LineWidth',4);
% %   hold on;
% %       plot(cal_60_end,tau_f_17c_food,'- m','LineWidth',4);;
% %   hold off;
% %   legend('3C','17C');
% %   legend('Location','NorthWest');
% %   title('Food');
% %   ax = gca;
% %   ax.FontSize = 30;
% %   ax.XLim = [1960 2015];
% %  
% % subplot(2,2,3);
% %   plot(cal_60_end,dis_tau_f_3c(:,3),'- g','LineWidth',4);
% %   hold on;
% %       plot(cal_60_end,tau_f_17c_energy,'- m','LineWidth',4);;
% %   hold off;
% %   legend('3C','17C');
% %   legend('Location','NorthWest');
% %   title('Energy');
% %   ax = gca;
% %   ax.FontSize = 30;
% %   ax.XLim = [1960 2015];
% 
% % Compute Impulse Responses for each Model by date
%   
%  % IRF for Headline
% %  ulabel = ['_PCEagg'];
% %  ydata = dp_agg_60_end;
% %  str_tmp = [matdir 'var_dtau_mean_pct' ulabel]; load(str_tmp);
% %  str_tmp = [matdir 'var_eps_mean_pct' ulabel]; load(str_tmp);
% %  str_tmp = [matdir 'vartotal_eps_mean_pct' ulabel]; load(str_tmp);
% %  
% %  var_dtau = var_dtau_mean_pct(:,4);
% %  var_eps = vartotal_eps_mean_pct(:,4);
% %  notim = dnobs_60_end; 
% %  % Compute Filtered Value
% %  u_tau_decomp = zeros(notim,notim);
% %  for t = 1:notim;
% %           ydata = zeros(notim,1);
% %           ydata(t) = 1;
% %           tmp = ufilt(ydata,var_dtau,var_eps);
% %           tmp = tmp(t:end);
% %           for j = 1:size(tmp,1);
% %               u_tau_decomp(t-1+j,j) = tmp(j);
% %           end;
% %  end;
% %  
% %  irf_h = u_tau_decomp;
% %  
% %  % IRF for XE
% %  ulabel = ['_PCEagg_xe'];
% %  ydata = dp_agg_xe_60_end;
% %  str_tmp = [matdir 'var_dtau_mean_pct' ulabel]; load(str_tmp);
% %  str_tmp = [matdir 'var_eps_mean_pct' ulabel]; load(str_tmp);
% %  str_tmp = [matdir 'vartotal_eps_mean_pct' ulabel]; load(str_tmp);
% %  
% %  var_dtau = var_dtau_mean_pct(:,4);
% %  var_eps = vartotal_eps_mean_pct(:,4);
% %  notim = dnobs_60_end; 
% %  % Compute Filtered Value
% %  u_tau_decomp = zeros(notim,notim);
% %  for t = 1:notim;
% %           ydata = zeros(notim,1);
% %           ydata(t) = 1;
% %           tmp = ufilt(ydata,var_dtau,var_eps);
% %           tmp = tmp(t:end);
% %           for j = 1:size(tmp,1);
% %               u_tau_decomp(t-1+j,j) = tmp(j);
% %           end;
% %  end;
% %  
% %  irf_xe = u_tau_decomp;
% %  
% %  % IRF for XFE
% %  ulabel = ['_PCEagg_xfe'];
% %  ydata = dp_agg_xfe_60_end;
% %  str_tmp = [matdir 'var_dtau_mean_pct' ulabel]; load(str_tmp);
% %  str_tmp = [matdir 'var_eps_mean_pct' ulabel]; load(str_tmp);
% %  str_tmp = [matdir 'vartotal_eps_mean_pct' ulabel]; load(str_tmp);
% %  
% %  var_dtau = var_dtau_mean_pct(:,4);
% %  var_eps = vartotal_eps_mean_pct(:,4);
% %  notim = dnobs_60_end; 
% %  % Compute Filtered Value
% %  u_tau_decomp = zeros(notim,notim);
% %  for t = 1:notim;
% %           ydata = zeros(notim,1);
% %           ydata(t) = 1;
% %           tmp = ufilt(ydata,var_dtau,var_eps);
% %           tmp = tmp(t:end);
% %           for j = 1:size(tmp,1);
% %               u_tau_decomp(t-1+j,j) = tmp(j);
% %           end;
% %  end;
% %  
% %  irf_xfe = u_tau_decomp;
% %  
% %  
% %  % Plot IRFs over for different time periods
% %  plot_date = [[1965 1];[1975 1];[1985 1];[1995 1];[2005 1];[2015 1]];
% %  n_irf = 7;
% %  ctvec = [0:n_irf-1];
% %  figure;
% %    for i_date = 1:size(plot_date,1);
% %      dvec = plot_date(i_date,:);
% %      ii = smpl(cal_60_end,dvec,dvec,4);
% %      i_h = irf_h(ii==1,1:n_irf);
% %      i_xe = irf_xe(ii==1,1:n_irf);
% %      i_xfe = irf_xfe(ii==1,1:n_irf);   
% %      str1 = [num2str(dvec(1,1),'%4i:Q') num2str(dvec(1,2),'%1i')];
% %      subplot(3,2,i_date);
% %        plot(ctvec,i_h,'- k','LineWidth',3);
% %        hold on;
% %          plot(ctvec,i_xe,'- b','LineWidth',3);
% %          plot(ctvec,i_xfe,'- r','LineWidth',3);
% %        hold off;
% %        legend('Headline','X-Energy','X-Food&Energy');
% %        legend('Location','NorthEast');
% %        title(str1);
% %        ax = gca;
% %        ax.FontSize = 20;
% %        ax.YLim = [0 1];
% %        ax.XTick = ctvec; 
% %    end;
%  
%  % ------------- Compute Univariate Forecast Errors -----------
%  % Compute Forecast Errors by Horizon
%  hmax = 12;  % Maximum Forecast Horizon;
%  dp_act = dp_agg_60_end;   % Target Series to forecast 
%  %dp_act = dp_agg_xfe_60_end;   % Target Series to forecast 
%  e_h = NaN*zeros(dnobs_60_end,hmax);
%  e_xe = NaN*zeros(dnobs_60_end,hmax);
%  e_xfe = NaN*zeros(dnobs_60_end,hmax);
%  e_3c = NaN*zeros(dnobs_60_end,hmax);
%  e_17c = NaN*zeros(dnobs_60_end,hmax);
%  e_ao_xe = NaN*zeros(dnobs_60_end,hmax);
%  e_ao_h = NaN*zeros(dnobs_60_end,hmax);
%  e_ao_xfe = NaN*zeros(dnobs_60_end,hmax);
%  e_l_h = NaN*zeros(dnobs_60_end,hmax);
%  e_l_xfe = NaN*zeros(dnobs_60_end,hmax);
%  e_l_xe = NaN*zeros(dnobs_60_end,hmax);
%  act_dph = NaN*zeros(dnobs_60_end,hmax);
%  for h = 1:hmax;
%      dph_act = NaN*zeros(dnobs_60_end,1);
%      dph_act(1:end-h) = dp_act(2:end-h+1);
%      for j = 2:h;
%          dph_act(1:end-h) = dph_act(1:end-h)+dp_act(j+1:end-h+j);
%      end;
%      dph_act = dph_act/h;
%      e_h(:,h) = dph_act - tau_f_h;
%      e_xe(:,h) = dph_act - tau_f_xe;
%      e_xfe(:,h) = dph_act - tau_f_xfe;
%      e_3c(:,h) = dph_act - tau_f_3c;
%      e_17c(:,h) = dph_act - tau_f_17c; 
%      e_l_h(:,h) = dph_act - dp_agg_60_end;
%      e_l_xfe(:,h) = dph_act - dp_agg_xfe_60_end;
%      e_l_xe(:,h) = dph_act - dp_agg_xe_60_end;
%      e_ao_h(:,h) = dph_act - dp4_agg_60_end;
%      e_ao_xe(:,h) = dph_act - dp4_agg_xe_60_end;
%      e_ao_xfe(:,h) = dph_act - dp4_agg_xfe_60_end;       
%      act_dph(:,h) = dph_act;
%  end;
%  
%  n_roll = 20;
%  r_e_h = NaN*zeros(dnobs_60_end,hmax);
%  r_e_xe = NaN*zeros(dnobs_60_end,hmax);
%  r_e_xfe = NaN*zeros(dnobs_60_end,hmax);
%  r_e_3c = NaN*zeros(dnobs_60_end,hmax);
%  r_e_17c = NaN*zeros(dnobs_60_end,hmax);
%  r_e_ao_xe = NaN*zeros(dnobs_60_end,hmax);
%  r_e_ao_h = NaN*zeros(dnobs_60_end,hmax);
%  r_e_ao_xfe = NaN*zeros(dnobs_60_end,hmax);
%  r_e_l_h = NaN*zeros(dnobs_60_end,hmax);
%  r_e_l_xfe = NaN*zeros(dnobs_60_end,hmax);
%  r_e_l_xe = NaN*zeros(dnobs_60_end,hmax);
%  
%  for t = 1:dnobs_60_end-n_roll+1;
%      r_e_h(t,:) = sqrt(mean(e_h(t:t+n_roll-1,:).^2)');
%      r_e_xe(t,:) = sqrt(mean(e_xe(t:t+n_roll-1,:).^2)');
%      r_e_xfe(t,:) = sqrt(mean(e_xfe(t:t+n_roll-1,:).^2)');
%      r_e_17c(t,:) = sqrt(mean(e_17c(t:t+n_roll-1,:).^2)');
%  end;
%  
% % Forecast Errors
%   h = 4;
% %   plot(cal_60_end,e_17c(:,h),'- k','LineWidth',3);
% %   hold on;
% %       plot(cal_60_end,e_h(:,h),'- b','LineWidth',1);
% %       plot(cal_60_end,e_xfe(:,h),'- r','LineWidth',1);
% %       plot(cal_60_end,e_xe(:,h),'- g','LineWidth',1);
% %   hold off;
% %   legend('Multivariate','Univ-Headline','Univ-XFE','Univ-XE');
% %   legend('Location','NorthEast');
% %   ax = gca;
% %   ax.FontSize = 30;
% %   ax.XLim = [1965 2015];
% %   error('tmp');
% %   
% %   plot(cal_60_end,r_e_17c(:,h),'- k','LineWidth',3);
% %   hold on;
% %       plot(cal_60_end,r_e_h(:,h),'- b','LineWidth',1);
% %       plot(cal_60_end,r_e_xfe(:,h),'- r','LineWidth',1);
% %       plot(cal_60_end,r_e_xe(:,h),'- g','LineWidth',1);
% %   hold off;
% %   legend('Multivariate','Univ-Headline','Univ-XFE','Univ-XE');
% %   legend('Location','NorthEast');
% %   ax = gca;
% %   ax.FontSize = 30;
% %   ax.XLim = [1965 2015];
% 
% % Forecast Errors in Order
% h = 4;
% e = [e_17c(:,h) e_3c(:,h) e_h(:,h) e_xe(:,h) e_xfe(:,h) e_l_h(:,h) e_l_xe(:,h) e_l_xfe(:,h) e_ao_h(:,h) e_ao_xe(:,h) e_ao_xfe(:,h)]; 
% 
% 
% 
% % Date for Computing forecast errors;
% % i_fcst = smpl(cal_60_end,[1965 1],[2015 1],4);
% % i_fcst = smpl(cal_60_end,[1965 1],[1989 4],4);
% % i_fcst = smpl(cal_60_end,[1990 1],[2015 1],4);
% % cal_fcst = cal_60_end(i_fcst==1);
% % 
% % e = e(i_fcst==1,:);
% % e = packr(e);
% % rmse_e = sqrt(mean(e.^2));
% % prtmat_comma_screen(rmse_e,'%4.2f','\n');
% % 
% % rslt_mat = NaN*zeros(4,size(e,2));
% % nma = 6;
% % for i = 1:size(e,2);
% %     tmp1 = e(:,1).^2;
% %     tmp = e(:,i).^2;
% %     [b,vb] = hacm(tmp,nma,1);
% %     rslt_mat(1,i) = b;
% %     rslt_mat(2,i) = sqrt(vb);
% %     if i > 1;
% %         tmpdif = tmp-tmp1;
% %         [b,vb] = hacm(tmpdif,nma,1);
% %         rslt_mat(3,i) = b;
% %         rslt_mat(4,i) = sqrt(vb);
% %     end;
% % end;
% % 
% % for i = 1:size(e,2);
% %     fprintf('%4.2f',rslt_mat(1,i));
% %     if i < size(e,2);
% %       fprintf('(%4.2f),',rslt_mat(2,i));
% %     else;
% %       fprintf('(%4.2f) \n',rslt_mat(2,i));
% %     end;
% % end;
% % 
% % for i = 1:size(e,2);
% %     fprintf('%4.2f ',rslt_mat(3,i));
% %     if i < size(e,2);
% %       fprintf('(%4.2f),',rslt_mat(4,i));
% %     else;
% %       fprintf('(%4.2f) \n',rslt_mat(4,i));
% %     end;
% % end;
% 
% tmp = [cal_60_end
% 
%     
% 
% 
% % % 
% % %  % Plot actual and forecast errors for particular horizon
% %  h = 4;
% %  
% %  % Headline 
% %  tau_f = tau_f_h;
% %  e = e_h;
% %  figure 
% %    subplot(2,1,1);
% %     plot(cal_60_end(i_fcst==1),act_dph(i_fcst==1,h),'- k','LineWidth',1);
% %     hold on;
% %      plot(cal_60_end(i_fcst==1),tau_f(i_fcst==1),'- k','LineWidth',3);
% %     hold off;
% %     ylim([-2 12]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Actual','Forecast');
% %     legend('Location','NorthEast');
% %    subplot(2,1,2);
% %     plot(cal_60_end(i_fcst==1),e(i_fcst==1,h),'- k','LineWidth',3);
% %     ylim([-5 5]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Forecast Error');
%    
% %  % XE
% %  tau_f = tau_f_xe;
% %  e = e_xe;
% %  figure 
% %    subplot(2,1,1);
% %     plot(cal_60_end(i_fcst==1),act_dph(i_fcst==1,h),'- k','LineWidth',1);
% %     hold on;
% %      plot(cal_60_end(i_fcst==1),tau_f(i_fcst==1),'- b','LineWidth',3);
% %     hold off;
% %     ylim([-2 12]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Actual','Forecast');
% %     legend('Location','NorthEast');
% %    subplot(2,1,2);
% %     plot(cal_60_end(i_fcst==1),e(i_fcst==1,h),'- b','LineWidth',3);
% %     ylim([-5 5]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Forecast Error');
%     
% %  % XFE
% %  tau_f = tau_f_xfe;
% %  e = e_xfe;
% %  figure 
% %    subplot(2,1,1);
% %     plot(cal_60_end(i_fcst==1),act_dph(i_fcst==1,h),'- k','LineWidth',1);
% %     hold on;
% %      plot(cal_60_end(i_fcst==1),tau_f(i_fcst==1),'- r','LineWidth',3);
% %     hold off;
% %     ylim([-2 12]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Actual','Forecast');
% %     legend('Location','NorthEast');
% %    subplot(2,1,2);
% %     plot(cal_60_end(i_fcst==1),e(i_fcst==1,h),'- r','LineWidth',3);
% %     ylim([-5 5]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Forecast Error');
% 
% %  % 3c
% %  tau_f = tau_f_3c;
% %  e = e_3c;
% %  figure 
% %    subplot(2,1,1);
% %     plot(cal_60_end(i_fcst==1),act_dph(i_fcst==1,h),'- k','LineWidth',1);
% %     hold on;
% %      plot(cal_60_end(i_fcst==1),tau_f(i_fcst==1),'- g','LineWidth',3);
% %     hold off;
% %     ylim([-2 12]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Actual','Forecast');
% %     legend('Location','NorthEast');
% %    subplot(2,1,2);
% %     plot(cal_60_end(i_fcst==1),e(i_fcst==1,h),'- g','LineWidth',3);
% %     ylim([-5 5]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Forecast Error');
% %     
% 
% %  % 17c
% %  tau_f = tau_f_17c;
% %  e = e_17c;
% %  figure 
% %    subplot(2,1,1);
% %     plot(cal_60_end(i_fcst==1),act_dph(i_fcst==1,h),'- k','LineWidth',1);
% %     hold on;
% %      plot(cal_60_end(i_fcst==1),tau_f(i_fcst==1),'- m','LineWidth',3);
% %     hold off;
% %     ylim([-2 12]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Actual','Forecast');
% %     legend('Location','NorthEast');
% %    subplot(2,1,2);
% %     plot(cal_60_end(i_fcst==1),e(i_fcst==1,h),'- m','LineWidth',3);
% %     ylim([-5 5]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Forecast Error');
% %     error('tmp');
% % 
% %  % l_xe
% %  tau_f = dp_agg_xe_60_end;
% %  e = e_l_xe;
% %  figure 
% %    subplot(2,1,1);
% %     plot(cal_60_end(i_fcst==1),act_dph(i_fcst==1,h),'- k','LineWidth',1);
% %     hold on;
% %      plot(cal_60_end(i_fcst==1),tau_f(i_fcst==1),'- b','LineWidth',4);
% %     hold off;
% %     ylim([-2 12]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Actual','Forecast');
% %     legend('Location','NorthEast');
% %    subplot(2,1,2);
% %     plot(cal_60_end(i_fcst==1),e(i_fcst==1,h),'- b','LineWidth',4);
% %     ylim([-5 5]);
% %     ax = gca;
% %     ax.FontSize = 30;
% %     ax.XLim = [1960 2015];
% %     legend('Forecast Error');
% %     error('tmp');
%     
% %  % l_xfe
% %  tau_f = dp_agg_xfe_60_end;
% %  e = e_l_xfe;
% %  figure 
% %    subplot(2,1,1);
% %     plot(cal_60_end(i_fcst==1),act_dph(i_fcst==1,h),'- k','LineWidth',1);
% %     hold on;
% %      plot(cal_60_end(i_fcst==1),tau_f(i_fcst==1),'- b','LineWidth',4);
% %     hold off;
% %     ylim([-2 12]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Actual','Forecast');
% %     legend('Location','NorthEast');
% %    subplot(2,1,2);
% %     plot(cal_60_end(i_fcst==1),e(i_fcst==1,h),'- b','LineWidth',4);
% %     ylim([-5 5]);
% %     ax = gca;
% %     ax.FontSize = 30;
% %     ax.XLim = [1960 2015];
% %     legend('Forecast Error');
% %     error('tmp');
%     
% %  % ao_h
% %  tau_f = dp4_agg_60_end;
% %  e = e_ao_h;
% %  figure 
% %    subplot(2,1,1);
% %     plot(cal_60_end(i_fcst==1),act_dph(i_fcst==1,h),'- k','LineWidth',1);
% %     hold on;
% %      plot(cal_60_end(i_fcst==1),tau_f(i_fcst==1),'- b','LineWidth',4);
% %     hold off;
% %     ylim([-2 12]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Actual','Forecast');
% %     legend('Location','NorthEast');
% %    subplot(2,1,2);
% %     plot(cal_60_end(i_fcst==1),e(i_fcst==1,h),'- b','LineWidth',4);
% %     ylim([-5 5]);
% %     ax = gca;
% %     ax.FontSize = 30;
% %     ax.XLim = [1960 2015];
% %     legend('Forecast Error');
% %     error('tmp');
%     
% % % ao_xe
% %  tau_f = dp4_agg_xe_60_end;
% %  e = e_ao_xe;
% %  figure 
% %    subplot(2,1,1);
% %     plot(cal_60_end(i_fcst==1),act_dph(i_fcst==1,h),'- k','LineWidth',1);
% %     hold on;
% %      plot(cal_60_end(i_fcst==1),tau_f(i_fcst==1),'- b','LineWidth',4);
% %     hold off;
% %     ylim([-2 12]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Actual','Forecast');
% %     legend('Location','NorthEast');
% %    subplot(2,1,2);
% %     plot(cal_60_end(i_fcst==1),e(i_fcst==1,h),'- b','LineWidth',4);
% %     ylim([-5 5]);
% %     ax = gca;
% %     ax.FontSize = 30;
% %     ax.XLim = [1960 2015];
% %     legend('Forecast Error');
% %     error('tmp');
%     
% % % ao_xfe
% %  tau_f = dp4_agg_xfe_60_end;
% %  e = e_ao_xfe;
% %  figure 
% %    subplot(2,1,1);
% %     plot(cal_60_end(i_fcst==1),act_dph(i_fcst==1,h),'- k','LineWidth',1);
% %     hold on;
% %      plot(cal_60_end(i_fcst==1),tau_f(i_fcst==1),'- b','LineWidth',4);
% %     hold off;
% %     ylim([-2 12]);
% %     ax = gca;
% %     ax.FontSize = 20;
% %     ax.XLim = [1960 2015];
% %     legend('Actual','Forecast');
% %     legend('Location','NorthEast');
% %    subplot(2,1,2);
% %     plot(cal_60_end(i_fcst==1),e(i_fcst==1,h),'- b','LineWidth',4);
% %     ylim([-5 5]);
% %     ax = gca;
% %     ax.FontSize = 30;
% %     ax.XLim = [1960 2015];
% %     legend('Forecast Error');
% %     error('tmp')nf
% 
% % nf = 10;
% % me_e = NaN*zeros(nf,1);
% % rmse_e = NaN*zeros(nf,1);
% % mae_e = NaN*zeros(nf,1);
% % min_e = NaN*zeros(nf,1);
% % max_e = NaN*zeros(nf,1);
% % nobs_e = NaN*zeros(5,1);
% % pct_vec = [0.01 0.05 0.167 0.50 0.833 0.95 0.99]';
% % pct_e = NaN*ones(nf,size(pct_vec,1));
% % 
% % figure 
% % edges = [-7 -5:0.5:5 7];
% %  j = 1;
% %  e = e_h(i_fcst==1,h);
% %  e = packr(e);
% %  subplot(5,2,j);
% %    histogram(e,edges);
% %    title('Headline');
% %    ax = gca;
% %    ax.FontSize = 20;
% %    me_e(j) = mean(e);
% %    rmse_e(j) = sqrt(mean(e.^2));
% %    min_e(j) = min(e);
% %    max_e(j) = max(e);
% %    mae_e(j) = mean(abs(e));
% %    nobs_e(j) = size(e,1);
% %    pct_vec = [0.01 0.05 0.167 0.50 0.833 0.95 0.99]';
% %    pct = pctile(e,pct_vec);
% %    pct_e(j,:) = pct';
% %  
% %  j = 2;
% %  e = e_xe(i_fcst==1,h);
% %  e = packr(e);
% %  subplot(5,2,j);
% %    histogram(e,edges);
% %    title('XE');
% %    ax = gca;
% %    ax.FontSize = 20;
% %    me_e(j) = mean(e);
% %    rmse_e(j) = sqrt(mean(e.^2));
% %    min_e(j) = min(e);
% %    max_e(j) = max(e);
% %    mae_e(j) = mean(abs(e));
% %    nobs_e(j) = size(e,1);
% %    pct_vec = [0.01 0.05 0.167 0.50 0.833 0.95 0.99]';
% %    pct = pctile(e,pct_vec);
% %    pct_e(j,:) = pct';
% %  
% %  j = 3;
% %  e = e_xfe(i_fcst==1,h);
% %  e = packr(e);
% %  subplot(5,2,j);
% %    histogram(e,edges);
% %    title('XFE');
% %    ax = gca;
% %    ax.FontSize = 20;
% %    me_e(j) = mean(e);
% %    rmse_e(j) = sqrt(mean(e.^2));
% %    min_e(j) = min(e);
% %    max_e(j) = max(e);
% %    mae_e(j) = mean(abs(e));
% %    nobs_e(j) = size(e,1);
% %    pct = pctile(e,pct_vec);
% %    pct_e(j,:) = pct';
% %    
% %  j = 4;
% %  e = e_3c(i_fcst==1,h);
% %  e = packr(e);
% %  subplot(5,2,j);
% %    histogram(e,edges);
% %    title('3c');
% %    ax = gca;
% %    ax.FontSize = 20;
% %    me_e(j) = mean(e);
% %    rmse_e(j) = sqrt(mean(e.^2));
% %    min_e(j) = min(e);
% %    max_e(j) = max(e);
% %    mae_e(j) = mean(abs(e));
% %    nobs_e(j) = size(e,1);
% %    pct = pctile(e,pct_vec);
% %    pct_e(j,:) = pct';
% %    
% %  j = 5;
% %  e = e_17c(i_fcst==1,h);
% %  e = packr(e);
% %  subplot(5,2,j);
% %    histogram(e,edges);
% %    title('17c');
% %    ax = gca;
% %    ax.FontSize = 20;
% %    me_e(j) = mean(e);
% %    rmse_e(j) = sqrt(mean(e.^2));
% %    min_e(j) = min(e);
% %    max_e(j) = max(e);
% %    mae_e(j) = mean(abs(e));
% %    nobs_e(j) = size(e,1);
% %    pct = pctile(e,pct_vec);
% %    pct_e(j,:) = pct';
% %    
% %  j = 6;
% %  e = e_l_xe(i_fcst==1,h);
% %  e = packr(e);
% %  subplot(5,2,j);
% %    histogram(e,edges);
% %    title('level xe');
% %    ax = gca;
% %    ax.FontSize = 20;
% %    me_e(j) = mean(e);
% %    rmse_e(j) = sqrt(mean(e.^2));
% %    min_e(j) = min(e);
% %    max_e(j) = max(e);
% %    mae_e(j) = mean(abs(e));
% %    nobs_e(j) = size(e,1);
% %    pct = pctile(e,pct_vec);
% %    pct_e(j,:) = pct';
% %    
% %  j = 7;
% %  e = e_l_xfe(i_fcst==1,h);
% %  e = packr(e);
% %  subplot(5,2,j);
% %    histogram(e,edges);
% %    title('level xfe');
% %    ax = gca;
% %    ax.FontSize = 20;
% %    me_e(j) = mean(e);
% %    rmse_e(j) = sqrt(mean(e.^2));
% %    min_e(j) = min(e);
% %    max_e(j) = max(e);
% %    mae_e(j) = mean(abs(e));
% %    nobs_e(j) = size(e,1);
% %    pct = pctile(e,pct_vec);
% %    pct_e(j,:) = pct';  
% %    
% %  j = 8;
% %  e = e_ao_h(i_fcst==1,h);
% %  e = packr(e);
% %  subplot(5,2,j);
% %    histogram(e,edges);
% %    title('AO h');
% %    ax = gca;
% %    ax.FontSize = 20;
% %    me_e(j) = mean(e);
% %    rmse_e(j) = sqrt(mean(e.^2));
% %    min_e(j) = min(e);
% %    max_e(j) = max(e);
% %    mae_e(j) = mean(abs(e));
% %    nobs_e(j) = size(e,1);
% %    pct = pctile(e,pct_vec);
% %    pct_e(j,:) = pct';  
% %    
% %  j = 9;
% %  e = e_ao_xe(i_fcst==1,h);
% %  e = packr(e);
% %  subplot(5,2,j);
% %    histogram(e,edges);
% %    title('AO xe');
% %    ax = gca;
% %    ax.FontSize = 20;
% %    me_e(j) = mean(e);
% %    rmse_e(j) = sqrt(mean(e.^2));
% %    min_e(j) = min(e);
% %    max_e(j) = max(e);
% %    mae_e(j) = mean(abs(e));
% %    nobs_e(j) = size(e,1);
% %    pct = pctile(e,pct_vec);
% %    pct_e(j,:) = pct'; 
% %    
% %  j = 10;
% %  e = e_ao_xfe(i_fcst==1,h);
% %  e = packr(e);
% %  subplot(5,2,j);
% %    histogram(e,edges);
% %    title('AO xfe');
% %    ax = gca;
% %    ax.FontSize = 20;
% %    me_e(j) = mean(e);
% %    rmse_e(j) = sqrt(mean(e.^2));
% %    min_e(j) = min(e);
% %    max_e(j) = max(e);
% %    mae_e(j) = mean(abs(e));
% %    nobs_e(j) = size(e,1);
% %    pct = pctile(e,pct_vec);
% %    pct_e(j,:) = pct'; 
% % 
% %   fprintf('nobs,'); prtmat_comma_screen(nobs_e','%3i','\n');
% %   fprintf('mean,'); prtmat_comma_screen(me_e','%4.2f','\n');
% %   fprintf('rmse,'); prtmat_comma_screen(rmse_e','%4.2f','\n');
% %   fprintf('mae,'); prtmat_comma_screen(mae_e','%4.2f','\n');
% %   fprintf('min,'); prtmat_comma_screen(min_e','%4.2f','\n');
% %   fprintf('max,'); prtmat_comma_screen(max_e','%4.2f','\n');
% %   for i = 1:size(pct_vec,1);
% %     fprintf('%4.2f,',pct_vec(i)); prtmat_comma_screen(pct_e(:,i)','%4.2f','\n');
% %   end;
% %   
% %   % Compute some test statistics for differences in Forecasting performance
% %   e = [e_h(:,h) e_xe(:,h) e_xfe(:,h) e_3c(:,h) e_17c(:,h) e_l_xe(:,h) e_l_xfe(:,h) e_ao_h(:,h) e_ao_xe(:,h) e_ao_xfe(:,h)]; 
% %   e = e(i_fcst==1,:);
% %   e = packr(e);
% %   e2 = (e.^2);
% %   rslt_mat1 = NaN*zeros(4,4);
% %   rslt_mat2 = NaN*zeros(4,4);
% %   ikern = 0;  % flat
% %   nma = h;    % Number of MA terms
% %   for i = 1:nf;
% %       [ybar,vybar] = hacm(e2(:,i),nma,ikern);
% %       rslt_mat1(i,i) = ybar;
% %       rslt_mat2(i,i) = sqrt(vybar);
% %       for j = i+1:nf;
% %         [ybar,vybar] = hacm(e2(:,i)-e2(:,j),nma,ikern);
% %         rslt_mat1(i,j) = ybar;
% %         rslt_mat2(i,j) = sqrt(vybar); 
% %         rslt_mat1(j,i) = -ybar;
% %         rslt_mat2(j,i) = sqrt(vybar);
% %       end;
% %   end
% %   
% %   for i = 1:nf;
% %       for j = 1:nf;
% %           fprintf('%4.2f (%4.2f)',[rslt_mat1(i,j) rslt_mat2(i,j)]);
% %           if j < nf;
% %               fprintf(',');
% %           else;
% %               fprintf('\n');
% %           end;
% %       end;
% %   end;
% %     
% %  
% % -- Carry out filtering and so forth to determine mean and SD from various
% % tau(t/t) estimators
% 
% % % Exercise 1: Use system parameters from 17-component model
% % alpha_eps = alpha_eps_17c;
% % alpha_tau = alpha_tau_17c;
% % 
% % % sigma_eps_common = sigma_eps_common_total_17c;
% % % sigma_dtau_common = sigma_dtau_common_17c;
% % % sigma_eps_unique = sigma_eps_unique_total_17c;
% % % sigma_dtau_unique = sigma_dtau_unique_17c;
% % 
% % sigma_eps_common = sqrt(var_eps_common_total_17c);
% % sigma_dtau_common = sqrt(var_dtau_common_17c);
% % sigma_eps_unique = sqrt(var_eps_unique_total_17c);
% % sigma_dtau_unique = sqrt(var_dtau_unique_17c);
% % 
% % n_y = 17;
% % share_h = share_avg_60_end;
% % share_xe = share_h;
% % share_xe(:,[ienergy]) = 0;
% % share_xfe = share_h;
% % share_xfe(:,[ifood ienergy]) = 0;
% % share_xe = share_xe./repmat(sum(share_xe,2),1,n_y);
% % share_xfe = share_xfe./repmat(sum(share_xfe,2),1,n_y);
% % 
% % 
% % % Set up State Vector
% %   % --- State Vector
% %   %     (1) eps(t)
% %   %     (2) tau(t)
% %   %     (3) tau_u(t)
% %   %     (4) eps_u(t)
% %   ns = 2+2*n_y;      % size of state
% %   Q = zeros(ns,ns);
% %   F = zeros(ns,ns);
% %   F(2,2) = 1;
% %   F(3:2+n_y,3:2+n_y) = eye(n_y);
% %   H = zeros(ns,n_y);
% %   H(3:2+n_y,1:n_y) = eye(n_y);
% %   H(2+n_y+1:2+2*n_y,1:n_y) = eye(n_y);
% %   
% %   % Set up KF to run
% %   % Initial conditions
% %   X1_init = zeros(ns,1);
% %   P1_init = zeros(ns,ns);
% %   var_dtau_unique = big*ones(n_y,1);
% %   P1_init(3:2+n_y,3:2+n_y) = diag(var_dtau_unique);
% %   
% %   tau_17c_mean = zeros(dnobs_60_end,3);
% %   tau_17c_var = zeros(dnobs_60_end,3);
% %   X1 = X1_init;
% %   P1 = P1_init;
% %   for t = 1:dnobs_60_end;
% %     H_t = H;
% %     H_t(1,:) = alpha_eps(t,:);
% %     H_t(2,:) = alpha_tau(t,:);
% %     Q_t = diag([(sigma_eps_common(t,1).^2) (sigma_dtau_common(t,1).^2) (sigma_dtau_unique(t,:).^2) (sigma_eps_unique(t,:).^2)]');
% %     S_Y_DP = eye(n_y);
% %     yt = S_Y_DP*dp_disagg_60_end(t,:)';
% %     Rt = zeros(size(yt,1),size(yt,1));
% %     Ht = H_t*S_Y_DP';
% %     yt = dp_disagg_60_end(t,:)';
% %     Rt = zeros(size(yt,1),size(yt,1));
% %     Ht = H_t;
% %     [X1,P1,X2,P2,e,ht,llf] = kfilt_llf(yt,X1,P1,Ht,F,Rt,Q_t);
% %     % Parameters mapping state into tau-headline 
% %     S = zeros(ns,1);
% %     S(2) = share_h(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_h(t,:)';
% %     tauh = S'*X1;
% %     varh = S'*P1*S;
% %     
% %     S = zeros(ns,1);
% %     S(2) = share_xe(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_xe(t,:)';
% %     tauxe = S'*X1;
% %     varxe = S'*P1*S;
% %     
% %     S = zeros(ns,1);
% %     S(2) = share_xfe(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_xfe(t,:)';
% %     tauxfe = S'*X1;
% %     varxfe = S'*P1*S;
% %     
% %     tau_17c_mean(t,:) = [tauh tauxe tauxfe];
% %     tau_17c_var(t,:) = [varh varxe varxfe];
% %   
% %   end;
% %   
% %   tau_3c_mean = zeros(dnobs_60_end,3);
% %   tau_3c_var = zeros(dnobs_60_end,3);
% %   X1 = X1_init;
% %   P1 = P1_init;
% %   for t = 1:dnobs_60_end;
% %     H_t = H;
% %     H_t(1,:) = alpha_eps(t,:);
% %     H_t(2,:) = alpha_tau(t,:);
% %     Q_t = diag([(sigma_eps_common(t,1).^2) (sigma_dtau_common(t,1).^2) (sigma_dtau_unique(t,:).^2) (sigma_eps_unique(t,:).^2)]');
% %     % Mapping from dp_disagg to observable
% %     S_Y_DP = zeros(3,n_y);
% %     S_Y_DP(1,icore) = share_h(t,icore);
% %     S_Y_DP(2,ifood) = share_h(t,ifood);
% %     S_Y_DP(3,ienergy) = share_h(t,ienergy);
% %     yt = S_Y_DP*dp_disagg_60_end(t,:)';
% %     Rt = zeros(size(yt,1),size(yt,1));
% %     Ht = H_t*S_Y_DP';
% %     [X1,P1,X2,P2,e,ht,llf] = kfilt_llf(yt,X1,P1,Ht,F,Rt,Q_t);
% %     % Parameters mapping state into tau-headline 
% %     S = zeros(ns,1);
% %     S(2) = share_h(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_h(t,:)';
% %     tauh = S'*X1;
% %     varh = S'*P1*S;
% %     
% %     S = zeros(ns,1);
% %     S(2) = share_xe(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_xe(t,:)';
% %     tauxe = S'*X1;
% %     varxe = S'*P1*S;
% %     
% %     S = zeros(ns,1);
% %     S(2) = share_xfe(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_xfe(t,:)';
% %     tauxfe = S'*X1;
% %     varxfe = S'*P1*S;
% %     
% %     tau_3c_mean(t,:) = [tauh tauxe tauxfe];
% %     tau_3c_var(t,:) = [varh varxe varxfe];
% %   
% %   end;
% %   
% %   tau_h_mean = zeros(dnobs_60_end,3);
% %   tau_h_var = zeros(dnobs_60_end,3);
% %   X1 = X1_init;
% %   P1 = P1_init;
% %   for t = 1:dnobs_60_end;
% %     H_t = H;
% %     H_t(1,:) = alpha_eps(t,:);
% %     H_t(2,:) = alpha_tau(t,:);
% %     Q_t = diag([(sigma_eps_common(t,1).^2) (sigma_dtau_common(t,1).^2) (sigma_dtau_unique(t,:).^2) (sigma_eps_unique(t,:).^2)]');
% %     % Mapping from dp_disagg to observable
% %     S_Y_DP = zeros(1,n_y);
% %     S_Y_DP(1,:) = share_h(t,:);
% %     yt = S_Y_DP*dp_disagg_60_end(t,:)';
% %     Rt = zeros(size(yt,1),size(yt,1));
% %     Ht = H_t*S_Y_DP';
% %     [X1,P1,X2,P2,e,ht,llf] = kfilt_llf(yt,X1,P1,Ht,F,Rt,Q_t);
% %     % Parameters mapping state into tau-headline 
% %     S = zeros(ns,1);
% %     S(2) = share_h(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_h(t,:)';
% %     tauh = S'*X1;
% %     varh = S'*P1*S;
% %     
% %     S = zeros(ns,1);
% %     S(2) = share_xe(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_xe(t,:)';
% %     tauxe = S'*X1;
% %     varxe = S'*P1*S;
% %     
% %     S = zeros(ns,1);
% %     S(2) = share_xfe(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_xfe(t,:)';
% %     tauxfe = S'*X1;
% %     varxfe = S'*P1*S;
% %     
% %     tau_h_mean(t,:) = [tauh tauxe tauxfe];
% %     tau_h_var(t,:) = [varh varxe varxfe];
% %   
% %   end;
% %   
% %   tau_xe_mean = zeros(dnobs_60_end,3);
% %   tau_xe_var = zeros(dnobs_60_end,3);
% %   X1 = X1_init;
% %   P1 = P1_init;
% %   for t = 1:dnobs_60_end;
% %     H_t = H;
% %     H_t(1,:) = alpha_eps(t,:);
% %     H_t(2,:) = alpha_tau(t,:);
% %     Q_t = diag([(sigma_eps_common(t,1).^2) (sigma_dtau_common(t,1).^2) (sigma_dtau_unique(t,:).^2) (sigma_eps_unique(t,:).^2)]');
% %     % Mapping from dp_disagg to observable
% %     S_Y_DP = zeros(1,n_y);
% %     S_Y_DP(1,:) = share_xe(t,:);
% %     yt = S_Y_DP*dp_disagg_60_end(t,:)';
% %     Rt = zeros(size(yt,1),size(yt,1));
% %     Ht = H_t*S_Y_DP';
% %     [X1,P1,X2,P2,e,ht,llf] = kfilt_llf(yt,X1,P1,Ht,F,Rt,Q_t);
% %     % Parameters mapping state into tau-headline 
% %     S = zeros(ns,1);
% %     S(2) = share_h(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_h(t,:)';
% %     tauh = S'*X1;
% %     varh = S'*P1*S;
% %     
% %     S = zeros(ns,1);
% %     S(2) = share_xe(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_xe(t,:)';
% %     tauxe = S'*X1;
% %     varxe = S'*P1*S;
% %     
% %     S = zeros(ns,1);
% %     S(2) = share_xfe(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_xfe(t,:)';
% %     tauxfe = S'*X1;
% %     varxfe = S'*P1*S;
% %     
% %     tau_xe_mean(t,:) = [tauh tauxe tauxfe];
% %     tau_xe_var(t,:) = [varh varxe varxfe];
% %   
% %   end;
% %   
% %   tau_xfe_mean = zeros(dnobs_60_end,3);
% %   tau_xfe_var = zeros(dnobs_60_end,3);
% %   X1 = X1_init;
% %   P1 = P1_init;
% %   for t = 1:dnobs_60_end;
% %     H_t = H;
% %     H_t(1,:) = alpha_eps(t,:);
% %     H_t(2,:) = alpha_tau(t,:);
% %     Q_t = diag([(sigma_eps_common(t,1).^2) (sigma_dtau_common(t,1).^2) (sigma_dtau_unique(t,:).^2) (sigma_eps_unique(t,:).^2)]');
% %     % Mapping from dp_disagg to observable
% %     S_Y_DP = zeros(1,n_y);
% %     S_Y_DP(1,:) = share_xfe(t,:);
% %     yt = S_Y_DP*dp_disagg_60_end(t,:)';
% %     Rt = zeros(size(yt,1),size(yt,1));
% %     Ht = H_t*S_Y_DP';
% %     [X1,P1,X2,P2,e,ht,llf] = kfilt_llf(yt,X1,P1,Ht,F,Rt,Q_t);
% %     % Parameters mapping state into tau-headline 
% %     S = zeros(ns,1);
% %     S(2) = share_h(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_h(t,:)';
% %     tauh = S'*X1;
% %     varh = S'*P1*S;
% %     
% %     S = zeros(ns,1);
% %     S(2) = share_xe(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_xe(t,:)';
% %     tauxe = S'*X1;
% %     varxe = S'*P1*S;
% %     
% %     S = zeros(ns,1);
% %     S(2) = share_xfe(t,:)*alpha_tau(t,:)';
% %     S(3:2+n_y) = share_xfe(t,:)';
% %     tauxfe = S'*X1;
% %     varxfe = S'*P1*S;
% %     
% %     tau_xfe_mean(t,:) = [tauh tauxe tauxfe];
% %     tau_xfe_var(t,:) = [varh varxe varxfe];
% %   
% %   end;
% %   
% %   h_mean = NaN*zeros(dnobs_60_end,3);
% %   h_var = NaN*zeros(dnobs_60_end,3);
% %   xe_mean = NaN*zeros(dnobs_60_end,3);
% %   xe_var = NaN*zeros(dnobs_60_end,3);
% %   xfe_mean = NaN*zeros(dnobs_60_end,3);
% %   xfe_var = NaN*zeros(dnobs_60_end,3);
% %   for t = 1:dnobs_60_end;
% %     S_Y_DP = zeros(1,n_y);
% %     S_Y_DP = share_h(t,:);
% %     h_mean(t,1) = S_Y_DP*dp_disagg_60_end(t,:)';
% %     tmp1 = ((S_Y_DP*alpha_eps(t,:)')*sigma_eps_common(t,1))^2;
% %     tmp2 = sum((S_Y_DP.*sigma_eps_unique(t,:)).^2,2);
% %     h_var(t,1) = tmp1 + tmp2;
% %     
% %     S_Y_DP = zeros(1,n_y);
% %     S_Y_DP = share_xe(t,:);
% %     xe_mean(t,2) = S_Y_DP*dp_disagg_60_end(t,:)';
% %     tmp1 = ((S_Y_DP*alpha_eps(t,:)')*sigma_eps_common(t,1))^2;
% %     tmp2 = sum((S_Y_DP.*sigma_eps_unique(t,:)).^2,2);
% %     xe_var(t,2) = tmp1 + tmp2;
% %     
% %     S_Y_DP = zeros(1,n_y);
% %     S_Y_DP = share_xfe(t,:);
% %     xfe_mean(t,3) = S_Y_DP*dp_disagg_60_end(t,:)';
% %     tmp1 = ((S_Y_DP*alpha_eps(t,:)')*sigma_eps_common(t,1))^2;
% %     tmp2 = sum((S_Y_DP.*sigma_eps_unique(t,:)).^2,2);
% %     xfe_var(t,3) = tmp1 + tmp2;
% %     
% %   end;
% %   
% %   
% %   ao_h_mean = NaN*zeros(dnobs_60_end,3);
% %   ao_h_var = NaN*zeros(dnobs_60_end,3);
% %   ao_xe_mean = NaN*zeros(dnobs_60_end,3);
% %   ao_xe_var = NaN*zeros(dnobs_60_end,3);
% %   ao_xfe_mean = NaN*zeros(dnobs_60_end,3);
% %   ao_xfe_var = NaN*zeros(dnobs_60_end,3);
% %   for t = 1:dnobs_60_end;
% %     % Use approximation, locally constant parameter values 
% %     S_Y_DP = zeros(1,n_y);
% %     S_Y_DP = share_h(t,:);
% %     ao_h_mean(t,1) = S_Y_DP*dp4_disagg_60_end(t,:)';
% %     tmp1 = ((S_Y_DP*alpha_eps(t,:)')*sigma_eps_common(t,1))^2;
% %     tmp2 = sum((S_Y_DP.*sigma_eps_unique(t,:)).^2,2);
% %     var_eps_comp = (tmp1 + tmp2)/4;
% %     tmp1 = ((S_Y_DP*alpha_tau(t,:)')*sigma_dtau_common(t,1))^2;
% %     tmp2 = sum((S_Y_DP.*sigma_dtau_unique(t,:)).^2,2);
% %     var_tau_comp = (tmp1+tmp2)*14/16;
% %     ao_h_var(t,1) = var_eps_comp+var_tau_comp; 
% %     
% %     S_Y_DP = zeros(1,n_y);
% %     S_Y_DP = share_xe(t,:);
% %     ao_xe_mean(t,2) = S_Y_DP*dp4_disagg_60_end(t,:)';
% %     tmp1 = ((S_Y_DP*alpha_eps(t,:)')*sigma_eps_common(t,1))^2;
% %     tmp2 = sum((S_Y_DP.*sigma_eps_unique(t,:)).^2,2);
% %     var_eps_comp = (tmp1 + tmp2)/4;
% %     tmp1 = ((S_Y_DP*alpha_tau(t,:)')*sigma_dtau_common(t,1))^2;
% %     tmp2 = sum((S_Y_DP.*sigma_dtau_unique(t,:)).^2,2);
% %     var_tau_comp = (tmp1+tmp2)*14/16;
% %     ao_xe_var(t,2) = var_eps_comp+var_tau_comp; 
% %     
% %     S_Y_DP = zeros(1,n_y);
% %     S_Y_DP = share_xfe(t,:);
% %     ao_xfe_mean(t,3) = S_Y_DP*dp4_disagg_60_end(t,:)';
% %     tmp1 = ((S_Y_DP*alpha_eps(t,:)')*sigma_eps_common(t,1))^2;
% %     tmp2 = sum((S_Y_DP.*sigma_eps_unique(t,:)).^2,2);
% %     var_eps_comp = (tmp1 + tmp2)/4;
% %     tmp1 = ((S_Y_DP*alpha_tau(t,:)')*sigma_dtau_common(t,1))^2;
% %     tmp2 = sum((S_Y_DP.*sigma_dtau_unique(t,:)).^2,2);
% %     var_tau_comp = (tmp1+tmp2)*14/16;
% %     ao_xfe_var(t,3) = var_eps_comp+var_tau_comp; 
% %     
% %   end;
% %   
% %   tmp = sqrt(mean(tau_17c_var(i_fcst==1,:),1));
% %   prtmat_comma_screen(tmp,'%5.3f ','\n');
% %   tmp = sqrt(mean(tau_3c_var(i_fcst==1,:),1));
% %   prtmat_comma_screen(tmp,'%5.3f ','\n');
% %   tmp = sqrt(mean(tau_h_var(i_fcst==1,:),1));
% %   prtmat_comma_screen(tmp,'%5.3f ','\n');
% %   tmp = sqrt(mean(tau_xe_var(i_fcst==1,:),1));
% %   prtmat_comma_screen(tmp,'%5.3f ','\n');
% %   tmp = sqrt(mean(tau_xfe_var(i_fcst==1,:),1));
% %   prtmat_comma_screen(tmp,'%5.3f ','\n');
% %   tmp = sqrt(mean(h_var(i_fcst==1,:),1));
% %   prtmat_comma_screen(tmp,'%5.3f ','\n');
% %   tmp = sqrt(mean(xe_var(i_fcst==1,:),1));
% %   prtmat_comma_screen(tmp,'%5.3f ','\n');
% %   tmp = sqrt(mean(xfe_var(i_fcst==1,:),1));
% %   prtmat_comma_screen(tmp,'%5.3f ','\n');
% %   tmp = sqrt(mean(ao_h_var(i_fcst==1,:),1));
% %   prtmat_comma_screen(tmp,'%5.3f ','\n');
% %   tmp = sqrt(mean(ao_xe_var(i_fcst==1,:),1));
% %   prtmat_comma_screen(tmp,'%5.3f ','\n');
% %   tmp = sqrt(mean(ao_xfe_var(i_fcst==1,:),1));
% %   prtmat_comma_screen(tmp,'%5.3f ','\n');
%   
% %   hmax = 12;
% %   rw_fac = zeros(hmax,1);
% %   for h = 1:hmax;
% %       tmp = (1:1:h)'/h;
% %       rw_fac(h) = sum(tmp.^2);
% %   end;
% %   fcst_var=NaN*zeros(dnobs_60_end,hmax);
% %   % Compute Variance of Forecast error, known tau
% %   for t = 1:dnobs_60_end;
% %     for h =1:hmax;
% %     % Use approximation, locally constant parameter values 
% %       S_Y_DP = zeros(1,n_y);
% %       S_Y_DP = share_h(t,:);
% %       tmp1 = ((S_Y_DP*alpha_eps(t,:)')*sigma_eps_common(t,1))^2;
% %       tmp2 = sum((S_Y_DP.*sigma_eps_unique(t,:)).^2,2);
% %       var_eps_comp = (tmp1 + tmp2)/h;
% %       tmp1 = ((S_Y_DP*alpha_tau(t,:)')*sigma_dtau_common(t,1))^2;
% %       tmp2 = sum((S_Y_DP.*sigma_dtau_unique(t,:)).^2,2);
% %       var_tau_comp = (tmp1+tmp2)*rw_fac(h);
% %       fcst_var(t,h) = var_eps_comp+var_tau_comp;
% %     end;
% %   end;
% %   tmp = sqrt(mean(fcst_var(i_fcst==1,:),1))';
% %   tmp = [(1:1:hmax)' tmp];
% %   prtmat_comma_screen(tmp,'%5.3f ','\n');
% %   
%   
% % % Read in Percentiles from Posterion of filtered values 
% % ulabel = ['_PCEagg'];
% % str_tmp = [matdir 'tau_f_mean_pct' ulabel]; load(str_tmp);
% % tau_f_mean_pct_h = tau_f_mean_pct;
% %   
% % ulabel = ['_PCEagg_xfe'];
% % str_tmp = [matdir 'tau_f_mean_pct' ulabel]; load(str_tmp);
% % tau_f_mean_pct_xfe = tau_f_mean_pct;
% % 
% %  mlabel = '_mucsv_fl_tvp'; 
% %  str_tmp = [matdir 'agg_f_tau_mean_pct' mlabel]; load(str_tmp);
% %  tau_f_mean_pct_17c = agg_f_tau_mean_pct;
% %  
% % ulabel = ['_PCEagg'];
% % str_tmp = [matdir 'tau_f_mean_sd' ulabel]; load(str_tmp);
% % tau_f_mean_sd_h = tau_f_mean_sd;
% %   
% % ulabel = ['_PCEagg_xfe'];
% % str_tmp = [matdir 'tau_f_mean_sd' ulabel]; load(str_tmp);
% % tau_f_mean_sd_xfe = tau_f_mean_sd;
% % 
% %  mlabel = '_mucsv_fl_tvp'; 
% %  str_tmp = [matdir 'agg_f_tau_mean_sd' mlabel]; load(str_tmp);
% %  tau_f_mean_sd_17c = agg_f_tau_mean_sd;
% %  
% % %  width_h = tau_f_mean_pct_h(:,5)-tau_f_mean_pct_h(:,3);
% % %  width_xfe = tau_f_mean_pct_xfe(:,5)-tau_f_mean_pct_xfe(:,3);
% % %  width_17c = tau_f_mean_pct_17c(:,5)-tau_f_mean_pct_17c(:,3);
% % %  
% % %  sd_h = tau_f_mean_sd_h(:,2);
% % %  sd_xfe = tau_f_mean_sd_xfe(:,2);
% % %  sd_17c = tau_f_mean_sd_17c(:,2);
% % %  
% % %  plot(cal_60_end,width_17c,'- k','LineWidth',3);
% % %  hold on;
% % %   plot(cal_60_end,width_h,'- g','LineWidth',3);
% % %   plot(cal_60_end,width_xfe,'- r','LineWidth',3);
% % %  hold off;
% % %  
% % %  plot(cal_60_end,sd_17c,'- k','LineWidth',3);
% % %  hold on;
% % %   plot(cal_60_end,sd_h,'- g','LineWidth',3);
% % %   plot(cal_60_end,sd_xfe,'- r','LineWidth',3);
% % %  hold off;
% % 
% % plot(cal_60_end,dp_agg_60_end,'- k','LineWidth',2);
% % hold on;
% %   plot(cal_60_end,dp_agg_xfe_60_end,'- r','LineWidth',2);
% %   plot(cal_60_end,tau_f_mean_pct_17c(:,4),'-- b','LineWidth',3);
% %   plot(cal_60_end,tau_f_mean_pct_xfe(:,4),'-- r','LineWidth',3);
% % hold off;
% % xlim([2000 2015]);
% % legend('Headline Inflation','Core(XFE) Inflation','Headline Trend (17c)','Core Trend (Univariate)');
% % legend('Location','NorthEast');
% % ax = gca;
% % ax.FontSize = 30;
% % waitforbuttonpress;
% % 
% % %plot(cal_60_end,dp_agg_60_end,'- k','LineWidth',1);
% % %plot(cal_60_end,tau_f_mean_pct_17c(:,4),'- b','LineWidth',3);
% % plot(cal_60_end,tau_f_mean_pct_17c(:,2),': b','LineWidth',2);
% % hold on;
% %   plot(cal_60_end,tau_f_mean_pct_17c(:,2),': b','LineWidth',2);
% %   plot(cal_60_end,tau_f_mean_pct_17c(:,6),': b','LineWidth',2);
% %  % plot(cal_60_end,tau_f_mean_pct_xfe(:,4),'- r','LineWidth',3);
% %   plot(cal_60_end,tau_f_mean_pct_xfe(:,2),': r','LineWidth',2);
% %   plot(cal_60_end,tau_f_mean_pct_xfe(:,6),': r','LineWidth',2);
% % hold off;
% % xlim([2000 2015]);
% % ax = gca;
% % ax.FontSize = 30; 