% Table 3 and Figures 3, 4, and 5
% 
%

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

% --- Parameters for Cyclical Adjustment and Cyclically adjusted Trend
  % 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_3_' this_date '.csv'] ];
  fileID = fopen(outfile_name,'w');
  
 % Misc 
 fprintf(fileID,'Dates Used for Okuns law regressions:');
 cal = calds(ismpl_ol==1,:);
 fprintf(fileID,'%4i:Q%1i - %4i:Q%1i \n',[cal(1,1) cal(1,2) cal(end,1) cal(end,2)]);
 fprintf(fileID,'Number of leads: %2i \n', du_nlead);
 fprintf(fileID,'Number of lags: %2i \n', du_nlag);
 
 % Sample Period for current expansion
  ismpl_ce = smpl(calvec,[2009 3],[2016 2],4);  % Truncate in 2016:Q2 because 2 lead of unemployment rate.
  fprintf(fileID,'Current Expansion Dates Used:  ');
  cal = calds(ismpl_ce==1,:);
  fprintf(fileID,'%4i:Q%1i - %4i:Q%1i \n',[cal(1,1) cal(1,2) cal(end,1) cal(end,2)]);
  fprintf(fileID,'  Number of quarters: %2i \n', sum(ismpl_ce));
  
  % Previous 3 recessions
  % Sample perior indicators for previous three expansions
  ismpl_p1 = smpl(calvec,[1983 1],[1989 4],4);
  ismpl_p2 = smpl(calvec,[1991 2],[1998 1],4);
  ismpl_p3 = smpl(calvec,[2002 1],[2007 4],4);
  fprintf(fileID,'Previous 3 expansions  \n');
  cal = calds(ismpl_p1==1,:);
  fprintf(fileID,'%4i:Q%1i - %4i:Q%1i \n',[cal(1,1) cal(1,2) cal(end,1) cal(end,2)]);
  fprintf(fileID,'Number of quarters: %2i \n', sum(ismpl_p1));
  cal = calds(ismpl_p2==1,:);
  fprintf(fileID,'%4i:Q%1i - %4i:Q%1i \n',[cal(1,1) cal(1,2) cal(end,1) cal(end,2)]);
  fprintf(fileID,'Number of quarters: %2i \n', sum(ismpl_p2));
  cal = calds(ismpl_p3==1,:);
  fprintf(fileID,'%4i:Q%1i - %4i:Q%1i \n',[cal(1,1) cal(1,2) cal(end,1) cal(end,2)]);
  fprintf(fileID,'Number of quarters: %2i \n', sum(ismpl_p3));
  ismpl_pe = ismpl_p1+ismpl_p2+ismpl_p3;
  
 % Construct Variables used in Supply-side analyis;
 supply_variables;  % This computes all of the data series used here 
 
% Variable to use for decomposition
  ymat = d_y_prod;
  ymat = [ymat d_y_gdo];
  ymat = [ymat d_bus_y_gdo];
  ymat = [ymat d_y_prod-d_cps_pop];
  ymat = [ymat d_y_gdo-d_cps_pop];
  ymat = [ymat d_bus_y_gdo-d_cps_pop];
  ymat = [ymat d_bus_tfp];
  ymat = [ymat kshare_d_bus_cap-kshare_d_cps_pop];
  ymat = [ymat lshare_d_bus_lq+lshare_d_bus_hours-lshare_d_cps_pop];
  ymat = [ymat d_bus_hours-d_cps_pop];
  ymat = [ymat d_bus_hours-d_bus_emp];
  ymat = [ymat d_bus_emp-d_cps_emp];
  ymat = [ymat d_cps_emp-d_cps_lf];
  ymat = [ymat d_cps_lf-d_cps_pop];
  ymat = [ymat d_bus_y_gdo - d_bus_hours];
  ymat = [ymat kshare2_d_bus_tfp];
  ymat = [ymat kshare3_d_bus_cap-kshare3_d_bus_y_gdo]; 
  ymat = [ymat d_bus_lq];
  ymat = [ymat d_bus_cap-d_bus_hours-d_bus_lq];
  ymat = [ymat d_bus_cap-d_bus_y_gdo];
  ymat = [ymat d_bus_cap-d_cps_pop];
 
  % 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 - Pop' ...
  }; 

