%% Replication code for generating WEI and WEI using alternative specifications
%  Generates information for tables for Tables 2 and 7 and generates
%  Figures 1, 2, 4, 5 and 8

% LICENSE FOR CODE:
% Copyright Federal Reserve Bank of New York and Federal Reserve Bank of Dallas.
% You may reproduce, use, modify, make derivative works of, and distribute this code in whole or in part 
% so long as you keep this notice in the documentation associated with any distributed works. 
% Neither the names of the Federal Reserve Bank of New York and Federal Reserve Bank of Dallas nor the names 
% of any of the authors may be used to endorse or promote works derived from this 
% code without prior written permission. Portions of the code attributed to third 
% parties are subject to applicable third party licenses and rights. By your 
% use of this code you accept this license and any applicable third party license.
% THIS CODE IS PROVIDED ON AN "AS IS" BASIS, WITHOUT ANY WARRANTIES OR CONDITIONS 
% OF ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY WARRANTIES
% OR CONDITIONS OF TITLE, NON-INFRINGEMENT, MERCHANTABILITY OR FITNESS FOR A 
% PARTICULAR PURPOSE, EXCEPT TO THE EXTENT THAT THESE DISCLAIMERS ARE HELD TO 
% BE LEGALLY INVALID. THE FEDERAL RESERVE BANK OF NEW YORK AND THE FEDERAL 
% RESERVE BANK OF DALLAS ARE NOT, UNDER ANY CIRCUMSTANCES, LIABLE TO YOU FOR
% DAMAGES OF ANY KIND ARISING OUT OF OR IN CONNECTION WITH USE OF OR INABILITY
% TO USE THE CODE, INCLUDING, BUT NOT LIMITED TO DIRECT, INDIRECT, INCIDENTAL,
% CONSEQUENTIAL, PUNITIVE, SPECIAL OR EXEMPLARY DAMAGES, WHETHER BASED ON BREACH
% OF CONTRACT, BREACH OF WARRANTY, TORT OR OTHER LEGAL OR EQUITABLE THEORY, EVEN
% IF THE FEDERAL RESERVE BANK OF NEW YORK OR THE FEDERAL RESERVE BANK OF DALLAS
% HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES OR LOSS AND REGARDLESS
% OF WHETHER SUCH DAMAGES OR LOSS IS FORESEEABLE.
%% Preliminaries
clear; close all; clc;
addpath(genpath(fullfile('..','Functions'))); % Add necessary functions

%% Loading in Input Series and adding percent change (year over year) series
%  Note: Percent change series are alternative transformation for d52l transformation for main specification 

% Reading in data from Data Sheet in Data_WEI_Alt_Spec.csv
rawData      = readtable(fullfile('..','Data','Data_WEI_Alt_Spec.csv'));
rawData.Date = cellfun(@datenum, rawData.Date); % Converting dates to MATLAB datenums

% Reading in series information (series name, id, transformation)
seriesInfo = readtable(fullfile('..','Data','series_info.csv'));

% Adding series for percent change transformation
whichSeriesPC = {'LICN','PSTL', 'ASASI', 'RSTOTL', 'LIUN', 'RAS', 'PELOUS', 'FUEL'};
pcSeriesIndex = find(ismember(seriesInfo.SeriesID, whichSeriesPC));

for iPC = 1:length(whichSeriesPC) % Looping through IDs and adding copies of series to seriesInfo/rawData
   
    % Adding series to rawData
    colName = [whichSeriesPC{iPC}, '_PC'];
    rawData.(colName) = rawData.(whichSeriesPC{iPC});
    
    % Adding series infor to seriesInfo and changing transformation to 'pc'
    seriesInfo = [seriesInfo; seriesInfo(pcSeriesIndex(iPC), :)];
    seriesInfo.Transformation(end) = {'pc52'};
    seriesInfo.SeriesID{end}       = colName;
    
end

%% Transforming WEI input series

% Initializing table with transformed data w/ raw data
transfData = rawData;

% Looping through each input series in series info and making
% transformation
for iSeries = 1:size(seriesInfo, 1)
   
    % Pulling correct values from rawData
    seriesID     = seriesInfo.SeriesID{iSeries};
    seriesValues = rawData.(seriesID);
    
    % Making transformation
    transf                = seriesInfo.Transformation{iSeries};
    transfSeries          = transformData(seriesValues, transf);
    transfData.(seriesID) = transfSeries;
