% Additional Panels for Figure 6
% 
%

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;


% 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;
d_y_prod_mu_fcst_supply = d_y_prod_mu_fcst;

% Demand Variables and trends
demand_variables;
demand_trends;
d_y_prod_mu_fcst = d_y_prod_mu_fcst_supply;
     
% Variables to use for decomposition .. entered as pairs, series and trend
% forecasts
 ylist = { ...
  'd_ces_emp' ...
  'd_bus_emp' ...
  'd_ces_emp_gov_sl' ...
  };
  
  % Labels for for output
  ymat_str = {...
  '(a) Employment - total' ...
  '(b) Employment - business'...
  '(c) Employment - state & local government' ...
  };
  
  % Carry out Analysis
 smpl_plot = smpl(calvec,[2005 1],[2016 2],4);
  axis_font_size = 20;
  title_font_size = 18;
  nplot_page = 4;  % Number of plots per page
  for iseries = 1:size(ylist,2);
      sname = char(ymat_str(iseries));
      sname_str = char(ylist(iseries));
      eval(['dy = ' sname_str ';']);
      eval(['dy_mu_fcst = ' sname_str '_mu_fcst;']);
     
      
      [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;
      dy_fcst = dy_fac_fcst+dy_mu_fcst;
      
      slabel = char(ymat_str(iseries));  
      
      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 ((iseries-1)/nplot_page) - floor((iseries-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_emp_' num2str(iseries) '_' this_date]);
      saveas(gcf,[figdir 'Fig_6_emp_' num2str(iseries) '.epsc']);
  end;
      
path(p_initial);  % Reset path