% xxx (fix description using conference draft table number) Table 4 for Discussants' draft
% 
%

clear all;
small = 1.0e-10;
big = 1.0e+6;
this_date = datestr(now,'yyyymmdd');

% -----  File Directories -- overhead, etc.
 outdir = 'out/';
 figdir = 'fig/';
 matdir = 'mat/';
 p_initial = path;
 vardef_dir = 'variable_definitions/';            % Directory for necessary matlab_procs and functions                         
 p1=path(vardef_dir,p_initial);
 utility_dir = 'm_utilities';
 path(utility_dir,p1);

% -- Load Data
load_data=0;         % Set = 1 first time to read data from Excel files (then they are saved as MATLAB file)
fhsw_bpea_data;      % fhsw_bpea_data reads in the full dataset .. all variables, etc. saved in datain.xx

% Bandwidth for BW used for Okun's Law Cyclical Adjustement
bw_bw = 60;

% Cyclcical Adjustment Parameters
ismpl_ol = smpl(calvec,[1981 3],[2016 2],4);  % Sample period for estimating Cyclical adjustment coefficients
du_nlead = 2;                   % number of leads
du_nlag = 2;                    % number of lags
% Cyclical adjusted trend bandwidth
  bw_bw = 60;  % Biweight bandwidth for 
  
% -- Set up regressors for cylical adjustment
  d_urate = dif(cps_urate,1);     % first difference of unemployment rate
  xca = lag(d_urate,-du_nlead);   % xca is cyclical adjustment variable
  for i = (-du_nlead+1):du_nlag;
    xca = [xca lag(d_urate,i)];
  end;

% Output file
  outfile_name = [outdir ['Table_4_' this_date '.csv'] ];
  fileID = fopen(outfile_name,'w');

% Read in Factors and Construct Forecasts
nfactors = 4;
% Read in factors
bw_bw_factors = 100;
i_demean_factors = 1;
str_tmp = [matdir 'factor_all_' num2str(i_demean_factors) '_' num2str(bw_bw_factors) '_' num2str(nfactors) ];
load(str_tmp);

T_lambda_var_first = [1984 1];             % Sample period for VAR and lambda estiamtion
T_lambda_var_end = [2009 2];
T_oos_first = [2009 3];                    % First out-of-sample forecast date
T_oos_end = [2016 2];                      % Last out-of-sample forecast date
ismpl_lambda = smpl(calvec,T_lambda_var_first,T_lambda_var_end,4); % Smpl period of lambda estimation
smpl_par.nfirst = T_lambda_var_first;                      % VAR estimation Sample
smpl_par.nlast = T_lambda_var_end;
ismpl_oos = smpl(calvec,T_oos_first,T_oos_end,4);

% Construct VAR for Factors that will used for the forecasts
% Estimate VAR and save parameter values
var_par.nlag   = 4;    % number of lags
var_par.iconst = 0;    % include constant (no constant because estimated factors have zero mean)
var_par.icomp  = 1;    % compute companion form of model .. (1) excluding constant .. (2) include constant
smpl_par.calvec = calvec;
smpl_par.nper = 4;
varout = varest(factor_all,var_par,smpl_par);  % This estimated VAR and computes various matrices
% VAR Companion Matrix
M = varout.coef.M;
Q = varout.coef.Q;
G = varout.coef.G;
% Set up state vector for forecasting factors
factor_fcst = NaN(dnobs,nfactors);
iT_first = find(calvec > (-0.001+T_oos_first(1) + (T_oos_first(2)-1)/4),1);
iT_last = find(calvec > (-0.001+T_oos_end(1) + (T_oos_end(2)-1)/4),1);
fstate = NaN(nfactors*var_par.nlag,1);  % State Vector for Companion form of VAR
% Set up initial value of state vector for forecasts
for i = 0:var_par.nlag-1;
   j = 1+nfactors*i;
   fstate(j:j+nfactors-1) = factor_all(iT_first-1-i,:)';
end;
% Carry out forecasting ..
for t = iT_first:iT_last;
   fstate = M*fstate;
   factor_fcst(t,:)=(Q*fstate)';
end;
% Compute average value of factor_fcst over OOS period
factor_fcst_avg = mean(factor_fcst(ismpl_oos==1,:));

% Construct Variables used in Supply-side analyis;
supply_variables;  % This computes all of the data series used here

% Construct mu-trend forecasts
supply_trends;
     