end

%% Setting Information for Each Specification

%  Specifications:
%  1) Baseline WEI
%  2) Baseline WEI rescaled using vintage of Real GDP as of 6/25/2020
%  3) WEI estimated using data from 1/1/2008 - today
%  4) WEI estimated using data from 1/1/2008 to 12/31/2014 
%  5) WEI estimated using data from 1/1/2015 - 2/2020
%  6) WEI + MBA
%  7) WEI no UI
%  8) WEI no UI 1/1/2008 - today
%  9) WEI + BCC
%  10) WEI public data (IC, CC, steel, fuel sales, unadjusted withholdings)
%  11) WEI-7 (not including tax withholdings, rail traffic, and
%  electricity)
%  12) WEI using sum of claims data instead of each claims individually
%  13) WEI using percent changes
%  14) WEI using ALS
%  15) WEI using DFM Filterted results (Note: Merged later in code, lines:
%  )


%  Which series in each specification?
nSpec    = 15; % Not including DFM specificaitons
cellStrs = cell(nSpec, 1);

cellStrs(1, 1)  = {{'Date','JRRSX', 'LICN','PSTL', 'ASASI', 'WITH', 'RSTOTL', 'LIUN', 'RAS', 'PELOUS', 'FUEL'}};
cellStrs(2, 1)  = {{'Date','JRRSX', 'LICN','PSTL','ASASI', 'WITH', 'RSTOTL', 'LIUN', 'RAS', 'PELOUS', 'FUEL'}};
cellStrs(3, 1)  = {{'Date','JRRSX', 'LICN','PSTL','ASASI', 'WITH', 'RSTOTL', 'LIUN', 'RAS', 'PELOUS', 'FUEL'}};
cellStrs(4, 1)  = {{'Date','JRRSX', 'LICN','PSTL','ASASI', 'WITH', 'RSTOTL', 'LIUN', 'RAS', 'PELOUS', 'FUEL'}};
cellStrs(5, 1)  = {{'Date','JRRSX', 'LICN','PSTL','ASASI', 'WITH', 'RSTOTL', 'LIUN', 'RAS', 'PELOUS', 'FUEL'}};
cellStrs(6, 1)  = {{'Date','JRRSX', 'LICN','PSTL','ASASI', 'WITH', 'RSTOTL', 'LIUN', 'RAS', 'PELOUS', 'FUEL', 'MBAMPN'}};
cellStrs(7, 1)  = {{'Date','JRRSX', 'PSTL', 'ASASI', 'WITH', 'RSTOTL', 'RAS', 'PELOUS', 'FUEL'}};
cellStrs(8, 1)  = {{'Date','JRRSX', 'PSTL', 'ASASI', 'WITH', 'RSTOTL', 'RAS', 'PELOUS', 'FUEL'}};
cellStrs(9, 1)  = {{'Date','JRRSX', 'LICN','PSTL','ASASI', 'WITH', 'RSTOTL', 'LIUN', 'RAS', 'PELOUS', 'FUEL', 'COMF'}};
cellStrs(10, 1) = {{'Date', 'LICN','PSTL', 'TDW', 'LIUN', 'FUEL'}};
cellStrs(11, 1) = {{'Date', 'JRRSX','LICN','PSTL', 'ASASI', 'LIUN', 'RAS', 'FUEL'}};
cellStrs(12, 1) = {{'Date','JRRSX','PSTL' 'ASASI', 'WITH', 'RSTOTL', 'CLAIMS', 'RAS', 'PELOUS', 'FUEL'}};
cellStrs(13, 1) = {{'Date','JRRSX', 'LICN_PC','PSTL_PC', 'ASASI_PC', 'WITH', 'RSTOTL_PC', 'LIUN_PC', 'RAS_PC', 'PELOUS_PC', 'FUEL_PC'}};
cellStrs(14, 1) = {{'Date','JRRSX', 'LICN','PSTL', 'ASASI', 'WITH', 'RSTOTL', 'LIUN', 'RAS', 'PELOUS', 'FUEL'}};
cellStrs(15, 1)  = {{'Date','JRRSX', 'LICN','PSTL', 'ASASI', 'WITH', 'RSTOTL', 'LIUN', 'RAS', 'PELOUS', 'FUEL'}};