% File Setip 
  fprintf(fileID,'\n\n Variable,,');
  fprintf(fileID,'Historical Values (not cycl. adj.),,,,Cylcically adjusted \n');
  fprintf(fileID,',,3Previous,PostCrisisRecovery,Annual Shortfall,,3Previous,PostCrisisRecovery,Annual Shortfall,mu,z,Cumuulative percent \n');

  % Carry out Analysis
  for iseries = 1:size(ymat,2);
      sname = char(ymat_str(iseries));
      dz = ymat(:,iseries);
      fprintf(fileID,[sname ',,']);
      
      [dz_cycle,dz_ca,dz_mu,dz_irregular,var_dz_mu,b_okun_sum,se_b_okun_sum ] = cyclical_adjustment(dz,xca,ismpl_ol,bw_bw);
      
      % Actuals
      dz_mean_previous = mean(dz(ismpl_pe==1));
      dz_mean_current =  mean(dz(ismpl_ce==1));
      dz_ca_mean_previous = mean(dz_ca(ismpl_pe==1));
      dz_ca_mean_current =  mean(dz_ca(ismpl_ce==1));
      dz_shortfall = dz_mean_previous-dz_mean_current;
      dz_ca_shortfall = dz_ca_mean_previous-dz_ca_mean_current;
      dz_ca_cum = 100*(exp(28*dz_ca_shortfall/400)-1);
      dz_mu_shortfall =  mean(dz_mu(ismpl_pe==1)) -  mean(dz_mu(ismpl_ce==1));
      dz_irregular_shortfall =  mean(dz_irregular(ismpl_pe==1)) -  mean(dz_irregular(ismpl_ce==1));
      
      fprintf(fileID,'%5.4f,',dz_mean_previous);
      fprintf(fileID,'%5.4f,',dz_mean_current); 
      fprintf(fileID,'%5.4f,,',dz_shortfall);
      fprintf(fileID,'%5.4f,',dz_ca_mean_previous);
      fprintf(fileID,'%5.4f,',dz_ca_mean_current); 
      fprintf(fileID,'%5.4f,',dz_ca_shortfall);
      fprintf(fileID,'%5.4f,',dz_mu_shortfall);
      fprintf(fileID,'%5.4f,',dz_irregular_shortfall);
      fprintf(fileID,'%5.4f \n',dz_ca_cum);
      
      % Save Series for plotting;
      if iseries == 1;
          dz_save = dz;
          dz_ca_save = dz_ca;
          dz_mu_save = dz_mu;
      else;
          dz_save = [dz_save dz];
          dz_ca_save = [dz_ca_save dz_ca];
          dz_mu_save = [dz_mu_save dz_mu];
      end;
      
  end;  

  % Add a version of John's Figures
  % Figure 4
  series_num = [5 6 10 15];
  plot_label = {
      '(a) GDO per capita' ...
      '(b) Business GDO per capita' ...
      '(c) Labor hours per capita, business' ...
      '(d) Ouput per hour. business' ...
  };

axis_font_size = 20;
title_font_size = 18;
for i = 1:length(series_num);
    slabel = char(plot_label(i));
    is = series_num(i);
    dz = dz_save(:,is)/400;
    dz_ca = dz_ca_save(:,is)/400;
    dz_mu = dz_mu_save(:,is)/400;
    z = cum_nan(dz);
    z_ca = cum_nan(dz_ca);
    z_mu = cum_nan(dz_mu);
    % Normalize values to zero in 1981:Q3
    z = z - mean(z(ismpl_ol==1));
    z_ca = z_ca - mean(z_ca(ismpl_ol==1));
    z_mu = z_mu - mean(z_mu(ismpl_ol==1));
    TFirst = find(calvec == (1983 + (3-1)/4));
    znorm = z(TFirst);
    z = z - znorm;
    z_ca = z_ca - znorm;
    z_mu = z_mu - znorm;
    figure('position',[200 200 550 400]);
    plot(calvec(ismpl_ol==1),z(ismpl_ol==1),'- k','LineWidth',3);
    hold on;
       plot(calvec(ismpl_ol==1),z_ca(ismpl_ol==1),': r','LineWidth',3);
       plot(calvec(ismpl_ol==1),z_mu(ismpl_ol==1),'- b','LineWidth',2);
    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;
     legend('actual','cyclically adjusted','trend');
    end;
    recession_bars(x_lim,y_lim);
    hold on;
     plot(calvec(ismpl_ol==1),z(ismpl_ol==1),'- k','LineWidth',3);
     plot(calvec(ismpl_ol==1),z_ca(ismpl_ol==1),': r','LineWidth',3);
     plot(calvec(ismpl_ol==1),z_mu(ismpl_ol==1),'- b','LineWidth',2);
    hold off;
    xlim([1981 2017]);
    ylim(y_lim);
    savefig([figdir 'Fig_3_' num2str(i) '_' this_date]);
    saveas(gcf,[figdir 'Fig_3_' num2str(i) '.epsc']);