% Variables to use for decomposition .. entered as pairs, series and trend
% forecasts
% Variable to use for decomposition
 
  ymat = [d_y_prod d_y_prod_mu_fcst];
  ymat = [ymat d_y_gdo d_y_gdo_mu_fcst];
  ymat = [ymat d_bus_y_gdo d_bus_y_gdo_mu_fcst];
  ymat = [ymat (d_y_prod-d_cps_pop) (d_y_prod_mu_fcst-d_cps_pop_mu_fcst)];
  ymat = [ymat (d_y_gdo-d_cps_pop) (d_y_gdo_mu_fcst-d_cps_pop_mu_fcst)];
  ymat = [ymat (d_bus_y_gdo-d_cps_pop) (d_bus_y_gdo_mu_fcst-d_cps_pop_mu_fcst)];
  ymat = [ymat d_bus_tfp d_bus_tfp_mu_fcst];
  ymat = [ymat kshare_d_bus_cap-kshare_d_cps_pop kshare_d_bus_cap_mu_fcst-kshare_d_cps_pop_mu_fcst];
  ymat = [ymat (lshare_d_bus_hours+lshare_d_bus_lq-lshare_d_cps_pop) (lshare_d_bus_hours_mu_fcst+lshare_d_bus_lq_mu_fcst-lshare_d_cps_pop_mu_fcst)];
  ymat = [ymat (d_bus_hours-d_cps_pop) (d_bus_hours_mu_fcst-d_cps_pop_mu_fcst)];
  ymat = [ymat (d_bus_hours-d_bus_emp) (d_bus_hours_mu_fcst-d_bus_emp_mu_fcst)];
  ymat = [ymat (d_bus_emp-d_cps_emp) (d_bus_emp_mu_fcst-d_cps_emp_mu_fcst)];
  ymat = [ymat (d_cps_emp-d_cps_lf) (d_cps_emp_mu_fcst-d_cps_lf_mu_fcst)];
  ymat = [ymat (d_cps_lf-d_cps_pop) (d_cps_lf_mu_fcst-d_cps_pop_mu_fcst)];
  ymat = [ymat (d_bus_y_gdo-d_bus_hours) (d_bus_y_gdo_mu_fcst-d_bus_hours_mu_fcst)];
  ymat = [ymat kshare2_d_bus_tfp kshare2_d_bus_tfp_mu_fcst];
  ymat = [ymat (kshare3_d_bus_cap-kshare3_d_bus_y_gdo) (kshare3_d_bus_cap_mu_fcst-kshare3_d_bus_y_gdo_mu_fcst)];
  ymat = [ymat d_bus_lq d_bus_lq_mu_fcst];
  ymat = [ymat (d_bus_cap-(d_bus_hours+d_bus_lq)) (d_bus_cap_mu_fcst-(d_bus_hours_mu_fcst+d_bus_lq_mu_fcst))];
  ymat = [ymat (d_bus_cap-d_bus_y_gdo) (d_bus_cap_mu_fcst-d_bus_y_gdo_mu_fcst)];
  %ymat = [ymat (d_bus_cap-d_cps_pop) (d_bus_cap_mu_fcst-d_cps_pop_mu_fcst)];
  %ymat = [ymat (d_cps_pop) (d_cps_pop_mu_fcst)];
  %ymat = [ymat (d_bus_cap) (d_bus_cap_mu_fcst)];

  % Labels for for output
  ymat_str = {...
  'GDP' ...
  'GDO' ...
  'GDO(Bus)' ...
  '\nGDP-Pop' ...
  'GDO-Pop' ...
  '\nGDO(bus)-Pop' ...
  'TFP' ...
  'kshare*(cap-pop)' ...
  'lshare(lq+hours-pop)' ...
  '\nHours-Pop' ...
  'Hours(bus)-Emp(bus)' ...
  'Emp(bus)-Emp(CPS)' ...
  'Emp(CPS)-LF(CPS)' ...
  'LF(CPS)-Pop(CPS)' ...
  '\nGDO(bus)-Hours' ...
  'TFP/lshare' ...
  '(kshare/lshare)*(Capital-Output)' ...
  'LQuality' ...
  '\n\nCapital-Hours-Quality' ...
  'Capital - output' ...
  %'Capital per capita' ...
  %'CPS population' ...
  %'Capital' ...
  }; 
 
  % File Setup
  % Print Results
   fprintf(fileID,',,Number of factors: %2i \n',nfactors);
   fprintf(fileID,',,VAR lags: %2i \n',var_par.nlag);
   fprintf(fileID,',,VAR and Lambda Sample period: %4i:Q%1i - %4i:Q%1i \n',[T_lambda_var_first T_lambda_var_end]);
   fprintf(fileID,',,OOS Sample period: %4i:Q%1i - %4i:Q%1i \n',[T_oos_first T_oos_end]);
   fprintf(fileID,'\n');
   fprintf(fileID,'Variable,,');
   fprintf(fileID,'Averages: Current Expansion \n');
   %fprintf(fileID,',,dy_mu_fcst,dy_cycl_fcst,dy_fcst_total_avg,dy_avg,shortfall,se_fac_fcst,,r2(for tab2) \n');
   fprintf(fileID,',,dy_fcst_total_avg,dy_avg,shortfall,se_fac_fcst,,r2(for tab2) \n');
  
  % Carry out Analysis
  for iseries = 1:size(ymat_str,2);
      sname = char(ymat_str(iseries));
      ii = 2*(iseries-1)+1;
      dy = ymat(:,ii);
      dy_mu_fcst = ymat(:,ii+1);
      fprintf(fileID,[sname ',,']);
      
      [dy_cycle,dy_ca,dy_mu,dy_irregular,var_dy_mu,b_okun_sum,se_b_okun_sum ] = cyclical_adjustment(dy,xca,ismpl_ol,bw_bw);
    
      % Detrended dy
      dy_dt = dy-dy_mu;
      
      % Estimate lambda over appropriate period
      [b_lambda,v_lambda,se_lambda,ser_lambda,rbarsq_lambda]=hac(dy_dt(ismpl_lambda==1),factor_all(ismpl_lambda==1,:),4,1);
      % Construct forecast
      dy_fac_fcst = factor_fcst*b_lambda;
      dy_fcst_error = dy-dy_mu_fcst-dy_fac_fcst;
      
      % Compute averages over expansion period
      dy_mean = mean(dy(ismpl_oos==1));
      dy_mu_fcst_mean = mean(dy_mu_fcst(ismpl_oos==1));
      dy_fac_fcst_mean = mean(dy_fac_fcst(ismpl_oos==1));
      se_dy_fac_avg_fcst = sqrt(factor_fcst_avg*v_lambda*factor_fcst_avg');
      dy_fcst_error_mean = mean(dy_fcst_error(ismpl_oos==1));
         
      % Print Results
      %fprintf(fileID,'%5.4f,',dy_mu_fcst_mean);
      %fprintf(fileID,'%5.4f,',dy_fac_fcst_mean);
      fprintf(fileID,'%5.4f,',dy_fac_fcst_mean+dy_mu_fcst_mean);
      fprintf(fileID,'%5.4f,',dy_mean);
      fprintf(fileID,'%5.4f,', -dy_fcst_error_mean);
      fprintf(fileID,'%5.4f,,', se_dy_fac_avg_fcst);
      fprintf(fileID,'%5.4f\n', rbarsq_lambda);
      
      % Save Series for plotting;
      if iseries == 1;
          dy_save = dy;
          dy_mu_fcst_save = dy_mu_fcst;
          dy_fac_fcst_save = dy_fac_fcst;
      else;
          dy_save = [dy_save dy];
          dy_mu_fcst_save = [dy_mu_fcst_save dy_mu_fcst];
          dy_fac_fcst_save = [dy_fac_fcst_save dy_fac_fcst];
      end;
  end;

% Add Figures
  % Figure 7
  %series_num = [5 6 10 15 11 12 13 14 7 18 19 20 21 22 23];
  series_num = [5 6 10 15 11 12 13 14 7 18 19 20];
  plot_label = {
      '(a) GDO per capita' ...
      '(b) Business GDO per capita' ...
      '(c) Labor hours per capita, business' ...
      '(d) Ouput per hour. business' ...
      '(e) Hours per worker, business' ...
      '(f) Ratio of business to CPS employment' ...
      '(g) CPS employment rate' ...
      '(h) Labor force participation rate' ...
      '(i) TFP' ...
      '(j) Labor quality' ...
      '(k) Capital-labor ratio' ...
      '(l) Capital-output ratio' ...
     % '(m) Capital per capita' ...
     % '(n) CPS population' ...
     % '(o) Capital' ...
  };

% Sample Period for Plots
smpl_plot = smpl(calvec,[2005 1],[2016 2],4);
axis_font_size = 20;
title_font_size = 18;
nplot_page = 6;  % Number of plots per page
for i = 1:length(series_num);
    slabel = char(plot_label(i));
    is = series_num(i);
    dy = dy_save(:,is);
    dy_mu_fcst = dy_mu_fcst_save(:,is);
    dy_fac_fcst = dy_fac_fcst_save(:,is);
    dy_fcst = dy_mu_fcst+dy_fac_fcst;
    figure('position',[200 200 550 400]);
    plot(calvec(smpl_plot==1),dy(smpl_plot==1),'- k','LineWidth',3);
    hold on;
       plot(calvec(smpl_plot==1),dy_mu_fcst(smpl_plot==1),'- b','LineWidth',2);
       plot(calvec(smpl_plot==1),dy_fcst(smpl_plot==1),': r','LineWidth',3);
    hold off;
    ax = gca;  % This gets axis limits for current plot (must be open)
    x_lim = ax.XLim;
    y_lim = ax.YLim;
    ax.FontSize = axis_font_size;
    title(slabel,'fontsize',title_font_size);
     if ((i-1)/nplot_page) - floor((i-1)/nplot_page) == 0;
     legend('actual','trend forecast','total forecast');
    end;
    recession_bars(x_lim,y_lim);
    hold on;
       plot(calvec(smpl_plot==1),dy(smpl_plot==1),'- k','LineWidth',3);
       plot(calvec(smpl_plot==1),dy_mu_fcst(smpl_plot==1),'- b','LineWidth',2);
       plot(calvec(smpl_plot==1),dy_fcst(smpl_plot==1),': r','LineWidth',3);
    hold off;
    xlim([2004 2017]);
    ylim(y_lim);
    savefig([figdir 'Fig_6_supply_' num2str(i) '_' this_date]);
    saveas(gcf,[figdir 'Fig_supply_' num2str(i) '.epsc']);
end;


path(p_initial);  % Reset path