%% Replication code for monthly and quarterly regression
%  Constructs output necessary for Table 4 and Table 5

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

%% Loading in WEI, GDP, and IP

% Loading in Baseline WEI
weiSpec  = readtable(fullfile('..','Data', 'Outdata', 'wei_alt_spec_series.csv'));
WEI      = weiSpec.WEI;
weiDates = datenum(weiSpec.Date);
weiDates = weiDates - 6; % Making into sunday (need for retiming)
% Loading in vintages of GDP/transforming to d4l (will only need latest
% vintage)
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 in and transforming IP into d12l (yoy % change)
IP      = readtable(fullfile('..','Data','IP.csv'));
pcIP    = ((IP.INDPRO_20210115 - lagmatrix(IP.INDPRO_20210115, 1)) ./ lagmatrix(IP.INDPRO_20210115, 1)) .* 100;
IPdates = datenum(IP.DATE);


%% Constructing table needed for quarterly regressions

% When do we want values in the table to start from and end?
startTableDate = datenum(2008, 1, 1);
endTableDate   = datenum(2019, 10, 1);
endMonthDate   = datenum(2019, 12, 1);

%  Constructing a quarterly average of WEI
weiTT = array2timetable(WEI, 'RowTimes', datetime(weiDates, 'ConvertFrom', 'datenum'));
weiTT = retime(weiTT, 'daily','previous'); 
weiTT = retime(weiTT, 'quarterly', 'mean');
WEIQ  = weiTT.WEI(weiTT.Time >= datetime(startTableDate, 'ConvertFrom', 'datenum') &...
                  weiTT.Time <= datetime(endTableDate,   'ConvertFrom', 'datenum'));

% Constructing a  monthly average of WEI
weiTT = array2timetable(WEI, 'RowTimes', datetime(weiDates, 'ConvertFrom', 'datenum'));
weiTT = retime(weiTT, 'daily','previous'); 
weiTT = retime(weiTT, 'monthly', 'mean');
WEIM  = weiTT.WEI(weiTT.Time >= datetime(startTableDate, 'ConvertFrom', 'datenum') &...
                  weiTT.Time <= datetime(endMonthDate,   'ConvertFrom', 'datenum'));

WEIM1 = WEIM(1:3:end); % First month in quarter
WEIM2 = WEIM(2:3:end); % Second month in quarter
WEIM3 = WEIM(3:3:end); % Third month in quarter

% Generating lags of GDP and subsetting
GDP       = lagmatrix(gdpVintages.GDPC1_20210128_d4l, 4:-1:0);
keepIndex = gdpVintages.DATE >= startTableDate & gdpVintages.DATE <= endTableDate;
GDP       = GDP(keepIndex, :);
qDates    = gdpVintages.DATE(keepIndex);
% Creating a table based on these values
qTable = array2table([qDates, WEIQ, WEIM1, WEIM2, WEIM3, GDP], 'VariableNames', {'Date', 'WEIQ', 'WEIM1', 'WEIM2', 'WEIM3', 'GDPL4', 'GDPL3', 'GDPL2', 'GDPL1','GDP'});

%% Constructing table needed for monthly regressions

% When do we want values in the table to start from and end?
startTableDate = datenum(2008, 1, 1);
endTableDate   = datenum(2020, 2, 1);

% Re-subsetting monthly aggregate of WEI
WEIM  = weiTT.WEI(weiTT.Time >= datetime(startTableDate, 'ConvertFrom', 'datenum') &...
                  weiTT.Time <= datetime(endTableDate,   'ConvertFrom', 'datenum'));
              
% Constructing pseudo-weekly measures
weiTT    = array2timetable(WEI, 'RowTimes', datetime(weiDates, 'ConvertFrom', 'datenum'));
weiTT    = retime(weiTT, 'daily','previous'); 

weiTable       = timetable2table(weiTT); 
weiTable.Time  = datenum(weiTable.Time);
weiTable.Year  = year(weiTable.Time);
weiTable.Month = month(weiTable.Time);

weiTable.PWEEK                     = ceil(day(weiTable.Time) / 7); % Pseudo week
weiTable.PWEEK(weiTable.PWEEK > 4) = 4; % Adjusting for months with more than 28 days
gMean = grpstats(weiTable, {'Year', 'Month', 'PWEEK'}, 'mean','DataVars',{'WEI'}); 

