%% Matlab code to combine data from various sources and write to .mat files

clear all
plotdata = 1 ;                              % Set to 1 to produce plots

%% Philadelphia Fed Real Time Output Data
% rows = quarters 1953Q2 - 2023Q4 (283 rows), columns = monthly vintages 1965M11 to 2024M4 (702 vintages)

DAT = xlsread('routputMvQd.xlsx') ;            

% R1 indicator based on "advance data" (middle month of quarter)
% Vintages start 1965M11, so v=1 is 1st vintage for 1965Q4
R1AD = [] ;                                 % Store recession indicators here
for v = 1:3:700 ;                           % Do loop to pick middle month of each quarter
    yv = DAT(:,v) ;                         % Output for vintage v
    T = length(yv) ;                        % Number of obs
    Rv = NaN(T,1) ;                         % Store R1 for vintage v here, 1953Q2 missing
    for t = 2:T ;                           % Begin do loop through time
        if yv(t) < yv(t-1) ;                % If output < previous period....
            Rv(t) = 1 ;                     % ....set R1 to 1.
        elseif yv(t) >= yv(t-1) ;           % If output >= previous period....
            Rv(t) = 0 ;                     % ....set R1 to 0
        end ;                           
    end ;                                   % End do loop through time
    R1AD = [ R1AD Rv ] ;                    % Store indicator for current vintage
    clear v yv T Rv t 
end ;                                       % End do loop through vintages
R1AD = [ NaN(283,50) R1AD ] ;               % Pad with missing values for vintages for 1953Q2-1965Q3

% R1 indicator based on "first final" (3rd month in each quarter)
R1FF = [] ;                                 % Store recession indicators here
for v = 2:3:701 ;                           % Do loop to pick last month of each quarter
    yv = DAT(:,v) ;                         % Output for vintage v
    T = length(yv) ;                        % Number of obs
    Rv = NaN(T,1) ;                         % Store R1 for vintage v here, 1953Q2 missing
    for t = 2:T ;                           % Begin do loop through time
        if yv(t) < yv(t-1) ;                % If output < previous period....
            Rv(t) = 1 ;                     % ....set recession indicator to 1.
        elseif yv(t) >= yv(t-1) ;           % If output >= previous period....
            Rv(t) = 0 ;                     % ....set R1 to 0.
        end ;                           
    end ;                                   % End do loop through time
    R1FF = [ R1FF Rv ] ;                    % Store indicator for current vintage
    clear v yv T Rv t 
end ;                                       % End do loop through vintages
clear v
R1FF = [ NaN(283,50) R1FF ] ;               % Pad with missing values for vintages for 1953Q2-1965Q3

% Create R1 indicator based on "final vintage" (April 2024)

yv = DAT(:,end) ;                           % Output for data vintage v
T = length(yv) ;                            % Number of obs
R1FV = NaN(T,1) ;                           % Store R1 here, 1953Q2 missing
for t = 2:T ;                               % Begin do loop through time
     if yv(t) < yv(t-1) ;                   % If output < previous period....
        R1FV(t) = 1 ;                       % ....set R1 to 1.
    elseif yv(t) >= yv(t-1) ;               % If output >= previous period....
        R1FV(t) = 0 ;                       % ....set R1 to 0
    end ;                           
end ;                                       % End do loop through time
clear yv T t DAT 

% Create vector of dates 1953Q2 - 2023Q4

qtr = linspace(1953.25,2023.75,size(R1AD,1))' ;   % Dates 1953Q2 - 2023Q4

% Plot R1 indicator for last vintage

if plotdata == 1 ;

    figure('position', [100, 100, 700, 1100])  

    subplot(3,1,1)
    bar(qtr,R1AD(:,end),'b')
    title('R1 Indicator (advance data)')

    subplot(3,1,2)
    bar(qtr,R1FF(:,end),'b')
    title('R1 Indicator (first final)')

    subplot(3,1,3)
    bar(qtr,R1FV,'b')
    title('R1 Indicator (final vintage)')

end ;

%% Load FRED Data 1934Q1 - 2023Q4

