%% Replication code for generating WEI Nowcasts in 2020 and Progression Plots
%  Constructs 2020 portion of Table 6 and generates subfigures for Figure 7

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

%% Vintage Data Path and Data Information

vintDataPath = fullfile('..','Data','WEI_Input_Vintage_Data'); % Where we load the vintage Data
lastDate     = datenum(2021, 1, 21); % Last vintage to consider when constructing real-time values of WEI
vintDates    = datenum(2020, 3, 12):7:lastDate; % List of vintage dates to consider (all Thursdays)

% Storage for objects of interest
targetQ   = [datenum(2020, 1, 1), datenum(2020, 4, 1), datenum(2020, 7, 1), datenum(2020, 10, 1)]; % Which quarters to consider forecasting
endTargQ  = eomdate(year(targetQ), month(targetQ) + 2); % End of period dates for quarters being forecasted
meth1Mat  = NaN(length(vintDates), length(targetQ)); % QoQ forecasts generated using method 1
meth2Mat  = NaN(length(vintDates), length(targetQ)); % QoQ forecasts generated using method 2
yoyMat    = NaN(length(vintDates), length(targetQ)); % YoY forecasts 

% When to start and end estimation for the WEI
estStart     = datenum(2008, 1, 1);
estEnd       = datenum(2020, 2, 29);

% Series that need to be transformed to d52l
seriesTransID = {'LICN','PSTL', 'ASASI', 'RSTOTL', 'LIUN', 'PELOUS', 'FUEL', 'RAS'};


%% Setting up vintage data for GDP

gdpVintageTable = readtable(fullfile('..', 'Data', 'gdp_vintage_data_full.csv'));

% Which vintages do we have?
gdpVintages = cellfun(@(x) datenum(x(end-7:end), 'yyyymmdd'), gdpVintageTable.Properties.VariableNames(2:end));
[gdpVintages, sortIndex] = sort(gdpVintages);
gdpVintageTable(:,2:end) = gdpVintageTable(:, sortIndex + 1);
gdpVintageTable.Properties.VariableNames(:,2:end) = gdpVintageTable.Properties.VariableNames(:, sortIndex + 1);

% Creating d4l version
gdpVintageMat = gdpVintageTable{:, 2:end};
d4lGDPMat     = (log(gdpVintageMat) - log(lagmatrix(gdpVintageMat, 4))) * 100;

% Cutting this off only for the quarters of interest (2008 to end of
% 2019)
gdpVintageTable.DATE = datenum(gdpVintageTable.DATE);
whichGDPIndex = gdpVintageTable.DATE >= datenum(2008, 1, 1) & gdpVintageTable.DATE <= datenum(2019, 10 , 1);
d4lGDPMat     = d4lGDPMat(whichGDPIndex, :);

% Creating a table for qoq percent change that we'll need to back out the
% Method 1 and Method 2 values of QoQ WEI Nowcasts
pcGDPMat   = gdpVintageTable{:, 2:end};
pcGDPMat   = ((pcGDPMat - lagmatrix(pcGDPMat, 1)) ./ lagmatrix(pcGDPMat, 1)) * 100;
pcGDPTable = gdpVintageTable;
pcGDPTable{:, 2:end} = pcGDPMat;

%% Looping through vintages and generating Nowcasts