gMean = gMean(gMean.Year >= year(startTableDate) & gMean.Year <= year(endTableDate), :);
WEIP1 = gMean.mean_WEI(1:4:end); % psuedo week 1
WEIP2 = gMean.mean_WEI(2:4:end); % pseudo week 2 
WEIP3 = gMean.mean_WEI(3:4:end); % pseudo week 3
WEIP4 = gMean.mean_WEI(4:4:end); % pseudo week 4

pDates = datenum(gMean.Year, gMean.Month, 1); % getting dates
pDates = pDates(1:4:end); % same dates every 4 weeks, only need first

% Subsetting pseudo week data
keepIndex = pDates >= startTableDate & pDates <= endTableDate;
pDates    = pDates(keepIndex);
WEIP      = [WEIP1, WEIP2, WEIP3, WEIP4];
WEIP      = WEIP(keepIndex, :);

% Generating lags of IP and susbetting
IP        = lagmatrix(pcIP, 4:-1:0);
keepIndex = IPdates >= startTableDate & IPdates <= endTableDate;
IP        = IP(keepIndex, :);

mTable = array2table([pDates, WEIP, WEIM, IP], 'VariableNames', {'Date', 'WEIP1', 'WEIP2', 'WEIP3', 'WEIP4', 'WEIM',  'IPL4', 'IPL3', 'IPL2', 'IPL1', 'IP'});


%% Running quarterly regressions

%  Which RHS Specifications to consider? (Note: each will inclued a
%  constant)

nSpec   = 4;
RHSSpec = cell(4, 1); 
RHSSpec{1,1} = {'WEIQ',  'GDPL4', 'GDPL3', 'GDPL2', 'GDPL1'};
RHSSpec{2,1} = {'WEIM1', 'GDPL4', 'GDPL3', 'GDPL2', 'GDPL1'};
RHSSpec{3,1} = {'WEIM1', 'WEIM2', 'GDPL4', 'GDPL3', 'GDPL2', 'GDPL1'};
RHSSpec{4,1} = {'WEIM1', 'WEIM2', 'WEIM3', 'GDPL4', 'GDPL3', 'GDPL2', 'GDPL1'};

QRegressStruct = struct; % Saving Results in a Struct

