% Compute and save factors
% MWW, 20170220
%

clear all;
small = 1.0e-10;
big = 1.0e+6;

% -----  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);

% ----------- Sample Period, Calendars and so forth for these data
[dnobs_m,calvec_m,calds_m] = calendar_make([1959 1],[2016 9],12);  % Monthly Calendar
[dnobs_q,calvec_q,calds_q] = calendar_make([1959 1],[2016 3],4);   % Quarterly Calendar
  
% All Parameters 
nfactors = 4;
fac_est_first = [1959 3];     % First period for factor + lambda estimation
fac_est_last = [2016 2];      % Truncation in 2016:Q1 because of 2 leads of du for detrending
    
% -- Load Data
  load_data=1;
  % Demeaning Parameters
  i_demean = 1;  % 0 Do Nothing
                 % 1 Eliminate low-frequency using local demeaning 
 
  bw_bw = 100;   % Bi-Weight Parameter for local demeaning

  datain_factor_bpea;      % datain_all reads in the full dataset .. all variables, etc. saved in datain.xx
  n_series = size(datain.bpdata,2);
  
  % Parameters for factor estimation program
  nfac.unobserved = nfactors;
  nfac.observed = 0;
  nfac.total = nfac.unobserved+nfac.observed;
  est_par.fac_par.nfac = nfac;
  est_par.fac_par.w = 1;                           % This contains data on the observed factors (scalar if none)
  est_par.smpl_par.nfirst = fac_est_first;         % start date
  est_par.smpl_par.nlast  = fac_est_last;    % end date for VAR and univariate AR 
  est_par.smpl_par.calvec = datain.calvec;         % calendar
  est_par.smpl_par.nper   = 4;                     % number of periods a year  
  est_par.fac_par.nt_min                  = 20;     % min number of obs for any series used to est factors
  est_par.lambda.nt_min                   = 40;     % min number of obs for any series used to estimate lamba, irfs, etc.
  est_par.fac_par.tol                     = 10^-10;  % precision of factor estimation (scaled by by n*t)  
  est_par.fac_par.lambda_constraints_est  = 1;  % no constraints on lambda
  est_par.fac_par.lambda_constraints_full = 1;  % no constraints on lambda
  
  % Estimate Factors and Lambda
  est_data = datain.bpdata(:,datain.bpinclcode==1);
  lsout = factor_estimation_ls(est_data, est_par);
  factor_estimate = lsout.fac;
  
  % Save Factors in Matlab Format for full sample
  % Calendars for problem
  nfirst = datain.calds(1,:);
  nlast = datain.calds(end,:);
  [dnobs,calvec,calds] = calendar_make([1947 1],[2016 4],4);  
  % -- Convert factor estimates to this calendar
  factor_all = NaN(dnobs,nfactors);
  ismpl = smpl(calvec,nfirst,nlast,4);
  factor_all(ismpl==1,:) = factor_estimate;    % Insample estimates of factors
  str_tmp = [matdir 'factor_all_' num2str(i_demean) '_' num2str(bw_bw) '_' num2str(nfactors) ];
  save(str_tmp,'factor_all'); 
    
path(p_initial);  % Reset path