%  Sample Estimation Start and End?
startDates    = repmat(datenum(2008, 1, 1), nSpec, 1);
startDates(5) = datenum(2015, 1, 1); 

endDates    = repmat(datenum(2020, 2, 29), nSpec, 1);
endDates(3) = datenum(2020, 12, 27); % WEI Baseline 2008 - latest
endDates(4) = datenum(2014, 12, 31); % WEI 2008 - 2014
endDates(8) = datenum(2020, 12, 27); % WEI No UI 2008 - latest

%  Which vintage of GDP?
whichVintageGDP    = repmat(datenum(2020, 12, 22), nSpec, 1); 
whichVintageGDP(2) = datenum(2020, 6, 25); % Before Annual Benchmark revisions

% AlS version?
isALS     = zeros(nSpec, 1);
isALS(14) = 1; 

%% Constructing WEI
%  Looping through each specifcation and constructing the WEI

% WEI start and end dates
weiStart = datenum(2007, 12, 30);
weiEnd   = datenum(2020, 12, 27);

% Storage for series/info each specifcation
WEICell      = cell(nSpec, 1); % Each alt. spec.
weightCell   = cell(nSpec, 1); % weights recovered from alt. spec.
varExplMat   = NaN(nSpec, 1);  % variance explained from each alt. spec.

% Setting and storing seed
for iSpec = 1:(nSpec - 1) % Loop through each specification thats not a DFM specification

    % Series  for the specification of interest
    seriesIDs    = cellStrs{iSpec};
    seriesSubset = transfData(:, seriesIDs);

    % WEI Estimation start and end dates
    estStart = startDates(iSpec);
    estEnd   = endDates(iSpec);
    
    % Subsetting data based on estimation start and end dates
    estKeepIndex = seriesSubset.Date >= estStart & seriesSubset.Date <= estEnd;
    XEst         = seriesSubset{estKeepIndex, 2:end}; 

    % Are we using ALS to in estimation?
    if isALS(iSpec) == 0


        % Estimating weights for WEI (Note: weights not identified to sign)
        meanEst      = nanmean(XEst, 1);
        stdEst       = nanstd(XEst, 0, 1);
        XNormEst     = (XEst - repmat(meanEst, size(XEst, 1), 1)) ./ repmat(stdEst, size(XEst, 1), 1); % Normalizing estimation data
        [loadings, ~, ~, ~, explained] = pca(XNormEst); % Using PCA to find loadings of series (factors -> series)
 

    else
        meanEst      = nanmean(XEst, 1);
        stdEst       = nanstd(XEst, 0, 1);
        XEst(isoutlier(XEst - lagmatrix(XEst, 1), 1)) = NaN; % Removing outliers based on change in value of series
        XNormEst     = (XEst - repmat(meanEst, size(XEst, 1), 1)) ./ repmat(stdEst, size(XEst, 1), 1); % Normalizing estimation data    
        rng(1, 'twister'); % Setting seed to produce identical results, since ALS initialized with random starting values
        [loadings, ~, ~, ~, explained] = pca(XNormEst, 'Algorithm', 'als'); % Using PCA to find loadings of series (factors -> series), als algorithm for missing data

    end
    
    % Storing the main specification, input estimation data
    if iSpec == 1
        mainSpecXNormEst = XNormEst; 
    end
    
    % Constructing weights
    weights      = loadings ^ -1; % Taking inverse to find weights of series (series -> factors)
    weights      = weights(1, :)'; % Only interested in weights of first principal component

    % Generating WEI (unscaled) for full sample period
    fullKeepIndex = seriesSubset.Date >= weiStart & seriesSubset.Date <= weiEnd;
    XFull         = seriesSubset{fullKeepIndex, 2:end}; % Subsetting data
    datesFull     = seriesSubset.Date(fullKeepIndex);   % Subsetting dates for subsequent plotting
    XNormFull     = (XFull - repmat(meanEst, size(XFull, 1), 1)) ./ repmat(stdEst, size(XFull, 1), 1); % Normalizing full sample period data based on estimation means and stds
    weiFactor     = XNormFull * weights;      
    
    % Assignign values for storage
    WEICell{iSpec,1}      = weiFactor;
    weightCell{iSpec, 1}  = weights;
    varExplMat(iSpec)     = explained(1);
