% 0.PCAoverTFPgrowth.m
%
% Asymptotic principal component analysis (Connor and Korajczyk, 1987) on
% firm-level TFP growth (Section 2.2 of the paper).
%
% Input : output\TFP.csv produced by 0.TFPestimation.do (firm-year panel).
% Output: output\TFP_pca.csv with columns {cyear, F1-F10}.
%         The main analysis uses the first six principal components (k=6);
%         F1-F10 are exported for the robustness checks in Appendix H.1.
diary 0.PCAoverTFPgrowth.txt
clear;

% Set the project folder before running.
projpath = 'D:\SAS\TFP\Dispersion_replication';

%Import firm-level TFP estimation.
%Keep {gvkey, cyear, d_TFPpool} only. Previously this three-column sheet was
%prepared manually in Excel; readtable now extracts the columns directly.
T = readtable(fullfile(projpath,'output','TFP.csv'));
M = [T.gvkey, T.cyear, T.d_TFPpool];
%Matrix M - Column{gvkey, year, Pooled TFP growth}
%Convert matrix M to Cell
C = num2cell(M,1);
%Make the unbalanced panel matrix
year = unique(C{1,2});
gvkey = unique(C{1,1});
%Make the tfp unbalanced matrix
tfp = zeros(length(gvkey), length(year));
for i = 1 : length(C{1,2})
    tfp(find(gvkey == C{1,1}(i)), C{1,2}(i)) = C{1,3}(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Our sample has 3,977 firms, so 1:3977 selects all of them.
%Equivalently, 1:length(gvkey) can be used.
tfp1 = tfp(1:3977,1971:2016);
%Replace the missing values with zeros
tfp1(isnan(tfp1))=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nrc, ncc] = size(tfp1);
notmissing = tfp1 ~= 0;
notmissingt = notmissing*ones(ncc,1);
cont = notmissingt >=11; %the minimum number of notmissing values for each firm.
%Reformat the matrix with the tfp obs. more than cont.
tfp1 = tfp1(cont,:);
notmissing = notmissing(cont,:);
% Asymptotic PCA
% Calculate Omega
Omega = (tfp1'*tfp1)./(real(notmissing)'*real(notmissing));
k=10;
[F, ev] = eig(Omega);
[nrg, ncg] = size(F);
F = (ncg^.5)*F(:,end:-1:end-k+1);
time = [1971:1:2016]';
Total=[time F];
tfp_pca = array2table(Total,'Variablenames',{'cyear','F1','F2','F3','F4','F5','F6','F7','F8','F9','F10'});
writetable(tfp_pca,fullfile(projpath,'output','TFP_pca.csv'));
diary off