clear
close all
clc

% =========================================================================
% DESCRIPTION 
% This script loads in cyclical component, factor and factor loading,
% then calculate the percentage of cyclical component explained
% by the first r factors. 

% -------------------------------------------------------------------------
% INPUT FILES 
%       hx_cyc_comp_'year'.mat - MAT file with cyclical component from OLS
%       detrending procedure
%       hxfactor_'year'.mat - MAT file with HX factors and factor loadings
% OUTPUT FILES
%       Rsquared'year'.xlsx: xlsx with R-squared explained by the first  
%       r HX factors. 
% -------------------------------------------------------------------------
% Control parameters

% Choose the data file. 
version = 2024;

% =========================================================================
% PART 1: Load data

% file name of cyclical component data
input_cyc = strcat('hx_cyc_comp_',string(version),'.mat'); 

% file name of factor data
input_hx = strcat('hxfactor_',string(version),'.mat'); 


% read data
load(input_cyc)

load(input_hx)


% =========================================================================
% PART 2: Calculate percentage of cyclical component being explained by 
%         factors.

% Demean and divide each cyclical component by standard deviation
cyc_hx = standardize(vals); 

% Calculate percentage explained by factors 
% perc_f: N by R matrix with R = number of factors. The rth column gives 
% the percentage of each variable explained by factor 1 to r. 

R2 = R_squared(fhat_hx,lambda_hx,cyc_hx);

% =========================================================================
% Save output
output = strcat('Rsquared',string(version),'.xlsx');
writecell(varnames', output, 'Sheet',1,'Range','A1')
writecell(des, output, 'Sheet',1,'Range','B1')
R2_cell = num2cell(R2);
writecell(R2_cell,output,'Sheet',1,'Range','C1')

%% ========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         SUBFUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Standardize cyclical component 

function y_standardized = standardize(y)

% =========================================================================
% DESCRIPTION
% This function demean and devide by standard deviation of each column
%
% -------------------------------------------------------------------------
% INPUTS
%           y                  = dataset (one series per column)  
%
% OUTPUTS
%           y_standardized     = each series demeaned and standardized by 
%                                standard deviation 
% =========================================================================
% FUNCTION

% Get unconditional mean of the non-missing values of each series
mu = repmat(nanmean(y),size(y, 1),1);

% calculate standard deviation of each series (in matrix)
sd = repmat(std(y),size(y, 1),1);


% Perform transformation based on type 

y_standardized = (y - mu)./ sd;
end

%% Subfunction for Out-of-Sample Inference    
 

function R2 = R_squared(Fhat,lamhat,cvals)

% =========================================================================
% DESCRIPTION
% This function calculates the percentage of y explained by factors
%
% -------------------------------------------------------------------------
% INPUT
        % Fhat    = estimated factors (one factor per column)
        % lamhat  = factor loadings (one factor per column)
% OUTPUT
        % perc_f  = sum of squared predicted cvals using up to r factors devided by
        % sum of squared cvals

% =========================================================================
% FUNCTION

% N = number of series, R = number of factors, T = number of observations
[N,R] = size(lamhat); 
T = size(Fhat,1);

% Initiate output matrix 
R2 = NaN*ones(N,R); 

% Initiate a matrix to keep estimated cyclical values by factors
chat = NaN*ones(T,R);

for r = 1:R
    for i = 1:N
        chat(:,i,r) = Fhat(:,1:r)* transpose(lamhat(i,1:r));  
    end   
end

chat2 = sum(chat.^2, 1);
cvals2 = sum(cvals.^2,1);
for r = 1:R
    R2(:,r) = transpose(chat2(:,:,r)./cvals2);
end 
    
end