end    
    
%% Merging DFM results

% Loading in and merging 
dfmWEItable               = readtable(['..' filesep 'Data' filesep 'Outdata' filesep 'wei_dfm.csv']);
WEICell{iSpec + 1, 1}     = dfmWEItable{:,2};

% Constructing weights based on loadings 
dfmOutput = readcell(fullfile('..','Data', 'Outdata', 'dfm_model_results.csv'));
loadings  = dfmOutput(7:3:34, 3);
loadings  = cellfun(@(x) str2double(x(1:end-4)), loadings);

weights     = loadings; 
signWeights = sign(weights); % Are they negative or positive?
weights     = abs(weights) ./ sum(abs(weights)); % Normalizing weights based on magnitude
weights     = weights .* sum(abs(weightCell{1,1})); % Adjusting to same scale as baseline weights 
weights     = weights .* signWeights;
weightCell{iSpec + 1, 1} = weights;
% Constructing explained variance (over estimation period)
keepIndex      =  find(datesFull >= startDates(end) & datesFull <= endDates(end)); % Same set of dates from PCA version and DFM version
dfmSeries      =  WEICell{iSpec + 1, 1};
dfmSeries      = dfmSeries(keepIndex);

fitValues = dfmSeries * loadings';
residual  = mainSpecXNormEst - fitValues;
residVar  = sum(sum(residual .^ 2, 1)); % Sum of Squared Residuals
actVar    = sum(sum(mainSpecXNormEst .^ 2, 1)); % Variance of actual adata

varExplMat(end) = (1 - (residVar ./ actVar)) .* 100; % How much variance is explained?
    
%% Rescaling each specification in GDP (YoY% Terms)

% Storage for rescaled WEI series
WEIGDPMat = NaN(size(WEICell{iSpec + 1, 1}, 1), nSpec);


% Loading in vintages of GDP/transforming to d4l
gdpVintages = readtable(['..' filesep 'Data' filesep 'gdp_vintage_data_full.csv']); 
nGDPVint    = size(gdpVintages, 2) - 1;
for iGDP = 2:size(gdpVintages, 2) % Looping through vintage to construct d4l version
    
    gdp     = gdpVintages{:,iGDP};
    gdp_d4l = [NaN(4,1); (log(gdp(5:end)) - log(gdp(1:end-4))) * 100];
    
    sName = gdpVintages.Properties.VariableNames{iGDP};
    sName = [sName '_d4l'];
    
    gdpVintages.(sName) = gdp_d4l;
    
    
end

gdpVintages.DATE = datenum(gdpVintages.DATE);

% Looping through and rescaling
for i = 1:nSpec
    
    
    % When does the reestimation period start and end?
    startEstim = startDates(i);
    endEstim   = endDates(i);
    
    % Construct a set of quarters    
    startQ = datenum(year(startEstim), (quarter(startEstim) - 1) * 3 + 1, 1); 
    endQ   = datenum(year(endEstim), (quarter(endEstim) - 1) * 3 + 1, 1);
    qList  = datenum(datetime(startQ, 'ConvertFrom', 'datenum'):calmonths(3):datetime(endQ, 'ConvertFrom', 'datenum'));
    qList  = eomdate(year(qList), month(qList) + 2); % End date of quarter
    
    % Finding which quarters to include (start and end) in estimation for rescaling
    rescaleStart = qList(find(qList >= startEstim, 1, 'first'));
    rescaleEnd   = qList(find(qList <= endEstim, 1, 'last'));  
    
    rescaleStart = datenum(year(rescaleStart), month(rescaleStart) - 2, 1); % Adjusting dates back to start of quarter
    rescaleEnd   = datenum(year(rescaleEnd),   month(rescaleEnd) - 2, 1);

    % Selecting correct specification
    weiFactor = WEICell{i, 1};
    
    % Converting unscaled WEI to a quarterly average
    weiTT = array2timetable(weiFactor, 'RowTimes', datetime(datesFull, 'ConvertFrom', 'datenum'));
    weiTT = retime(weiTT, 'daily','previous'); 
    weiTT = retime(weiTT, 'quarterly', 'mean');

    % Subsetting quarterly WEI to period used for rescaling
    weiRescale = weiTT.weiFactor(weiTT.Time >= datetime(rescaleStart, 'ConvertFrom', 'datenum') & ...
                                 weiTT.Time <= datetime(rescaleEnd, 'ConvertFrom', 'datenum'));

    
    % Which values of GDP will we use to rescale WEI?
    vintName      = ['GDPC1_', datestr(whichVintageGDP(i), 'yyyymmdd'), '_d4l'];
    d4lGDPFull    = gdpVintages.(vintName);                            
    whichGDPIndex = gdpVintages.DATE >= rescaleStart & gdpVintages.DATE <= rescaleEnd;
    gdpRescale    = d4lGDPFull(whichGDPIndex);

    % Generating rescaling coefficients (regressing quarterly WEI on GDP)
    rescaleCoeff = regress(gdpRescale, [ones(size(weiRescale)), weiRescale]);
    
    % Rescaling WEI
    WEI             = [ones(size(weiFactor)), weiFactor] * rescaleCoeff;
    WEIGDPMat(:, i) = WEI;