end;

% Figure 5
  series_num = [11 12 13 14];
  plot_label = {
      '(a) Hours per worker, business' ...
      '(b) Ratio of business to CPS employment' ...
      '(c) CPS employment rate' ...
      '(d) Labor force participation rate' ...
  };
for i = 1:length(series_num);
    slabel = char(plot_label(i));
    is = series_num(i);
    dz = dz_save(:,is)/400;
    dz_ca = dz_ca_save(:,is)/400;
    dz_mu = dz_mu_save(:,is)/400;
    z = cum_nan(dz);
    z_ca = cum_nan(dz_ca);
    z_mu = cum_nan(dz_mu);
    % Normalize values to zero in 1981:Q3
    z = z - mean(z(ismpl_ol==1));
    z_ca = z_ca - mean(z_ca(ismpl_ol==1));
    z_mu = z_mu - mean(z_mu(ismpl_ol==1));
    TFirst = find(calvec == (1983 + (3-1)/4));
    znorm = z(TFirst);
    z = z - znorm;
    z_ca = z_ca - znorm;
    z_mu = z_mu - znorm;
    %subplot(2,2,i);
    figure('position',[200 200 550 400]);
    plot(calvec(ismpl_ol==1),z(ismpl_ol==1),'- k','LineWidth',3);
    hold on;
       plot(calvec(ismpl_ol==1),z_ca(ismpl_ol==1),': r','LineWidth',3);
       plot(calvec(ismpl_ol==1),z_mu(ismpl_ol==1),'- b','LineWidth',2);
    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;
     legend('actual','cyclically adjusted','trend');
    end;
    recession_bars(x_lim,y_lim);
    hold on;
     plot(calvec(ismpl_ol==1),z(ismpl_ol==1),'- k','LineWidth',3);
     plot(calvec(ismpl_ol==1),z_ca(ismpl_ol==1),': r','LineWidth',3);
     plot(calvec(ismpl_ol==1),z_mu(ismpl_ol==1),'- b','LineWidth',2);
    hold off;
    xlim([1981 2017]);
    ylim(y_lim);
    savefig([figdir 'Fig_4_' num2str(i) '_' this_date]);
    saveas(gcf,[figdir 'Fig_4_' num2str(i) '.epsc']);
end;


% Figure 3
  series_num = [7 18 19 20 21];
  plot_label = {
      '(a) TFP' ...
      '(b) Labor quality' ...
      '(c) Capital-labor ratio' ...
      '(d) Capital-output ratio' ...
      '(e) Capital per capita' ...
  };
for i = 1:length(series_num);
    slabel = char(plot_label(i));
    is = series_num(i);
    dz = dz_save(:,is)/400;
    dz_ca = dz_ca_save(:,is)/400;
    dz_mu = dz_mu_save(:,is)/400;
    z = cum_nan(dz);
    z_ca = cum_nan(dz_ca);
    z_mu = cum_nan(dz_mu);
    % Normalize values to zero in 1981:Q3
    z = z - mean(z(ismpl_ol==1));
    z_ca = z_ca - mean(z_ca(ismpl_ol==1));
    z_mu = z_mu - mean(z_mu(ismpl_ol==1));
    TFirst = find(calvec == (1983 + (3-1)/4));
    znorm = z(TFirst);
    z = z - znorm;
    z_ca = z_ca - znorm;
    z_mu = z_mu - znorm;
    %subplot(2,2,i);
    figure('position',[200 200 550 400]);
    plot(calvec(ismpl_ol==1),z(ismpl_ol==1),'- k','LineWidth',3);
    hold on;
       plot(calvec(ismpl_ol==1),z_ca(ismpl_ol==1),': r','LineWidth',3);
       plot(calvec(ismpl_ol==1),z_mu(ismpl_ol==1),'- b','LineWidth',2);
    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;
     legend('actual','cyclically adjusted','trend');
    end;
    recession_bars(x_lim,y_lim);
    hold on;
     plot(calvec(ismpl_ol==1),z(ismpl_ol==1),'- k','LineWidth',3);
     plot(calvec(ismpl_ol==1),z_ca(ismpl_ol==1),': r','LineWidth',3);
     plot(calvec(ismpl_ol==1),z_mu(ismpl_ol==1),'- b','LineWidth',2);
    hold off;
    xlim([1981 2017]);
    ylim(y_lim);
    savefig([figdir 'Fig_5_' num2str(i) '_' this_date]);
    saveas(gcf,[figdir 'Fig_5_' num2str(i) '.epsc']);
end;
  
% ... overheard .. reset path to initial path  
path(p_initial);  % Reset path
