% Figure 1.b
% 
%

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;
  
  supply_variables;
  
  dz = d_bus_y_gdo - d_cps_pop;
  [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);
  dz = dz/400;
  dz_ca = dz_ca/400;
  z = cum_nan(dz);
  z_ca = cum_nan(dz_ca);
  % Normalize z values to zero in 1991:Q1
  TFirst = find(calvec == (1991 + (1-1)/4));
  znorm = z(TFirst);
  z = z - znorm;
  % Normalize z_ca to equal z in 2007:3
  iT = find(calvec == (2007 + (3-1)/4));
  znorm = z_ca(iT)-z(iT);
  z_ca = z_ca - znorm;
  
  axis_font_size = 20;
  title_font_size = 18;
  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);
  hold off;
    ax = gca;  % This gets axis limits for current plot (must be open)
    x_lim = ax.XLim;
    y_lim = ax.YLim;
    xlim([1989.999 2017]);
    ylim([-0.1 0.5]);
    ax.FontSize = axis_font_size;
    title('Business output per capita','fontsize',title_font_size);
    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);
    hold off;
    legend('actual','cyclically adjusted');
    savefig([figdir 'Fig_1b_' this_date]);
    saveas(gcf,[figdir 'Fig_1b.epsc']);

  
% ... overheard .. reset path to initial path  
path(p_initial);  % Reset path