end

%% Aggregating weights into a single matrix

% Which series should have a negative sign for interpretability?
signFlip = {'LICN', 'LIUN', 'CLAIMS', 'LICN_PC', 'LIUN_PC'};

% Storage for weights
weightMat     = NaN(size(seriesInfo, 1), nSpec);
seriesIDs     = seriesInfo.SeriesID; % ID and Names needed for matching
seriesNames   = seriesInfo.SeriesName;

for i = 1:nSpec % Looping through specification and assinging weights to correct rows of weightMat
    
    whichIDs = cellStrs{i,1};
    whichIDs = whichIDs(2:end); % Which series including in specification?
    
    weights  = weightCell{i,1}; % Weights for a given specification
    
    % Checking weights for interpretation
    claimIndex = find(ismember(whichIDs, signFlip));
    if ~isempty(claimIndex) % Is there a claims series?
        
       claimIndex = claimIndex(1);
       if weights(claimIndex) > 0 % If its greater than 0 flip the signs for interpretation
          weights = weights * -1; 
       end
       
    elseif weights(1) < 0 % If no claims series and weights < 0, flip all the signs fo the weights for interpretation
       
       weights = weights * -1;  
        
    end
    
    
    % Including pc and d52l transformation in the same row to make easier
    % to read
    [~, wIndex] = ismember(whichIDs, seriesIDs);
    wNames      = seriesNames(wIndex);
    [~, wIndex] = ismember(wNames, seriesNames); % ismember finds the first instance
    weightMat(wIndex, i) = weights;
end

% Dropping percent change variable names and weights b/c they are in same row as d52l
keep        = sum(~isnan(weightMat), 2) > 0;
weightMat   = weightMat(keep, :);
seriesNames = seriesNames(keep);

%% Correlation matrix

corrMat  = corr(WEIGDPMat);

%% Writing results to csv

% Column Names for Each Specification
specIDs = {'WEI', 'WEI_GDP_0625', 'WEI_full', 'WEI_pre_2015', 'WEI_post_2015',...
                    'WEI_MBA', 'WEI_NoUI', 'WEI_NoUI_full', 'WEI_BCC', 'WEI_public', ...
                    'WEI_7', 'WEI_sum_claims', 'WEI_pc', 'WEI_ALS', 'WEI_dfm_filter'};
                
% Table for full history of each series
dates             = datesFull + 6; % Adjusting to Saturday
weiSpecMat        = [dates, WEIGDPMat];
weiSpecTable      = array2table(weiSpecMat, 'VariableNames', [{'Date'}, specIDs]);
weiSpecTable.Date = datestr(weiSpecTable.Date, 'mm/dd/yyyy');
% Table for Weights
weightTable = array2table([num2cell(weightMat), seriesNames], 'VariableNames', [specIDs, {'Series'}]);

% Table for Correlations
colNames  = [{'Index'}, specIDs];
rowNames  = colNames(2:end)';
corrTable = array2table([rowNames, num2cell(corrMat)], 'VariableNames', colNames);

