%% Replication code for constructing WEI Nowcast Table and Figures
%  Note: Generates Table 6 panel a and panel b for 2018-2019 and RMSE
%  results (comptued using quarters in 2010-2019).

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


%% Loading in Data

%  Loading in WEI and previous vintage of WEI
WEISpec      = readtable(fullfile('..','Data', 'Outdata', 'wei_alt_spec_series.csv'));
WEI_latest   = WEISpec.WEI; % Baseline WEI rescaled using vintage of GDP of 01/28/2021
WEI_previous = WEISpec.WEI_GDP_0625; % WEI rescaled using vintage of GDP of 06/25/2020
WEIDates     = datenum(WEISpec.Date) - 6; % MATLAB datenums of WEI starting, marked as Sunday not Saturday

%  Loading in real time GDP estimates in levels
GDP_level = readtable(fullfile('..','Data','GDP_level_RT.csv')); 
for iCol = 2:size(GDP_level, 2) % Looping through and converting strings to numbers

    colName = GDP_level.Properties.VariableNames{iCol};
    GDP_level.(colName) = cellfun(@str2double, GDP_level.(colName)); 
end 

GDP_level = table2timetable(GDP_level); % Converting to a time table

%  Loading in latest GDP value
GDP      = readtable(fullfile('..','Data','gdp_vintage_data_full.csv'));
GDP      = timetable(GDP.GDPC1_20210128, 'RowTimes', datetime(GDP.DATE));

%% Computing 4-Quarter Growth, 2010 - 2019

% When to start and end each sample
startSamp = datetime(2010, 1, 1);
endSamp   = datetime(2019, 10, 1);

%  Computing 4-Q growth for each version of WEI
WEItt = timetable([WEI_latest, WEI_previous], 'RowTimes', datetime(WEIDates, 'ConvertFrom', 'datenum'));
WEIQ  = retime(WEItt, 'daily', 'previous');
WEIQ  = retime(WEIQ, 'quarterly', 'mean');
WEIQ.('Lag1')  = lagmatrix(WEIQ{:,:}, 1); % Incorprating lagged values

keepIndex = WEIQ.Time >= startSamp & WEIQ.Time <= endSamp;
WEIQ      = WEIQ(keepIndex, :);

% Computing 4-Q growth for realtime values
GDP_4Q_rt = ((GDP_level.Adv - GDP_level.L4Adv) ./ GDP_level.L4Adv) .* 100;

keepIndex = GDP_level.Date >= startSamp & GDP_level.Date <= endSamp;
GDP_4Q_rt = GDP_4Q_rt(keepIndex);

% Computing 4-Q growth for latest values
GDP_4Q_lt = ((GDP.Var1 - lagmatrix(GDP.Var1, 4)) ./ lagmatrix(GDP.Var1, 4)) .* 100;

keepIndex = GDP.Time >= startSamp & GDP.Time <= endSamp;
GDP_4Q_lt = GDP_4Q_lt(keepIndex);
%% Computing QoQ Growth Real Time, 2010 - 2019

%  Method 1
qoqGrowthRate = ((GDP_level{:, 2:5} - GDP_level{:, 1:4}) ./ GDP_level{:, 1:4}) + 1; % Computing 1 + QoQ growth rate (not annualized) for each lag
keepIndex     = GDP_level.Date >= startSamp & GDP_level.Date <= endSamp; % Subsetting for 2008-2019
qoqGrowthRate = qoqGrowthRate(keepIndex, :);

WEIQrt            = (WEIQ.Var1(:, 2) ./ 100) + 1; % Real time value of WEI as 4Q growth rate + 1
WEI_method1       = WEIQrt ./ (qoqGrowthRate(:, end) .* qoqGrowthRate(:, end-1) .* qoqGrowthRate(:, end-2)); % Dividing 4Q of WEI by lag 1, lag2, and lag 3 of GDP
WEI_method1_rt    = (((WEI_method1) .^ 4) - 1) .* 100; % Annualizing method 1 value and computing as percent

%  Method 2
GDP_lag_log_diff = log(GDP_level{:, 2}) - log(GDP_level{:, 1}); % ln GDP_q-4 - ln GDP_q-5
GDP_lag_log_diff = GDP_lag_log_diff(keepIndex);