for iVint = 1:length(vintDates)
    
    vintDate = vintDates(iVint); % Which vintage
    
    % Reading in vintage data
    vintData = readtable(fullfile(vintDataPath, ['Data_', datestr(vintDate, 'yyyy_mmdd'), '.csv']));
    
    % Transforming Series
    transIndex             = find(ismember(vintData.Properties.VariableNames, seriesTransID));
    vintData{:,transIndex} = (log(vintData{:,transIndex}) - log(lagmatrix(vintData{:,transIndex}, 52))) * 100;
    vintData.Date          = datenum(vintData.Date);
    
    % WEI start and end dates
    weiStart = datenum(2007, 12, 30);
    weiEnd   = vintDate - weekday(vintDate) - 6; % Starting Sunday of the last week

    % Estimating weights for WEI (Note: weights not identified to sign)
    estKeepIndex = vintData.Date >= estStart & vintData.Date <= estEnd;
    XEst         = vintData{estKeepIndex, 2:end}; % Subsetting data
    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     = pca(XNormEst); % Using PCA to find loadings of series (factors -> series)
    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 = vintData.Date >= weiStart & vintData.Date <= weiEnd;
    XFull         = vintData{fullKeepIndex, 2:end}; % Subsetting data
    datesFull     = vintData.Date(fullKeepIndex);   % Subsetting dates 
    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;
    
    % WEI over period of estimation for partial 
    keepIndex  = datesFull >= weiStart & datesFull <= estEnd;
    XNormEstim = XNormFull(keepIndex, :);
    weiEstim   = weiFactor(keepIndex, :);
    
    % Constructing WEI for dates with missing data
    sIndex = find(~isnan(weiFactor), 1, 'last') + 1; % Last week where there is a full set of data + 1
    eIndex = length(weiFactor); % Last week to consider
    
    for iRef = sIndex:eIndex % Looping through and running regression to generate partial data WEI
        
        dataSlice              = XNormFull(iRef, :);      % Week of data with missing data 
        nonMissingIndex        = find(~isnan(dataSlice)); % Series that are not missing this week
        
        % RHS - Constant, 2 lags of WEI, current value of nonmissing
        % series, LHS - WEI
        X = [ones(size(XNormEstim, 1), 1), lagmatrix(weiEstim, 1:2), XNormEstim(:, nonMissingIndex)];
        Y = weiEstim;
        
        B = regress(Y, X); % Running regression

        % Forecasting WEI for the week
        input                = [1, weiFactor([iRef - 1,iRef - 2])', XNormFull(iRef, nonMissingIndex)];
        weiFactor(iRef)      = input * B;       

    end
    
    % Rescaling WEI to GDP (using correct vintage)
    gdpIndex     = find(gdpVintages <= vintDate, 1, 'last'); % Which vintage of GDP to use
    gdpRescale   = d4lGDPMat(:, gdpIndex);
    
    weiTT = timetable(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(datenum(2008, 1, 1), 'ConvertFrom', 'datenum') & ...
                                 weiTT.Time <= datetime(datenum(2019, 10, 1), 'ConvertFrom', 'datenum'));

    % Generating rescaling coefficients (regressing quarterly WEI on GDP)
    rescaleCoeff = regress(gdpRescale, [ones(size(weiRescale)), weiRescale]);
    weiGDP       = [ones(size(weiFactor, 1), 1), weiFactor] * rescaleCoeff;
    
    % Generating quarterly average of WEI
    if datesFull(end) < eomdate(year(targetQ(end)), month(targetQ(end)) + 2) % Extend WEI if last date for the last quarter not present
       
        weiGDP    = [weiGDP; NaN];
        datesFull = [datesFull; eomdate(year(targetQ(end)), month(targetQ(end)) + 2)];
    end
    weiGDPTT  = timetable(weiGDP, 'RowTimes', datetime(datesFull, 'ConvertFrom', 'datenum'));
    weiGDPTT  = retime(weiGDPTT, 'daily', 'previous');
    weiGDPTT  = retime(weiGDPTT, 'quarterly', 'mean');
    
    % Generating Nowcasts using WEI (yoy, QoQ method 1, QoQ method 2)
    pcValues   = pcGDPTable{:, gdpIndex + 1}; % Value of GDP from current vintage in percent change
    pcDates    = pcGDPTable.DATE;             % Dates of GDP
    
    for iTarg = 1:length(targetQ)  % Looping through and generating forecasts for each quarter
        
        if isnan(pcValues(pcDates == targetQ(iTarg))) % Only run if the value for this quarter does not exist

            
           targetDT = datetime(targetQ(iTarg), 'ConvertFrom', 'datenum'); % Target Quarter as a datetime
           
           yoyMat(iVint, iTarg) = weiGDPTT.weiGDP(weiGDPTT.Properties.RowTimes == targetDT);  % Assinging annalized forecast  
           
           % Constructing Method 1 Value
           wei4Q       = weiGDPTT.weiGDP(weiGDPTT.Properties.RowTimes == targetDT) ./ 100;
           notNaNIndex = find(~isnan(pcValues), 1, 'last');              % Used to take first 3 lags of pc values (either actual or forecast)
           last3Q      = pcValues((notNaNIndex - 2):notNaNIndex) ./ 100; % Note: pc values are not annualized (no need to raise to 1/4)
           last3Q      = last3Q + 1;                                     % Adding back one to compute final number
           
           weiMeth1 = (1 + wei4Q) ./ (last3Q(1) * last3Q(2) * last3Q(3));    % Backing out QoQ value, using equation (9)
           weiMeth1 = (weiMeth1 ^ 4) - 1;
           weiMeth1 = 100 * weiMeth1; 
           
           meth1Mat(iVint, iTarg)    = weiMeth1; 
           pcValues(notNaNIndex + 1) = (((1 + (weiMeth1 ./ 100)) .^ 0.25) - 1) * 100; % Assinging forecast value of WEI in consistent units (will need for Method 1 computation of multiple horizons ahead)
                                                                                      % Note: Adding this value will not affect computation of Method 2
           
           % Constructing Method 2 Value
           tgQIndex = find(weiGDPTT.Properties.RowTimes == targetDT);
           lQIndex  = tgQIndex - 1;
           
           wei4Q_T   = weiGDPTT.weiGDP(tgQIndex);     % Current Value of WEI (quarterly average)
           wei4Q_TL1 = weiGDPTT.weiGDP(tgQIndex - 1); % Lagged  Value of WEI (quarterly average)
           
           gdpVals   = table2array(gdpVintageTable(:, gdpIndex + 1)); % GDP values in levels
           currIndex = find(gdpVintageTable.DATE == targetQ(iTarg));
           GDP_TL4   = gdpVals(currIndex - 4); % GDP_t-4
           GDP_TL5   = gdpVals(currIndex - 5); % GDP_t-5
           
           weiMeth2                = 1 + ((wei4Q_T - wei4Q_TL1) ./ 100) + log(GDP_TL4) - log(GDP_TL5); % Backing out QoQ value, using equation (10)
           weiMeth2                = 100 * ((weiMeth2 ^ 4) - 1);
           meth2Mat(iVint, iTarg)  = weiMeth2;
           
        end

    end

end

%% Constructing Nowcast Table for 2020
%  Note: creating rows in Table 6 for the quarters of 2020

%%% Constructing Latest Value of GDP (YoY and QoQ)
yoyLatest = ((gdpVintageTable.GDPC1_20210128 - lagmatrix(gdpVintageTable.GDPC1_20210128, 4)) ./ lagmatrix(gdpVintageTable.GDPC1_20210128, 4)) .* 100;
qoqLatest = ((((gdpVintageTable.GDPC1_20210128) ./ lagmatrix(gdpVintageTable.GDPC1_20210128, 1)) .^ 4) - 1) .* 100;
keepIndex = gdpVintageTable.DATE >= datenum(2020, 1, 1) & gdpVintageTable.DATE <= datenum(2020, 10, 1);
yoyLatest = yoyLatest(keepIndex);
qoqLatest = qoqLatest(keepIndex);


%%% Constructing Advance Value of GDP (YoY and QoQ)
advColNames = {'GDPC1_20200429', 'GDPC1_20200730' , 'GDPC1_20201029', 'GDPC1_20210128'};
yoyAdvance  = NaN(4, 1);
qoqAdvance  = NaN(4, 1);
for iQ = 1:length(advColNames)
   
    gdpVals   = gdpVintageTable.(advColNames{iQ});
    lastIndex = find(~isnan(gdpVals), 1, 'last');
    yoy       = ((gdpVals - lagmatrix(gdpVals, 4)) ./ lagmatrix(gdpVals, 4)) .* 100;
    qoq       = ((((gdpVals) ./ lagmatrix(gdpVals, 1)) .^ 4) - 1) .* 100;
    
    yoyAdvance(iQ) = yoy(lastIndex);
    qoqAdvance(iQ) = qoq(lastIndex);
    
end

%%% Constructing WEI using Method 1/2 using the latest vintage of GDP
weiMeth1Latest = NaN(4, 1);
weiMeth2Latest = NaN(4, 1);
gdpLatest      = gdpVintageTable.GDPC1_20210128;

for iQ = 1:length(weiMeth1Latest)
    
    targetQuarter = targetQ(iQ);
    tqIndex       = find(targetQuarter == datenum(weiGDPTT.Time));
    wei4Q         = weiGDPTT.weiGDP(tqIndex);
    wei4QLag      = weiGDPTT.weiGDP(tqIndex - 1);
    
    % Constructing Method 1 Value
    pcVal  = gdpLatest ./ lagmatrix(gdpLatest, 1); % 1 + percent change in GDP (not annualized and not multiplied by 100)
    gIndex = find(gdpVintageTable.DATE == targetQuarter);   
    pcVal  = pcVal(gIndex-3:gIndex-1); % First 3 lags of GDP
    m1     = (1 + (wei4Q ./ 100)) ./ (prod(pcVal)); % Give 1 + implied percent change in GDP (not annualized and not mulitplied by 100)
    m1     = ((m1 .^ 4) - 1) .* 100; % Annualizing, subtracting by 1, and multiplying by 100
    
    % Constructing Method 2 Value
    gdp_TL4 = gdpLatest(gIndex - 4);
    gdp_TL5 = gdpLatest(gIndex - 5);
    m2      = 1 + ((wei4Q - wei4QLag) ./ 100) + log(gdp_TL4) - log(gdp_TL5); % RHS of equation (10)
    m2      = ((m2 .^ 4) - 1) .* 100; % Solving for g_hat in equation 10 and multiplying by 100
    
    % Assigning values
    weiMeth1Latest(iQ) = m1;
    weiMeth2Latest(iQ) = m2;
    
end

%%% Subsetting correct values of real time values of yoy, method 1, and
%%% method 2

yoyRT            = NaN(4, 1);
weiMeth1RealTime = NaN(4, 1);
weiMeth2RealTime = NaN(4, 1);

for iQ = 1:length(yoyRT)
    
    yoy = yoyMat(:, iQ);
    m1  = meth1Mat(:, iQ);
    m2  = meth2Mat(:, iQ);

    selectIndex = find(~isnan(yoy), 1, 'last');
    
    yoyRT(iQ) = yoy(selectIndex);
    weiMeth1RealTime(iQ)  = m1(selectIndex);
    weiMeth2RealTime(iQ)  = m2(selectIndex);
    
end

%%% Constructing and writing table
dates2020 = gdpVintageTable.DATE(keepIndex); 
yoyTable  = array2table([dates2020, yoyLatest ,yoyAdvance, yoyRT], 'VariableNames', {'Date', 'Latest', 'Advance', 'WEI'});
qoqTable  = array2table([dates2020, qoqLatest ,qoqAdvance, weiMeth1RealTime, weiMeth2RealTime, weiMeth1Latest, weiMeth2Latest],...
            'VariableNames', {'Date','Latest', 'Advance', 'Method1_RealTime', 'Method2_RealTime', 'Method1_Latest','Method2_Latest'});

yoyTable.Date = datestr(yoyTable.Date, 'mm/dd/yyyy');
qoqTable.Date = datestr(qoqTable.Date, 'mm/dd/yyyy');
        
writetable(yoyTable, fullfile('..', 'Data', 'Outdata', 'wei_nowcast_2020_yoy.csv'));
writetable(qoqTable, fullfile('..', 'Data', 'Outdata', 'wei_nowcast_2020_qoq.csv'));


%% Loading in data for progression plots

%%% Setting up struct to store each forecaster's forecast
for iTarg = 1:length(targetQ)
   ForStruct(iTarg) = struct;
end

%%% Loading in Atlanta Fed GDP Now

% Parsing in data
gdpnowTable     = readtable(fullfile('..', 'Data', 'GDPTrackingModelDataAndForecasts_TrackingArchives.csv'));
forecastDates   = datenum(gdpnowTable.ForecastDate); % When the forecast was made
quarterForecast = datenum(gdpnowTable.QuarterBeingForecasted); % Which quarter was being forecasted

% Looping through and adding dates and forecasts for 2020Q1-2020Q3
for iTarg = 1:(length(targetQ) - 1)
        
    qIndex = quarterForecast == endTargQ(iTarg);
    ForStruct(iTarg).GDPNow        = [forecastDates(qIndex), gdpnowTable.GDPNowcast(qIndex)]; 
       
end

% Constructing forecasts for 2020Q4
gdpnowCurrCell      = readcell(fullfile('..', 'Data', 'GDPTrackingModelDataAndForecasts_TrackingHistory.csv'));
gdpnowDates         = gdpnowCurrCell(1, 3:end); % Cleaning dates for forecast dates for GDP now
forecastDates       = NaN(length(gdpnowDates), 1);
for iDate = 1:length(gdpnowDates) % Dates formatted just dd-mmm, no year, so adding these
   
    gdpnowDateStr = gdpnowDates{iDate};
    monthStr      = gdpnowDateStr(end-2:end);
    if strcmp(monthStr, 'Jan')
       yearStr = '2021';  
    else
       yearStr = '2020';
    end
    forecastDates(iDate) = datenum([gdpnowDateStr, '-', yearStr]);
    
end                
forecastVals        = cell2mat(gdpnowCurrCell(31, 3:end));
ForStruct(4).GDPNow = [forecastDates, forecastVals'];


%%% Loading in New York Fed Staff Nowcast

% Reading in and cleaning excel sheet into a table
nowcastTable  = readcell(fullfile('..', 'Data', 'New-York-Fed-Staff-Nowcast_data_2002-present_ForecastsByQuarter.csv'));
nowcastTable  = array2table(nowcastTable(2:end, :), 'VariableNames', [{'ForecastDate'}, cellfun(@(x) ['Nowcast_', x], nowcastTable(1, 2:end), 'UniformOutput', false)]);

% Correct nowcast excel sheet dating convention
correctDates = [datenum(2020, 4, 28), datenum(2020, 7, 29), datenum(2020, 10, 28), datenum(2021, 1, 27)]; % Actual last date of the nowcast, not date of publication

% Looping through and adding forecasts/dates
for iTarg = 1:length(targetQ)
    
    colStr               = ['Nowcast_', num2str(year(targetQ(iTarg))), 'Q', num2str(quarter(targetQ(iTarg)))];
    nowcastVal           = nowcastTable.(colStr);
    nanIndex             = cellfun(@ismissing, nowcastVal);
    nowcastVal(nanIndex) = cellfun(@(x) NaN, nowcastVal(nanIndex), 'UniformOutput', false);
    nowcastVal           = cell2mat(nowcastVal);
    
    nowcastDates                                                       = cellfun(@(x) x + calyears(2000), nowcastTable.ForecastDate); % Making sure year is correct
    nowcastDates                                                       = datenum(nowcastDates);
    nowcastDates(find(nowcastDates > correctDates(iTarg), 1, 'first')) = correctDates(iTarg); % Correcting nowcast excel sheet dating convention for last forecast
    
    ForStruct(iTarg).Nowcast = [nowcastDates(~isnan(nowcastVal)), nowcastVal(~isnan(nowcastVal))];
    
end

%%% Loading in Bluechip

forecast = readtable(fullfile('..','Data','bluechip.csv'));
forecast.ForecastDate  = datenum(forecast.ForecastDate);
forecast.ReferenceDate = datenum(forecast.ReferenceDate);

% Looping through and adding bluechip forecastss/dates
gdpReleaseDates = [datenum(2020, 4, 29), datenum(2020, 7, 30), datenum(2020, 10, 29), datenum(2021, 1, 28)]; % Release Dates
for iTarg = 1:length(targetQ)
   
    bcVal   = forecast.Value(forecast.ReferenceDate == endTargQ(iTarg)); 
    bcDates = forecast.ForecastDate(forecast.ReferenceDate == endTargQ(iTarg));
    
    [bcDates,dateIndex] = sort(bcDates);
    bcVal               = bcVal(dateIndex);
    bcIndex             = bcDates >= datenum(datetime(targetQ(iTarg), 'ConvertFrom', 'datenum') - calmonths(1)) &  bcDates < gdpReleaseDates(iTarg);
    
    ForStruct(iTarg).BlueChip = [bcDates(bcIndex), bcVal(bcIndex)];
end

%%% Incorporating WEI Forecasts from Method 1
for iTarg = 1:length(targetQ) % Looping through and adding forecasts/dates
   
    weiVals  = meth1Mat(:,iTarg);
    weiDates = vintDates';
    
    % Finding the first date in the quarter (ignoring forecasts from before)
    startQIndex = find(weiDates >= targetQ(iTarg), 1, 'first');
    lastQIndex  = find(~isnan(weiVals), 1, 'last');
    
    ForStruct(iTarg).WEI = [weiDates(startQIndex:lastQIndex), weiVals(startQIndex:lastQIndex)];
    
end

%%% Incorporating Advance, Latest, and date for each GDP release
advance     = qoqTable.Advance;
latest      = qoqTable.Latest;
latest(end) = NaN; % Because first revision has not come out yet for 2020:Q4 GDP
for iTarg = 1:length(targetQ)
    
    ForStruct(iTarg).Advance = advance(iTarg); 
    ForStruct(iTarg).Latest  = latest(iTarg);
    ForStruct(iTarg).GDPDate = gdpReleaseDates(iTarg);
    
end

%% Creating Progression Plots

% What colors to use for each forecaster?
colorMat = [
228,26,28
55,126,184
77,175,74
152,78,163] ./ 255;

%%% Progression Plot for 2020Q1
iTarg = 1;
Fig   = figure;
hold on;
plot(ForStruct(iTarg).WEI(:,1), ForStruct(iTarg).WEI(:,2), 'LineWidth', 1.5, 'Marker', 'diamond', 'Color', colorMat(1, :), 'MarkerFaceColor', colorMat(1, :));
plot(ForStruct(iTarg).BlueChip(:,1), ForStruct(iTarg).BlueChip(:,2), 'LineWidth', 1.5, 'Marker', 'diamond', 'Color', colorMat(2, :), 'MarkerFaceColor', colorMat(2, :));
plot(ForStruct(iTarg).Nowcast(:,1), ForStruct(iTarg).Nowcast(:,2), 'LineWidth', 1.5, 'Marker', 'diamond', 'Color', colorMat(3, :), 'MarkerFaceColor', colorMat(3, :));
plot(ForStruct(iTarg).GDPNow(:,1), ForStruct(iTarg).GDPNow(:,2), 'LineWidth', 1.5, 'Marker', 'diamond', 'Color', colorMat(4, :), 'MarkerFaceColor', colorMat(4, :));
scatter(ForStruct(iTarg).GDPDate, ForStruct(iTarg).Advance ,120, 'Marker', 'o', 'MarkerEdgeColor', 'k');
scatter(ForStruct(iTarg).GDPDate, ForStruct(iTarg).Latest, 120, 'Marker', 's', 'MarkerEdgeColor', 'k');

xlim([datenum(2019, 11, 28), datenum(2020, 5, 5)]);
datetick('x', 'mmm-yyyy', 'keeplimits');
legLoc = 'southwest';
legend({'WEI', 'Blue Chip', 'NY Fed Staff Nowcast', 'GDPNow', 'Advance Release', 'Latest'}, ...
        'Location', legLoc, 'Orientation', 'vertical', 'EdgeColor', [0, 0, 0]);
box on;
grid on;
set(gca, 'FontSize', 13);
ylabel('Percent (QoQ Annualized Percent Change)')
title([num2str(year(targetQ(iTarg))), ':Q', num2str(quarter(targetQ(iTarg))), ' Real GDP Forecasts'], 'FontSize', 16);
Fig.Units = 'inches';
Fig.Position(1:4) = [0, 0, 10, 6];
tightfig;
saveas(Fig, fullfile('..','Figures', 'Figure7_Panel1.png'));
print(Fig, fullfile('..','Figures', 'Figure7_Panel1.png'), '-dpdf', '-r300');


%%% Progression Plot for 2020Q2 - 2020Q4
xlimMat = [datenum(2020, 2, 28), datenum(2020, 8, 2);... % x limits of each forecasting quarter
           datenum(2020, 5, 30), datenum(2020, 11, 5);...
           datenum(2020, 8, 30), datenum(2021, 2, 5)];
ylimMat = [-60, 10;...
           -20, 40;...
           -1,  13];
legLocStrs = {'northeast','northwest', 'southwest'};     % Legend Locations
for iTarg = 2:4
    Fig   = figure;
    hold on;
        
    WEIDates = ForStruct(iTarg).WEI(:, 1); % Dividing values of WEI into pre-2020Qx adv. release and post
    WEIVals  = ForStruct(iTarg).WEI(:, 2);
    WEIVals  = WEIVals(2:end);  % Removing values from before we had any information for the Qx
    WEIDates = WEIDates(2:end); 

    pQReleaseDate = gdpReleaseDates(iTarg - 1); % Date of advance gdp release date of prior quarter
    cutIndex      = find(WEIDates < pQReleaseDate, 1, 'last') + 1;
    WEIDatesP1    = WEIDates(1:cutIndex);
    WEIValsP1     = WEIVals(1:cutIndex); % Pre advance GDP release date
    WEIDatesP2    = WEIDates(cutIndex:end);
    WEIValsP2     = WEIVals(cutIndex:end); % Post advance GDP release date

    p0 = plot(WEIDatesP1, WEIValsP1, 'LineStyle', '--', 'LineWidth', 1.5, 'Marker', 'diamond', 'Color', colorMat(1, :), 'MarkerFaceColor', colorMat(1, :));
    p1 = plot(WEIDatesP2, WEIValsP2, 'LineStyle', '-', 'LineWidth', 1.5, 'Marker', 'diamond', 'Color', colorMat(1, :), 'MarkerFaceColor', colorMat(1, :));
    p2 = plot(ForStruct(iTarg).BlueChip(:,1), ForStruct(iTarg).BlueChip(:,2), 'LineWidth', 1.5, 'Marker', 'diamond', 'Color', colorMat(2, :), 'MarkerFaceColor', colorMat(2, :));
    p3 = plot(ForStruct(iTarg).Nowcast(:,1), ForStruct(iTarg).Nowcast(:,2), 'LineWidth', 1.5, 'Marker', 'diamond', 'Color', colorMat(3, :), 'MarkerFaceColor', colorMat(3, :));
    p4 = plot(ForStruct(iTarg).GDPNow(:,1), ForStruct(iTarg).GDPNow(:,2), 'LineWidth', 1.5, 'Marker', 'diamond', 'Color', colorMat(4, :), 'MarkerFaceColor', colorMat(4, :));
    p5 = scatter(ForStruct(iTarg).GDPDate, ForStruct(iTarg).Advance ,120, 'Marker', 'o', 'MarkerEdgeColor', 'k');
    p6 = scatter(ForStruct(iTarg).GDPDate, ForStruct(iTarg).Latest, 120, 'Marker', 's', 'MarkerEdgeColor', 'k');

    xlim(xlimMat(iTarg - 1, :));
    ylim(ylimMat(iTarg - 1, :));
    datetick('x', 'mmm-yyyy', 'keeplimits');
    legLoc = legLocStrs{iTarg - 1};
    legend([p1, p2, p3, p4, p5, p6], {'WEI', 'Blue Chip', 'NY Fed Staff Nowcast', 'GDPNow', 'Advance Release', 'Latest'},...
           'Location', legLoc, 'Orientation', 'vertical', 'EdgeColor', [0, 0, 0]);
    box on;
    grid on;
    ylabel('Percent (QoQ Annualized Percent Change)')
    set(gca, 'FontSize', 13);
    title([num2str(year(targetQ(iTarg))), ':Q', num2str(quarter(targetQ(iTarg))), ' Real GDP Forecasts'], 'FontSize', 16);
    Fig.Units = 'inches';
    Fig.Position(1:4) = [0, 0, 10, 6];
    tightfig;
    saveas(Fig, fullfile('..','Figures', ['Figure7_Panel', num2str(iTarg), '.png']));
    print(Fig, fullfile('..','Figures', ['Figure7_Panel', num2str(iTarg), '.pdf']), '-dpdf', '-r300');

end
