function [lambda,factor,eigenvalue]= fhat(x,k, type)
% =========================================================================
% DESCRIPTION
% This code was adapted from McCracken and Ng, JBES, 2016.
% The program estimates a set of factors for a given dataset using
% principal component analysis. The number of factors estimated is
% determined by an information criterion specified by the user. 
%
% -------------------------------------------------------------------------
% INPUTS
%           x       = dataset (one series per column)
%           k       = number of factors to be estimated
%            
%                        
%           type  = an integer indicating the type of normalization on each 
%                    series in x before the factors are estimated;
%                    it can take on the following values:
%                           0 (no transformation)
%                           1 (demean only)
%                           2 (demean and standardize)
%                           
%
% OUTPUTS
%           lamhat  = factor loadings
%           factor  = set of factors
%           
%        eigenvalue = eigenvalues of x'*3 (where x is the dataset x post
%                     transformation and with missing values filled in)
% -------------------------------------------------------------------------
% SUBFUNCTIONS
%
% standardize() - standardize each series based on the chosen type of
% transformation
%
% pc() - runs principal component analysis
% -------------------------------------------------------------------------
% BREAKDOWN OF THE FUNCTION
%
% Part 1: Check that inputs are specified correctly.
%
% Part 2: Replace missing values with unconditional mean and standardize 
%         data 
%
% Part 3: Perform principal component analysis
%
% =========================================================================
% PART 1: CHECKS

% Check that x is not missing values for an entire row
if sum(sum(isnan(x),2)==size(x,2))>0
    error('Input x contains entire row of missing values.');
end

% Check that x is not missing values for an entire column
if sum(sum(isnan(x),1)==size(x,1))>0
    error('Input x contains entire column of missing values.');
end

% =========================================================================
% PART 2: SETUP

% Standardize based on the input value of type:  
%                  type  =  0 (no transformation)
%                           1 (demean only)
%                           2 (demean and standardize)


x_standardized = standardize(x, type);


% =========================================================================
% PART e: Perform principal component analysis on standardized data

[lambda,factor,eigenvalue]=pc(x_standardized,k);

end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         SUBFUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% =========================================================================
%                           standardize()
% =========================================================================
% Description
% This subfunction standardize based on the type specified:  
%                  type  =  0 (no transformation)
%                           1 (demean only)
%                           2 (demean and standardize)


function y_standardized = standardize(y, type)
% 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 
switch type
% ---------------------------------------------------------------------
    % Each series is demeaned only
    case 0
        y_standardized = y;
        
% ---------------------------------------------------------------------
    % Each series is demeaned only
    case 1
        y_standardized = y - mu;
        
% ---------------------------------------------------------------------
    % Each series is demeaned then diveded by standard deviation
    case 2
        y_standardized = (y - mu)./ sd;
end 
end


% =========================================================================
%                           pc()
% =========================================================================
% Description
% This function runs principal component analysis.
% -------------------------------------------------------------------------
% INPUTS
%           X      = standardized dataset (one series per column)
%           k      = number of factors 
%
% OUTPUTS
%  
%           factor = k by 1 factors estimates
%           lambda = N by k factor loadings 
%           chat = prediction of matrix x_standard (standardized x)


% =========================================================================
% FUNCTION
function [lambda,factor,eigenvalue]=pc(X,k)
[~,N] = size(X);

% Perform singular value decomposition
[U,S,~] = svd(X'*X);

% Factor loadings scaled by sqrt(N)
lambda = U(:,1:k)*sqrt(N);

% Factors scaled by 1/sqrt(N)
factor = X*lambda/N;

% eigenvalues of X'*X
eigenvalue = diag(S);
end