% Table for Variance Explained
varTable = array2table(varExplMat', 'VariableNames', specIDs);
    
% Writing Results
writetable(weiSpecTable, fullfile('..','Data','Outdata','wei_alt_spec_series.csv'));
writetable(weightTable,  fullfile('..','Data','Outdata','wei_alt_spec_weights.csv'));
writetable(corrTable,    fullfile('..','Data','Outdata','wei_alt_spec_corr.csv'));
writetable(varTable,     fullfile('..','Data','Outdata','wei_alt_spec_explvar.csv'));


%% Creating Figures

% Path for figures
figPath = fullfile('..','Figures');


%%% Generating Figure 1 - Plotting Baseline WEI and WEI from Jan-2020 to
%%% Jan-2021
WEI       = weiSpecTable.WEI; % Baseline
plotDates = datenum(weiSpecTable.Date);

Fig = figure; % Overall Figure
plot(plotDates, WEI, 'LineWidth', 2);
grid on; box on;
ylim([-12, 6]);
xlim([plotDates(1) - 28, plotDates(end) + 56]);
xticks(datenum(datetime(2008, 1, 1):calyears(2):datetime(2020, 1, 1)))
datetick('x','yyyy','keeplimits', 'keepticks');
ylabel('Percent');
set(gca, 'FontSize', 16);
Fig.Units = 'Inches';
Fig.OuterPosition = [0 0 12 6];
tightfig;
print(Fig, fullfile(figPath, 'Figure1_Panel1.png'), '-dpng', '-r300');
print(Fig, fullfile(figPath, 'Figure1_Panel1.pdf'), '-dpdf', '-r300');

keepIndex = plotDates >= datenum(2020, 1, 1); % Subset for WEI after start of 2020
Fig = figure;
plot(plotDates(keepIndex), WEI(keepIndex), 'LineWidth', 2);
grid on; box on;
ylim([-12, 4]);
xlim([datenum(2020, 1, 1) - 14, datenum(2021, 1, 2) + 14])
xticks(datenum(datetime(2008, 1, 1):calyears(2):datetime(2020, 1, 1)))
datetick('x','yyyy','keeplimits', 'keepticks');
ylabel('Percent');
set(gca, 'FontSize', 16);
Fig.Units = 'Inches';
Fig.OuterPosition = [0 0 12 6];
tightfig;
print(Fig, fullfile(figPath, 'Figure1_Panel2.png'), '-dpng', '-r300');
print(Fig, fullfile(figPath, 'Figure1_Panel2.pdf'), '-dpdf', '-r300');

%%% Generating Figure 2 - Compare Baseline WEI against constituent series
consSeries = transfData(:, cellStrs{1, 1}); % Input series used by WEI
keepIndex  = consSeries.Date >= datenum(2007, 12, 30); % Subsetting for post 2008
consSeries = consSeries(keepIndex, :);
consSeries = consSeries{:, 2:end}; % Turning table output into a matrix

orderIndex   = [1 8 2 7 4 5 3 10 6 9]; % Re-ordering series to match order of plots
figureTitles = {'Redbook Retail' ,'Rasmussen Consumer Index, 52-week % change',...
                'Initial Claims, -52 week % change', 'Continuing Claims, -52 % change', ...
                'Staffing Index, 52-week % change', 'Tax Withholding',...
                'Steel Production, 52-week % change', 'Fuel Sales, 52-week % change',...
                'Rail Traffic, 52-week % change', 'Electricity Usage, 52-week % change'};
            
% Computing mean and standard deviation of WEI from start - 2/2020
weiMean = nanmean(WEI(plotDates >= datenum(2008, 1, 1) & plotDates <= datenum(2020, 2, 29)));
weiStd  = nanstd(WEI(plotDates >= datenum(2008, 1, 1) & plotDates <= datenum(2020, 2, 29)));
            

ylimMat        = repmat([-15, 10], length(orderIndex), 1); % y-axis limits for figures
ylimMat(3, :)  = [-30, 5];  % adjust Initial Claims
ylimMat(4, :)  = [-20, 5];  % adjust Cont. Claims
ylimMat(8, :)  = [-20, 10]; % adjust Fuel Sales

flipSign  = zeros(length(orderIndex)); % Flipping sign for IC and CC
flipSign(3) = 1;
flipSign(4) = 1;

for iFig = 1:length(orderIndex)

    % Picking correct series and restandardizing it in baseline WEI terms
    orderI      = orderIndex(iFig);
    compWEI     = consSeries(:, orderI);
    if flipSign(iFig) == 1
       compWEI = compWEI .* -1; 
    end
    compWEIMean = nanmean(compWEI(plotDates >= datenum(2008, 1, 1) & plotDates <= datenum(2020, 2, 29)));
    compWEIStd  = nanstd(compWEI(plotDates >= datenum(2008, 1, 1) & plotDates <= datenum(2020, 2, 29)));
    compWEI     = (compWEI - compWEIMean) ./ compWEIStd; % Normalize
    compWEI     = (compWEI .* weiStd) + weiMean; % Rescale in terms of baseline WEI
    
    
    % Generating Figure
    Fig = figure;
    hold on;
    plot(plotDates, WEI, 'LineWidth', 2);
    plot(plotDates, compWEI, 'LineWidth', 2);
    grid on; box on;
    title(figureTitles{iFig})
    ylim(ylimMat(iFig, :));
    xlim([plotDates(1) - 28, plotDates(end) + 56]);
    xticks(datenum(datetime(2008, 1, 1):calyears(2):datetime(2020, 1, 1)))
    datetick('x','yyyy','keeplimits', 'keepticks');
    ylabel('Percent');
    if iFig == 2 % Only produce legend for DFM, filtered estimates
       legend({'WEI', 'Weekly Indicator'}, 'Location', 'southwest', 'Box', 'off') 
    end
    set(gca, 'FontSize', 16);
    Fig.Units = 'Inches';
    Fig.OuterPosition = [0 0 9 6];
    tightfig;
    print(Fig, fullfile(figPath, ['Figure2_Panel', num2str(iFig), '.png']), '-dpng', '-r300');
    print(Fig, fullfile(figPath, ['Figure2_Panel', num2str(iFig), '.pdf']), '-dpdf', '-r300');
    

end

%%% Generating Figure 4 - Comparing Baseline WEI against other
%%% specifications

orderIndex   = [14, 15, 4, 5, 3, 11, 10, 6, 9, 7, 12, 13]; % Order of Transformation to Comapre
figureTitles = {'Trimmed (ALS)', 'DFM, filtered estimates',...
                '1/2008-12/2014', '1/2015-2/2020', '1/2008-12/2020',...
                'WEI-7', 'Public Data Only','MBA Mortgages', ...
                'Bloomberg Consumer Comfort', 'No UI claims',...
                'Total Claims', 'Percent Change Transformation'};
            
ylimMat        = repmat([-12, 6], length(orderIndex), 1); % y-axis limits for figures
ylimMat(1, :)  = [-20, 5];  % adjust for ALS
ylimMat(2, :)  = [-20, 5];  % adjust for DFM 
ylimMat(6, :)  = [-15, 5];  % adjust for WEI-7
ylimMat(7, :)  = [-20, 5];  % adjust for WEI public data only
ylimMat(12, :) = [-50, 10]; % adjust for WEI percent change transformation


% Looping through and generating figures
for iFig = 1:length(orderIndex)
   
    
    % Picking correct series and restandardizing it in baseline WEI terms
    orderI      = orderIndex(iFig);
    compWEI     = table2array(weiSpecTable(:, orderI + 1));
    compWEIMean = nanmean(compWEI(plotDates >= datenum(2008, 1, 1) & plotDates <= datenum(2020, 2, 29)));
    compWEIStd  = nanstd(compWEI(plotDates >= datenum(2008, 1, 1) & plotDates <= datenum(2020, 2, 29)), 0, 1);
    compWEI     = (compWEI - compWEIMean) ./ compWEIStd; % Normalize
    compWEI     = (compWEI .* weiStd) + weiMean; % Rescale in terms of baseline WEI
    
    % Generating Figure
    Fig = figure;
    hold on;
    plot(plotDates, WEI, 'LineWidth', 2);
    plot(plotDates, compWEI, 'LineWidth', 1.5);
    grid on; box on;
    title(figureTitles{iFig})
    ylim(ylimMat(iFig, :));
    xlim([plotDates(1) - 28, plotDates(end) + 56]);
    xticks(datenum(datetime(2008, 1, 1):calyears(2):datetime(2020, 1, 1)))
    datetick('x','yyyy','keeplimits', 'keepticks');
    ylabel('Percent');
    if iFig == 2 % Only produce legend for DFM, filtered estimates
       legend({'Baseline WEI', 'Alternative Specification'}, 'Location', 'southwest', 'Box','off') 
    end
    set(gca, 'FontSize', 16);
    Fig.Units = 'Inches';
    Fig.OuterPosition = [0 0 12.5 6];
    tightfig;
    print(Fig, fullfile(figPath, ['Figure4_Panel', num2str(iFig), '.png']), '-dpng', '-r300');
    print(Fig, fullfile(figPath, ['Figure4_Panel', num2str(iFig), '.pdf']), '-dpdf', '-r300');
    
    
end


%%% Creating Figure 5 - Plots WEI against 13-week moving average of WEI and
%%% GDP growth
WEIMA  = movmean(WEI, [12 0]); % 13 week moving average

d4lGDP    = gdpVintages.GDPC1_20210128_d4l;
gdpDates  = gdpVintages.DATE;
keepIndex = gdpDates >= datenum(2008, 1, 1);
gdpDates  = eomdate(year(gdpDates), month(gdpDates) + 2); % Plotting GDP at the end of quarter
gdpDates  = gdpDates(keepIndex); % Subsetting to only include post-2008 data
d4lGDP    = d4lGDP(keepIndex);

Fig = figure;
hold on;
plot(plotDates, WEI, 'LineWidth', 2);
plot(gdpDates,  d4lGDP, 'LineWidth', 2);
plot(plotDates, WEIMA, 'LineWidth', 2);
box on; grid on;
ylim([-12, 6]);
xlim([plotDates(1) - 28, plotDates(end) + 56]);
xticks(datenum(datetime(2008, 1, 1):calyears(2):datetime(2020, 1, 1)))
datetick('x','yyyy','keeplimits', 'keepticks');
ylabel('Percent');
legend({'WEI','GDP growth', '"WEI Nowcast"'}, 'Location',  'northeast', 'FontSize', 10);
set(gca, 'FontSize', 16);
Fig.Units = 'Inches';
Fig.OuterPosition = [0 0 15.42 6];
tightfig;
print(Fig, fullfile(figPath, 'Figure5.png'), '-dpng', '-r300');
print(Fig, fullfile(figPath, 'Figure5.pdf'), '-dpdf', '-r300');


            
%%% Creating Figure 8 - Plots WEI against WEI estimated on Full Sample and
%%% WEI using no claims data as an input series
Fig = figure;
hold on;
plot(plotDates, WEI, 'LineWidth', 1.5, 'LineStyle', '-', 'Color', [228,26,28] ./ 255);
plot(plotDates, weiSpecTable.WEI_full, 'LineWidth', 1.5, 'LineStyle', '--', 'Color', [228,26,28] ./ 255);
plot(plotDates, weiSpecTable.WEI_NoUI, 'LineWidth', 1.5, 'LineStyle', '-', 'Color', [55,126,184] ./ 255);
plot(plotDates, weiSpecTable.WEI_NoUI_full, 'LineWidth', 1.5, 'LineStyle', '--', 'Color', [55,126,184] ./ 255);
xline(datenum(2020, 3, 21), 'LineStyle','--', 'Color', 'k');
xlim([datenum(2020, 1, 1), datenum(2021, 1, 21)]);
xticks(datenum(datetime(2020, 1, 1):calmonths(2):datetime(2021, 1, 1)));
datetick('x', 'mmm-yyyy', 'keeplimits', 'keepticks');
ylabel('Percent');
box on;
legend({'WEI', 'WEI full sample', 'WEI No UI', 'WEI No UI full sample'}, 'Location', 'southoutside', 'NumColumns', 2)

Fig.Units = 'inches';
Fig.Position = [0 0 8 5];
tightfig;
print(Fig, fullfile(figPath,'Figure8.png'), '-dpng', '-r300');
print(Fig, fullfile(figPath,'Figure8.pdf'), '-dpdf', '-r300');


