% Table 5 of Conference 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_5_' 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;
d_y_prod_mu_fcst_supply = d_y_prod_mu_fcst;

% Demand Variables and trends
demand_variables;
demand_trends;
d_y_prod_mu_fcst_demand = d_y_prod_mu_fcst;
dif_mu_gdp = d_y_prod_mu_fcst_supply-d_y_prod_mu_fcst_demand;
     
% Variables to use for decomposition .. entered as pairs, series and trend
% forecasts
 ylist = { ...
  'd_y_prod' ...
  'sd_pce' ...
  'sd_pce_goods' ...
  'sd_pce_goods_dur' ...
  'sd_pce_dg_mvp' ...
  'sd_pce_dg_fdhe' ...
  'sd_pce_dg_rgv' ...
  'sd_pce_dg_odg' ...
  'sd_pce_goods_nondur' ... 
  'sd_pce_ng_fb' ...
  'sd_pce_ng_cf' ...
  'sd_pce_ng_ge' ...
  'sd_pce_ng_ong' ...
  'sd_pce_services' ...
  'sd_pce_s_hu' ...
  'sd_pce_s_hc' ...
  'sd_pce_s_ts' ...
  'sd_pce_s_rs' ...
  'sd_pce_s_fsa' ...
  'sd_pce_s_fsi' ...
  'sd_pce_s_os' ...
  'sd_pce_s_npish' ...
  'sd_gpdi' ...
  'sd_fpi' ...
  'sd_fpi_nonres' ...
  'sd_gpdi_nonres_structures' ...
  'sd_gpdi_nonres_equip' ...
  'sd_gpdi_nonres_ipp' ...
  'sd_fpi_res' ...
  'sd_gpdi_res_structures' ...
  'sd_gpdi_res_equip' ...
  'sd_gov' ...
  'sd_gov_fed' ...
  'sd_gov_sl' ...
  'sd_exports' ...
  'sd_imports' ...
  'd_gov_cons_tp' ...
  };
  
  % Labels for for output
  ymat_str = {...
  'GDP' ...
  '\nshare pce'...
  'share pce goods' ...
  'share pce goods dur' ...
  'share pce dg mvp' ...
  'share pce dg fdhe' ...
  'share pce dg rgv' ...
  'share pce dg odg' ...
  'share pce goods nondur' ... 
  'share pce ng fb' ...
  'share pce ng cf' ...
  'share pce ng ge' ...
  'share pce ng ong' ...
  'share pce services' ...
  'share pce s hu' ...
  'share pce s hc' ...
  'share pce s ts' ...
  'share pce s rs' ...
  'share pce s fsa' ...
  'share pce s fsi' ...
  'share pce s os' ...
  'share pce s npish' ...
  '\nshare gpdi' ...
  'share fpi' ...
  'share fpi nonres' ...
  'share fpi nonres structures' ...
  'share fpi nonres equip' ...
  'share fpi nonres ipp' ...
  'share fpi res' ...
  'share fpi res structures' ...
  'share fpi res equip' ...
  '\nshare gov' ...
  'share gov fed' ...
  'share gov sl' ...
  '\nshare exports' ...
  'share imports' ...
  '\n\nAddendum\ngov cons and transfers'...
  };
  
  % 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_avg,share,Okun,OkunSE,3Previous,PostCrisisRecovery,Annual Shortfall,mu,z,dy_fcst_total_avg,dy_avg,shortfall,se_fac_fcst\n');
   
   % 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\n\n', sum(ismpl_p3));
   ismpl_pe = ismpl_p1+ismpl_p2+ismpl_p3;
   
  
  % Carry out Analysis
  for iseries = 1:size(ylist,2);
      sname = char(ymat_str(iseries));
      fprintf(fileID,[sname ',,']);
      sname_str = char(ylist(iseries));
      eval(['dy = ' sname_str ';']);
      eval(['dy_mu_fcst = ' sname_str '_mu_fcst;']);
      
      tmp = strfind(sname_str,'d_'); 
      share_str = sname_str(tmp(end)+2:end);
      eval(['share_dy = s_' share_str ';']);
      
      [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;
      
      % Actuals
      dy_ca_mean_previous = mean(dy_ca(ismpl_pe==1));
      dy_ca_mean_current =  mean(dy_ca(ismpl_ce==1));
      dy_ca_shortfall = dy_ca_mean_previous-dy_ca_mean_current;
      dy_mu_shortfall =  mean(dy_mu(ismpl_pe==1)) -  mean(dy_mu(ismpl_ce==1));
      dy_irregular_shortfall =  mean(dy_irregular(ismpl_pe==1)) -  mean(dy_irregular(ismpl_ce==1));
      
      
      % 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_cycle_mean = mean(dy_cycle(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');
      share_mean = mean(share_dy(ismpl_oos==1));
      
      % Adjustment for LFPR change in GDP trend
      dif_mu_gdp_mean = mean(dif_mu_gdp(ismpl_oos==1));
      if isnan(share_mean) == 0;
       dy_mu_fcst_mean = dy_mu_fcst_mean + share_mean*dif_mu_gdp_mean;
      end;
      dy_fcst_error_mean = dy_mean-dy_mu_fcst_mean-dy_fac_fcst_mean;
         
      % Print Results
      fprintf(fileID,'%5.4f,',dy_mean);
      fprintf(fileID,'%5.2f,',share_mean);
      fprintf(fileID,'%5.4f,',b_okun_sum/4);       % Divide by 4 Percentage points per quarter
      fprintf(fileID,'%5.4f,',se_b_okun_sum/4);
      fprintf(fileID,'%5.4f,',dy_ca_mean_previous);
      fprintf(fileID,'%5.4f,',dy_ca_mean_current); 
      fprintf(fileID,'%5.4f,',dy_ca_shortfall);
      fprintf(fileID,'%5.4f,',dy_mu_shortfall);
      fprintf(fileID,'%5.4f,',dy_irregular_shortfall);
      fprintf(fileID,'%5.4f,',dy_fac_fcst_mean+dy_mu_fcst_mean);
      fprintf(fileID,'%5.4f,', -dy_fcst_error_mean);
      fprintf(fileID,'%5.4f\n', se_dy_fac_avg_fcst);
      
      
  end;

  
path(p_initial);  % Reset path