WEI_method2    = 1 + ((WEIQ.Var1(:, 2) - WEIQ.Lag1(:, 2)) ./ 100) + GDP_lag_log_diff;
WEI_method2_rt = ((WEI_method2 .^ 4) - 1) * 100;

%  Advance
advGDP = (((GDP_level.Adv ./ GDP_level.L1Adv) .^ 4) .* 100) - 100;
advGDP = advGDP(keepIndex);

%% Computing QoQ Growth Latest, 2010 - 2019 
%  Using Latest vintage of GDP 01/28/2021

% Method 1
qoqGrowthRate = ((GDP.Var1 - lagmatrix(GDP.Var1, 1)) ./ lagmatrix(GDP.Var1, 1)) + 1; % Computing 1 + QoQ growth rate (not annualized)
qoqGrowthRate = lagmatrix(qoqGrowthRate, 3:-1:1); % Setting up necessary lags

keepIndex = GDP.Time >= startSamp & GDP.Time <= endSamp;
qoqGrowthRate = qoqGrowthRate(keepIndex, :);

WEIQrt            = (WEIQ.Var1(:, 1) ./ 100) + 1; % Latest value of WEI as 4Q growth rate + 1
WEI_method1       = WEIQrt ./ (qoqGrowthRate(:, end) .* qoqGrowthRate(:, end-1) .* qoqGrowthRate(:, end-2)); % Dividing 4Q of WEI by lag 1, lag2, and lag 3 of GDP
WEI_method1_lt    = (((WEI_method1) .^ 4) - 1) .* 100; % Annualizing method 1 value and computing as percent

% Method 2
GDP_lag_log_diff = log(lagmatrix(GDP.Var1, 4)) - log(lagmatrix(GDP.Var1, 5)); % ln GDP_q-4 - ln GDP_q-5
GDP_lag_log_diff = GDP_lag_log_diff(keepIndex);

WEI_method2        = 1 + ((WEIQ.Var1(:, 1) - WEIQ.Lag1(:, 1)) ./ 100) + GDP_lag_log_diff;
WEI_method2_lt     = ((WEI_method2 .^ 4) - 1) * 100;

% Latest
latestGDP  = ((GDP.Var1 ./ lagmatrix(GDP.Var1, 1)) .^ 4) * 100 - 100;
latestGDP  = latestGDP(keepIndex);

%% Constructing into tables and computing RMSEs

% YoY Table
yoyTable      = array2table([datenum(WEIQ.Time), GDP_4Q_lt, GDP_4Q_rt, WEIQ.Var1(:,2)], 'VariableNames', {'Date','Latest','Advance', 'WEI'});
yoyTable.Date = datestr(yoyTable.Date, 'mm/dd/yyyy');

% QoQ Table 
qoqTable  = array2table([datenum(WEIQ.Time), latestGDP, advGDP, WEI_method1_rt, WEI_method2_rt, WEI_method1_lt, WEI_method2_lt], 'VariableNames', ...
            {'Date','Latest', 'Advance', 'Method1_RealTime', 'Method2_RealTime', 'Method1_Latest','Method2_Latest'});
qoqTable.Date = datestr(qoqTable.Date, 'mm/dd/yyyy');

% RMSE YoY
rmseYoY = mean((repmat(yoyTable{:, 2}, 1, 2) - yoyTable{:, 3:end}) .^ 2, 1) .^ (1/2);
rmseYoY = array2table(rmseYoY, 'VariableNames', yoyTable.Properties.VariableNames(3:end));

% RMSE QoQ
rmseQoQ = mean((repmat(qoqTable{:, 2}, 1, 5) - qoqTable{:, 3:end}) .^ 2, 1) .^ (1/2);
rmseQoQ = array2table(rmseQoQ, 'VariableNames', qoqTable.Properties.VariableNames(3:end));

%% Writing tables

writetable(yoyTable, fullfile('..', 'Data', 'Outdata', 'wei_nowcast_2010-2019_yoy.csv'));
writetable(qoqTable, fullfile('..', 'Data','Outdata', 'wei_nowcast_2010-2019_qoq.csv'));
writetable(rmseYoY,  fullfile('..', 'Data','Outdata', 'wei_nowcast_2010-2019_rmse_yoy.csv'));
writetable(rmseQoQ,  fullfile('..', 'Data','Outdata', 'wei_nowcast_2010-2019_rmse_qoq.csv'));