DAT = xlsread('FRED.xlsx','Data') ;            
TB3 = DAT(:,2) ;                                % 3 Month TBill Yield 
GS10 = DAT(:,3) ;                               % 10 Year Bond Yield - starts 1953Q2
CFNFCI = DAT(:,4) ;                             % Chicago Fed NFCI - starts 1971Q1
disp('Check indicies used to identify dates')
disp(['For 1953Q2 GS10 & TB3 should be 2.9967 and 2.1533.              Check: ' num2str([GS10(78) TB3(78)]) ])
disp(['For 2023Q4 GS10, TB3 & CFNFCI should be 4.44, 5.28 & -0.39074.  Check: ' num2str([GS10(360) TB3(360) CFNFCI(360)]) ])
Spread = GS10(78:360) - TB3(78:360) ;           % Spread 1953Q2 - 2023Q4
NFCI = CFNFCI(78:360) ;                         % NFCI 1953Q2 - 2023Q4 (Missing = NaN before 1971Q1)
clear DAT TB3 GS10 CFNFCI

% Create negative Term Structure Dummy
T = size(Spread,1) ;
SpreadNeg = zeros(T,1) ;                        % Vector of zeros (Spread > 0)
for t = 1:T ;
    if Spread(t) < 0 ;                          % If Spread < 0
        SpreadNeg(t) = 1 ;                      % Set SpreadNeg = 1
    end ;
end  ;
clear T t

% Plot Spread, Inversion Dummy and NFCI

if plotdata == 1 ;

    figure('position', [100, 100, 700, 1100])  
    subplot(2,1,1)
    yyaxis left ;   ax = gca ;  ax.YColor = 'r' ;
    bar(qtr,SpreadNeg,'r')
    yyaxis right ;  ax.YColor = 'k' ;
    plot(qtr,Spread,'linewidth',2,'Color','k')
    title('Term Spread and Inversion Dummy')

    subplot(2,1,2)
    plot(qtr(72:end),NFCI(72:end),'linewidth',2)
    axis([1971 2024 -2 5])
    title('NFCI')

end ;

%% Load realtime NFCI data from McCracken's Website
% Weekly data from Jan 8, 1971 to Jul 23, 2021 (2638 weeks) 
% 1971 & 1972 missing until Sep 7, 2017 vintage
% Weekly vintages from Jan 6, 1988 to July 28, 2021 (1750 vintages)
% Take the 7th week of each quarter to match advance data release of GDP
% And 11th week from each quarter to match first final release of GDP
% Then average the 13 weeks for each quarter to get quarterly observations
% Note real time NFCI data ends in 2021Q2

DAT = xlsread('NFCI_vintages.xlsx','B2:BOI2636') ;
[ NQ, NV ] = size(DAT) ;

% Create Quarterly Advance Data series -- use 7th week of each quarter
WeeklyFCI = [] ;                                        % Store quarterly vintages of weekly data here (1988Q1-2021Q2)
for v = 7:13:NV-4 ;                                     % Do loop to pick 7th week of each quarter
    WeeklyFCI = [ WeeklyFCI DAT(:,v) ] ;                % Store current 7th week as current quarter vintage of weekly data
end ;
[ ~, NVQ ] = size(WeeklyFCI) ;
NFCIAD = [] ;                                           % Store quarterly real time NFCI data here (will be 202 quarters, 1971Q1 - 2021Q2)
for v = 1:NVQ ;                                         % Do loop through quarterly vintages
    VFCI = [] ;                                         % Store quarterly data for current vintage here
    for t = 13:13:NQ ;                                  % Do loop through quarters                           
        VFCI = [ VFCI ; mean(WeeklyFCI(t-12:t,v)) ] ;   % Current quarterly obs is average of weekly obs
    end ;
    NFCIAD = [ NFCIAD VFCI ] ;                          % Store current vintage
end ;
NFCIAD = [ NFCIAD ; NaN(10,134) ] ;                     % Pad with missing values for 2021Q3-2023Q4 
NFCIAD = [ NaN(212,60) NFCIAD NaN(212,11)] ;            % Pad with missing values for vintages 1973Q1 - 1987Q4 and 2021Q3 - 2024Q1
clear WeeklyFCI v NVQ VFCI t 

% Create Quarterly First Final series -- use 11th week of each quarter
WeeklyFCI = [] ;                                        % Store quarterly vintages of weekly data here (1988Q1-2021Q2)
for v = 11:13:NV ;                                      % Do loop to pick 11th week of each quarter
    WeeklyFCI = [ WeeklyFCI DAT(:,v) ] ;                % Store current 11th week as current quarter vintage of weekly data
