%% Replication code for plotting WEI against alternative indicators
%  Constructs Figure 3 and Figure 6
%  NOTE: Relies on output generated from wei_alt_spec.m (Generates Baseline
%  WEI Specification in wei_alt_spec_series.csv)

% 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 WEI and GDP

% Loading in Baseline WEI
weiSpec  = readtable(fullfile('..','Data','Outdata', 'wei_alt_spec_series.csv'));
WEI      = weiSpec.WEI;
weiDates = datenum(weiSpec.Date);

% Loading in vintages of GDP/transforming to d4l, will need for rescaling
% of ADS and CFNAI
gdpVintages = readtable(fullfile('..', 'Data','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);


%% Loading and transforming ADS

% Loading in and transforming ADS
ADS           = readtable(fullfile('..','Data','ADS_Index_Most_Current_Vintage.csv'));
ADS.Var1      = datetime(ADS.Var1);
ADS           = table2timetable(ADS, 'RowTimes','Var1'); % Converting to time-table
ADS.ADS_Index = movmean(ADS.ADS_Index, [364, 0]); % 365 moving average
ADS_quarter   = retime(ADS, 'quarterly','mean'); % Converting to quarterly series

% Subsetting ADS for regression 
rescaleStart = datenum(2008, 1, 1);
rescaleEnd   = datenum(2019, 10, 1);

qDates    = datenum(ADS_quarter.Var1);
keepIndex = qDates >= rescaleStart & qDates <= rescaleEnd;
ADSQ      = ADS_quarter.ADS_Index(keepIndex);

% Subsetting GDP for regression
GDP       = gdpVintages.GDPC1_20210128_d4l;
gdpDates  = gdpVintages.DATE;
keepIndex = gdpDates >= rescaleStart & gdpDates <= rescaleEnd;
GDP       = GDP(keepIndex);

% Regressing GDP on quarterly ADS (2008 through 2019)
ADScoeff = regress(GDP, [ones(size(ADSQ)), ADSQ]);

% Rescaling and subsetting ADS
ADSrescale = [ones(size(ADS.ADS_Index)), ADS.ADS_Index] * ADScoeff;
keepIndex  = ADS.Var1 >= datetime(2008, 1, 1) & ADS.Var1 <= datetime(weiDates(end), 'ConvertFrom', 'datenum');
ADSrescale = ADSrescale(keepIndex);
ADSdates   = datenum(ADS.Var1(keepIndex)); % Need for subsequent plotting

%% Loading and transforming CFNAI

% Loading in and transforming CFNAI
CFNAI         = readtable(fullfile('..','Data','CFNAI.csv'));

CFNAI         = CFNAI(11:end, 1:2); % Cleaning data
CFNAI.Properties.VariableNames = {'observation_date', 'CFNAI'};
CFNAI.observation_date = cellfun(@datetime, CFNAI.observation_date);
CFNAI.CFNAI            = cellfun(@str2double, CFNAI.CFNAI);

CFNAI         = table2timetable(CFNAI, 'RowTimes','observation_date');
CFNAI.CFNAI   = movmean(CFNAI.CFNAI, [11, 0]) .* 12; % Sum over past 12 months
CFNAI_quarter = retime(CFNAI, 'quarterly','mean');

% Subsetting quarterly CFNAI for rescaling regression
qDates      = datenum(CFNAI_quarter.observation_date);
keepIndex   = qDates >= rescaleStart & qDates <= rescaleEnd;
CFNAIQ      = CFNAI_quarter.CFNAI(keepIndex);

% Regressing GDP on quarterly CFNAI (2008 through 2019)
CFNAIcoeff = regress(GDP, [ones(size(CFNAIQ)), CFNAIQ]); % Note: GDP already subsetted correctly

% Rescaling and subsetting ADS
CFNAIrescale = [ones(size(CFNAI.CFNAI)), CFNAI.CFNAI] * CFNAIcoeff;
keepIndex    = CFNAI.observation_date >= datetime(2008, 1, 1) & CFNAI.observation_date <= datetime(weiDates(end), 'ConvertFrom', 'datenum');
CFNAIrescale = CFNAIrescale(keepIndex);
CFNAIdates   = datenum(CFNAI.observation_date(keepIndex)) + 15; % Need for subsequent plotting, adding 15 days to get roughly middle of the month

%% Loading and transforming IP

% Loading in and transforming IP into d12l (yoy % change)
IP      = readtable(fullfile('..','Data','IP.csv'));
d12l_IP = (log(IP.INDPRO_20210115) - log(lagmatrix(IP.INDPRO_20210115, 12)) ) * 100; % Transfomring to YoY % change
IPdates = datenum(IP.DATE);

% Putting IP into WEI units
weiKeep = weiDates >= datenum(2008, 1, 1) & weiDates <= datenum(2020, 2, 29);
weiMean = nanmean(WEI(weiKeep));
weiStd  = nanstd(WEI(weiKeep), 0, 1);

IPKeep  = IPdates >= datenum(2008, 1, 1) & IPdates <= datenum(2020, 2, 1);
IPMean  = nanmean(d12l_IP(IPKeep));
IPStd   = nanstd(d12l_IP(IPKeep));

d12l_IP = (d12l_IP - IPMean) ./ IPStd;
d12l_IP = (d12l_IP .* weiStd) + weiMean;

% Subsetting IP for post 2008
IPKeep  = IPdates >= datenum(2008, 1, 1) & IPdates <= datenum(2021, 1, 1);
d12l_IP = d12l_IP(IPKeep);
IPdates = IPdates(IPKeep);
IPdates = eomdate(IPdates);
%% Creating Figures

figPath = fullfile('..','Figures');

%%% Generating Figure 3 - WEI vs ADS vs CFNAI
Fig = figure;
hold on;
plot(weiDates, WEI, 'LineWidth', 2);
plot(ADSdates, ADSrescale, 'LineWidth', 2);
plot(CFNAIdates, CFNAIrescale, 'LineWidth', 2);
grid on; box on;
ylim([-12, 6]);
xlim([weiDates(1) - 28, weiDates(end) + 56]);
xticks(datenum(datetime(2008, 1, 1):calyears(2):datetime(2020, 1, 1)))
datetick('x','yyyy','keeplimits', 'keepticks');
ylabel('Percent');
legend({'WEI','ADS', 'CFNAI'}, 'Location',  'northwest', 'FontSize', 12);
set(gca, 'FontSize', 16);
Fig.Units = 'Inches';
Fig.OuterPosition = [0 0 15 6];
tightfig;
print(Fig, fullfile(figPath, 'Figure3.png'), '-dpng', '-r300');
print(Fig, fullfile(figPath, 'Figure3.pdf'), '-dpdf', '-r300');

%%% Generating Figure 6
Fig = figure;
hold on;
plot(weiDates, WEI, 'LineWidth', 2);
plot(IPdates, d12l_IP, 'LineWidth', 2);
grid on; box on;
ylim([-12, 6]);
xlim([weiDates(1) - 28, weiDates(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 15 6];
tightfig;
print(Fig, fullfile(figPath, 'Figure6.png'), '-dpng', '-r300');
print(Fig, fullfile(figPath, 'Figure6.pdf'), '-dpdf', '-r300');