T      = size(qTable, 1); % Number of observations
Lag = 4; % How many lags of GDP included
for iSpec = 1:nSpec
    
    spec = RHSSpec{iSpec, 1};
    QRegressStruct(iSpec).XVarName = [spec, {'Const.'}]; 
    
    % Running regression
    XTable = qTable(:, spec);
    X      = [XTable{:,:}, ones(size(qTable, 1), 1)];
    Y      = qTable{:,end};
    [coeff, ~, resid, ~, stats] = regress(Y, X);
    
    % Information from regression    
    weiCoeffIndex = 1:(length(spec) - Lag);
    m             = weiCoeffIndex(end); % Last index of WEI RHS 
    V     = (X'*X)^-1 * (X'.*resid') * (resid .* X) * (X'*X)^-1; % Computing variance (Heteroscedastic Corrected)
    SE    = sqrt(diag(V(1:m, 1:m))); % Standard error of WEI coeff
    pval  = 2 * (1 - tcdf(abs(coeff(1:m )./ sqrt(diag(V(1:m,1:m)))), T - (m + Lag + 1)));
    R2    = stats(1);
    ADJR2 = 1 - (1-R2) * (T - 1)/(T - (m + Lag) - 1); % Adjusted R2
    SER   = std(resid); % Standard deviation of residuals
    if  iSpec > 1 % Should we compute F-Test of monthly coeff = 0
        
        L      = [eye(m) zeros(m, Lag + 1)];
        F1stat = (coeff' * L' * ((L*V*L')^-1) * L * coeff) / m;
        F1pval = fcdf(F1stat, m , T - (m + Lag + 1), 'upper');
        
    else
       
        F1stat = NaN;
        F1pval = NaN;
        
    end
    if iSpec > 2 % Should we compute F-Test of monthly coeff being equal
        
        L      = [ones(m - 1,1) -eye(m - 1) zeros(m - 1,Lag + 1)];
        F2stat = (coeff' * L' * ((L*V*L')^-1) * L * coeff) / m;
        F2pval = fcdf(F2stat, m , T - (m + Lag + 1),'upper');

    else
       
        F2stat = NaN;
        F2pval = NaN;
        
    end
    
    % Allocating value to struct
    QRegressStruct(iSpec).Coeff  = coeff;
    QRegressStruct(iSpec).SE     = SE;
    QRegressStruct(iSpec).pval   = pval;
    QRegressStruct(iSpec).ADJR2  = ADJR2;
    QRegressStruct(iSpec).SER    = SER;
    QRegressStruct(iSpec).F1stat = F1stat;
    QRegressStruct(iSpec).F1pval = F1pval;
    QRegressStruct(iSpec).F2stat = F2stat;
    QRegressStruct(iSpec).F2pval = F2pval;
    
end



%% Running monthly regressions


%  Which RHS Specifications to consider? (Note: each will inclued a
%  constant)

nSpec   = 5;
RHSSpec = cell(nSpec, 1); 
RHSSpec{1,1} = {'WEIM',  'IPL4',  'IPL3',  'IPL2', 'IPL1'};
RHSSpec{2,1} = {'WEIP1', 'IPL4',  'IPL3',  'IPL2', 'IPL1'};
RHSSpec{3,1} = {'WEIP1', 'WEIP2', 'IPL4',  'IPL3', 'IPL2', 'IPL1'};
RHSSpec{4,1} = {'WEIP1', 'WEIP2', 'WEIP3', 'IPL4', 'IPL3', 'IPL2',  'IPL1'};
RHSSpec{5,1} = {'WEIP1', 'WEIP2', 'WEIP3', 'WEIP4', 'IPL4', 'IPL3', 'IPL2', 'IPL1'};

MRegressStruct = struct; % Saving Results in a Struct

T      = size(mTable, 1); % Number of observations
Lag    = 4; % How many lags of GDP included
for iSpec = 1:nSpec
    
    spec = RHSSpec{iSpec, 1};
    MRegressStruct(iSpec).XVarName = [spec, {'Const.'}]; 
    
    % Running regression
    XTable = mTable(:, spec);
    X      = [XTable{:,:}, ones(size(mTable, 1), 1)];
    Y      = mTable{:,end};
    [coeff, ~, resid, ~, stats] = regress(Y, X);
    
    % Information from regression    
    weiCoeffIndex = 1:(length(spec) - Lag);
    m             = weiCoeffIndex(end); % Last index of WEI RHS 
    V     = (X'*X)^-1 * (X'.*resid') * (resid .* X) * (X'*X)^-1; % Computing variance (Heteroscedastic Corrected)
    SE    = sqrt(diag(V(1:m, 1:m))); % Standard error of WEI coeffs
    pval  = 2 * (1 - tcdf(abs(coeff(1:m )./ sqrt(diag(V(1:m,1:m)))), T - (m + Lag + 1)));
    R2    = stats(1);
    ADJR2 = 1 - (1-R2) * (T - 1)/(T - (m + Lag) - 1); % Adjusted R2
    SER   = std(resid); % Standard deviation of residuals
    if  iSpec > 1 % Should we compute F-Test of monthly coeff = 0
        
        L      = [eye(m) zeros(m, Lag + 1)];
        F1stat = (coeff' * L' * ((L*V*L')^-1) * L * coeff) / m;
        F1pval = fcdf(F1stat, m , T - (m + Lag + 1), 'upper');
        
    else
       
        F1stat = NaN;
        F1pval = NaN;
        
    end
    if iSpec > 2 % Should we compute F-Test of monthly coeff being equal
        
        L      = [ones(m - 1,1) -eye(m - 1) zeros(m - 1,Lag + 1)];
        F2stat = (coeff' * L' * ((L*V*L')^-1) * L * coeff) / m;
        F2pval = fcdf(F2stat, m , T - (m + Lag + 1),'upper');

    else
       
        F2stat = NaN;
        F2pval = NaN;
        
    end
    
    % Allocating value to struct
    MRegressStruct(iSpec).Coeff  = coeff;
    MRegressStruct(iSpec).SE     = SE;
    MRegressStruct(iSpec).pval   = pval;
    MRegressStruct(iSpec).ADJR2  = ADJR2;
    MRegressStruct(iSpec).SER    = SER;
    MRegressStruct(iSpec).F1stat = F1stat;
    MRegressStruct(iSpec).F1pval = F1pval;
    MRegressStruct(iSpec).F2stat = F2stat;
    MRegressStruct(iSpec).F2pval = F2pval;
    
end