end ;
[ ~, NVQ ] = size(WeeklyFCI) ;
NFCIFF = [] ;                                           % Store quarterly real time NFCI data here (will be 202 quarters, 1971Q1 - 2021Q2)
for v = 1:NVQ ;                                         % Do loop through quarterly vintages
    VFCI = [] ;                                         % Store quarterly data for current vintage here
    for t = 13:13:NQ ;                                  % Do loop through quarters                           
        VFCI = [ VFCI ; mean(WeeklyFCI(t-12:t,v)) ] ;   % Current quarterly obs is average of weekly obs
    end ;
    NFCIFF = [ NFCIFF VFCI ] ;                          % Store current vintage
end ;
NFCIFF = [ NFCIFF ; NaN(10,134) ] ;                      % Pad with missing values for 2021Q3-2023Q4                
NFCIFF = [ NaN(212,60) NFCIFF NaN(212,11)] ;             % Pad with missing values for vintages 1973Q1 - 1987Q4 and 2021Q3 - 2024Q1
clear WeeklyFCI v NVQ VFCI t DAT NQ NV

if plotdata == 1 ;

    % Plot FRED NFCI & Real Time NFCI 1971Q1-2021Q2 for 2021Q2 vintage 
    figure('position', [100, 100, 700, 1100])  
    subplot(2,1,1)
    plot(NFCIFF(:,end-11),'linewidth',1)
    title('NFCI - First Final (2021Q2 Vintages) 1971 - 2022')
    subplot(2,1,2)
    plot(NFCIAD(:,end-11),'linewidth',1)
    title('NFCI - Advance Data (2021Q2 Vintages) 1971 - 2022')

end ; 

%% Load SPF Data from Philadelphia Fed 1968Q4 - 2023Q4

DAT = xlsread('SPF.xlsx','Data') ;
spf0 = DAT(:,3)/100  ;            % Prob recession current quarter 
spf1 = DAT(:,4)/100  ;            % Prob recession one quarter ahead 
spf2 = DAT(:,5)/100  ;            % Prob recession two quarters ahead 
spf3 = DAT(:,6)/100  ;            % Prob recession three quarters ahead 
spf4 = DAT(:,7)/100  ;            % Prob recession four quarters ahead 
SPF = [ spf0 spf1 spf2 spf3 spf4 ] ;
clear spf0 spf1 spf2 spf3 spf4 DAT

if plotdata == 1 ;
    % Plot SPFs
    figure('position', [100, 100, 1100, 600])   
    plot(qtr(63:end),SPF,'linewidth',2)
    legend('h=0','h=1','h=2','h=3','h=4')
    title('SPF Recession Forecasts')
    axis([1968 2024 0 1])
end ;

%% Write Data to Files

% Full Sample

save('YSRep.mat', 'qtr', 'R1AD', 'R1FF', 'R1FV', 'Spread', 'SpreadNeg', 'SPF')

% 1971 - 2023: add NFCI (final vintage) & inversion dummy
qtr = qtr(72:end) ;                             % Data starts in 1953Q2, so 1971Q1 is obs #72
Spread = Spread(72:end) ;                       % Spread
NFCI = NFCI(72:end) ;                           % Final Vintage NFCI
SPF = SPF(10:end,:) ;                           % SPF forecasts start 1968Q4 so 1971Q1 is obs #10
[NQ,NV] = size(R1FF) ;                          % Size of real time data matrices
R1AD = R1AD(72:NQ,72:NV) ;                      % Initial release
R1FF = R1FF(72:NQ,72:NV) ;                      % First final
R1FV = R1FV(72:NQ) ;                            % Final vintage
save('YSRep1971.mat', 'qtr', 'R1AD', 'R1FF', 'R1FV', 'Spread', 'NFCI', 'SPF')

% 1973 - 2023: add real time NFCI
qtr = qtr(9:end) ;                              % Previously shortened data starts in 1971Q1, so 1973Q1 is obs #9
Spread = Spread(9:end) ;                        % Spread
NFCI = NFCI(9:end) ;                            % Final Vintage NFCI
NFCIAD = NFCIAD(9:end,:) ;                      % Real Time NFCI - initial release
NFCIFF = NFCIFF(9:end,:) ;                      % Real Time NFCI - first final
SPF = SPF(9:end,:) ;                            % SPF forecasts
[NQ,NV] =size(R1FF) ;                           % Size of real time data matrices
R1AD = R1AD(9:NQ,9:NV) ;                        % Initial release
R1FF = R1FF(9:NQ,9:NV) ;                        % First final
R1FV = R1FV(9:NQ) ;                             % Final vintage
save('YSRep1973.mat', 'qtr', 'R1AD', 'R1FF', 'R1FV', 'Spread', 'NFCI', 'NFCIAD', 'NFCIFF', 'SPF')
