clear
clc
close all

addpath functions

set(0,'DefaultAxesFontName', 'Times New Roman')
set(0,'DefaultAxesFontSize', 24)
set(0,'DefaultTextFontname', 'Times New Roman')
set(0,'DefaultTextFontSize', 24)

load data.mat

userets={'wrds','ff25','ind','op_inv','me_var','me_ni','me_beta','me_ac','me_mom','br_bonds','br_fx'};

R=[];
for j=1:length(userets)
    R=[R returns.(userets{j})];
end
R=R(:,~isnan(sum(R)));

[~,~,eigret]=pca(bsxfun(@minus,R,mean(R)));
 
figure
subplot(1,2,1);
plot(1:20,eigret(1:20)','r*');
title('Eigenvalues 1 to 20')
xlabel('Eigenvalue');
xlim([1 20]);
subplot(1,2,2);
plot(5:20,eigret(5:20)','r*');
xlim([5 20])
xlabel('Eigenvalue');
title('Eigenvalues 5 to 20')
set(gcf, 'units', 'inches', 'position', [2 2 12 7])
set(gcf, 'Color', 'w');
save2pdf('output/eigenvalues')
close

% Construct factor models. Recall all models are added the FF3 factors later on.
factors.ff3=[factors.rmrf factors.smb factors.hml];
factors.ff5=[factors.rmw factors.cma];
factors.intermediary=[factors.intermediarycapitalriskfactor factors.aemleveragefactor];
factors.novimarx=[factors.tmax factors.glbt factors.nino factors.ssn];

% -----------------------------------------------------------------
% Analysis by factor
factorlist={'rmrf','smb','hml','mom','ff5','bab','qmj','liquidity','intermediary','novimarx','ipgrowth','jln','cons','consstock'};
lab={{'Market'},{'SMB'},{'HML'},{'MOM'},{'RMW','CMA'},{'BAB'},{'QMJ'},{'Liquidity'},{'Interm. (He)','Interm. (Adrian)'},...
    {'NY temp.','Global temp.','El Ni\~no','Sunspots'},{'IP growth'},{'Macro PC 1','Macro PC 2','Macro PC 3'},{'Cons. growth'},{'Stockholder cons.'}};

tradedfactors={'rmrf','smb','hml','mom','ff5','bab','qmj'};

pmin=5;
pmax=20;
nwlags=3;
formatSpec = '%10.2f';

TABLE1=cell(0);
TABLE2=cell(0);

for j=1:length(factorlist)

    factorlist{j}
    
    % -----------------------------------------
    % Results using equity + bonds + currencies
    % -----------------------------------------
    
    userets={'wrds','ff25','ind','op_inv','me_var','me_ni','me_beta','me_ac','me_mom','br_bonds','br_fx'};
    
    R=[];
    group=[];
    for i=1:length(userets)
        R=[R returns.(userets{i})];
        group=[group i*ones(1,size(returns.(userets{i}),2))];
    end
    group=group(1,~isnan(sum(R)));
    R=R(:,~isnan(sum(R)));
    
    % Number of factors
    nf=size(factors.(factorlist{j}),2);
    
    % Controls and test portfolios
    R1=factors.rmrf;
    R2=[factors.rmrf factors.smb factors.hml];
    
    % Avoid duplication of factors when adding the market or the FF3 factors
    if strcmp(factorlist{j},'rmrf')
        G1=factors.rmrf;
        G2=factors.ff3;
    elseif strcmp(factorlist{j},'smb')
        G1=[factors.smb R1] ;
        G2=[factors.smb factors.rmrf factors.hml];
    elseif strcmp(factorlist{j},'hml')
        G1=[factors.hml R1] ;
        G2=[factors.hml factors.rmrf factors.smb];
    else
        G1=[factors.(factorlist{j}) R1];
        G2=[factors.(factorlist{j}) R2];
    end
    
    rf=factors.rf;
    
    nonmissingobs=find(~isnan(sum([R G1 G2],2)));
    R=R(nonmissingobs,:);
    G1=G1(nonmissingobs,:);
    G2=G2(nonmissingobs,:);
    R1=R1(nonmissingobs,:);
    R2=R2(nonmissingobs,:);
    
    G1_baseline=G1;
    
    % Re are excess returns
    Re=R-rf(nonmissingobs);
    
    % Mimicking portfolio
    resmp1=MP(R1',G1',pmax,nwlags); % MP w/ RmRf
    resmp2=MP(R2',G2',pmax,nwlags); % MP w/ FF3
    
    % Uses returns, case with zero-beta rate
    res0=PCAFM(R',G1(:,1:nf)',pmax,nwlags); % No controls, for FM
    res1=PCAFM(R',G1',pmax,nwlags); % All assets, w/ market
    res2=PCAFM(R',G2',pmax,nwlags); % All assets, w/ ff3
    
    % Uses excess returns, without zero-beta rate
    res0_nozero=PCAFM(Re',G1(:,1:nf)',pmax,nwlags); % No controls, for FM
    res1_nozero=PCAFM(Re',G1',pmax,nwlags); % All assets, w/ market
    res2_nozero=PCAFM(Re',G2',pmax,nwlags); % All assets, w/ ff3
    res1_nozero_ridge=PCAFM(Re',G1',7,nwlags); % Ridge
    
    % Run with weighted PCA
    res_weighted_2=PCA_CK(Re',G2',pmax,nwlags,1); % All assets, w/ ff3
    
    % Avg return
    meanret=mean(G1(:,1:nf));
    
    % Sparsity test and figure
    p_baseline=7;
    V=res1_nozero.vres;
    V=squeeze(V(:,:,p_baseline));
    C=V./sqrt(diag(V)*diag(V)');
    C(abs(C)<0.25)=0;
    
    V1=res1_nozero.vres_nofactor;
    C1=V1./sqrt(diag(V1)*diag(V1)');
    C1(abs(C1)<0.25)=0;
    
    if j==1
        figure
        subplot(1,2,1);
        spy(C1)
        title('Sparsity plot, returns')
        xlabel('Assets');
        ylabel('Assets');
        subplot(1,2,2);
        spy(C)
        title('Sparsity plot, residuals')
        xlabel('Assets');
        ylabel('Assets');
        set(gcf, 'units', 'inches', 'position', [2 2 12 7])
        set(gcf, 'Color', 'w');
        save2pdf('output/sparsity')
        close
    end
    
    
    % -----------------------------------------
    % Results using equity only
    % -----------------------------------------
    
    userets={'wrds','ff25','ind','op_inv','me_var','me_ni','me_beta','me_ac','me_mom'};
    
    % Extract matrices
    R=[];
    for i=1:length(userets)
        R=[R returns.(userets{i})];
    end
    R=R(:,~isnan(sum(R)));
    
    G1=factors.(factorlist{j});
    
    rf=factors.rf;
    
    % Remove missing
    nonmissingobs=find(~isnan(sum([R G1],2)));
    R=R(nonmissingobs,:);
    G1=G1(nonmissingobs,:);
    
    % Re are excess returns
    Re=R-rf(nonmissingobs);
    
    % Uses returns, case with zero-beta rate
    res1_equityonly=PCAFM(R',G1',pmax,nwlags); 
    R_equityonly=R;
    G1_equityonly=G1;
    
    % Uses excess returns, without zero-beta rate
    res1_nozero_equityonly=PCAFM(Re',G1',pmax,nwlags); 
    
    
    % ------------------------------
    % Tables
    % ------------------------------
    
    p_baseline=7;
    p_robustness1=10;
    p_robustness2=13;
    
    % Table 1: Results without zero-beta rate
    % --------------------------------------------

    block=cell(0);
    for jj=1:nf
        
        row=cell(0);
        
        row=[row lab{j}{jj} ' '];
        
        % Average return
        if sum(strncmp(factorlist{j},tradedfactors,20))
            row=[row cellblock(meanret(jj)*100,sqrt(res1_nozero.avarghat(jj))*100,2) ' '];
        else
            row=[row ' ' ' ' ' ' ' '];
        end
    
        % FM no controls
        row=[row cellblock(res0_nozero.GammaFM_nozero(jj)*100,sqrt(res0_nozero.avarhatFM_nozero(jj))*100,2) ' '];
    
        % FM w/ rmrf
        row=[row cellblock(res1_nozero.GammaFM_nozero(jj)*100,sqrt(res1_nozero.avarhatFM_nozero(jj))*100,2) ' '];
    
        % FM w/ FF3
        row=[row cellblock(res2_nozero.GammaFM_nozero(jj)*100,sqrt(res2_nozero.avarhatFM_nozero(jj))*100,2) ' '];
    
        % MP w/ Rmrf
        row=[row cellblock(resmp1.gtMP(jj)*100,sqrt(resmp1.avarMP(jj))*100,2) ' '];
    
        % MP w/ FF3
        row=[row cellblock(resmp2.gtMP(jj)*100,sqrt(resmp2.avarMP(jj))*100,2) ' '];
    
        % FAFM, pp=p_baseline
        pp=p_baseline;
        row=[row cellblock(res1_nozero.Gammahat_nozero(jj+1,pp)*100,sqrt(res1_nozero.avarhat_nozero(jj,pp))*100,2) ' '];
    
        R2g=['$' num2str(res2.R2G(jj,pp)*100,'%2.2f') '$'];
        row=[row R2g ' '];
        
        ctest=['$' num2str(res1.chi2test(jj,6),'%2.2f') '$'];
        row=[row ctest];
        
        block=[block;row];
        
        res2.R2F(pp)
        
    end
        
    TABLE1=[TABLE1;block;cell(1,size(block,2))];  
    TABLE1;
     
    
    % Appendix Table III.9
    % ----------------------------------------------------------

    block=cell(0);
    for jj=1:nf
        
        row=cell(0);
        
        row=[row lab{j}{jj} ' '];
        
        % pp=p_baseline
        pp=p_baseline;
        row=[row cellblock(res1_nozero.Gammahat_nozero(jj+1,pp)*100,sqrt(res1_nozero.avarhat_nozero(jj,pp))*100,2) ' '];
        
        % p_robustness1
        pp=p_robustness1;
        row=[row cellblock(res1_nozero.Gammahat_nozero(jj+1,pp)*100,sqrt(res1_nozero.avarhat_nozero(jj,pp))*100,2) ' '];
    
        % p_robustness2
        pp=p_robustness2;
        row=[row cellblock(res1_nozero.Gammahat_nozero(jj+1,pp)*100,sqrt(res1_nozero.avarhat_nozero(jj,pp))*100,2) ' '];
    
        % pp=p_baseline, zero-beta rate
        pp=p_baseline;
        row=[row cellblock(res1.Gammahat(jj+1,pp)*100,sqrt(res1.avarhat(jj+1,pp))*100,2) ' '];
        
        % CK
        row=[row num2str(res_weighted_2.Gammahat(jj+1,p_baseline)*100,'%2.2f') ' '];
        
        % Ridge
        row=[row num2str(res1_nozero_ridge.Gammahat_nozero_ridge(jj)*100,'%2.2f') ' '];
        
        % Only equity
        pp=p_baseline;
        row=[row cellblock(res1_nozero_equityonly.Gammahat_nozero(jj+1,pp)*100,sqrt(res1_nozero_equityonly.avarhat_nozero(jj,pp))*100,2)];
    
        block=[block;row];
        
    end
        
    TABLE2=[TABLE2;block;cell(1,size(block,2))]; 
    
 
    % Plot of gt vs gthat
    % ----------------------------
    Gtemp=G1_baseline-mean(G1_baseline); % Demean original factor
    
    for jx=1:nf
        
        if strcmp(factorlist{j},'rmrf')
            figure;
            title('RmRf');
        elseif strcmp(factorlist{j},'smb')
            figure;
            title('SMB');
        elseif strcmp(factorlist{j},'ipgrowth')
            figure;
            title('IP growth');
        elseif strcmp(factorlist{j},'hml')
            figure;
            title('HML');
        else
            close
            continue;
        end
        
        rclean=res1_nozero.gthat(jx,:,p_baseline)';
        rorig=Gtemp(:,jx);
        
        hold on
        plot(cumsum(rorig),'k','LineWidth',1);
        plot(cumsum(rclean),'k','LineWidth',3);
        legend('Original factor','Cleaned factor','location','northwest')
        plotaxis(masterdates,5)
        ylabel('Cumulative sum of factor')
        
        set(gca,'FontSize',28);
        set(gcf, 'units', 'inches', 'position', [2 2 10 8])
        set(gcf, 'Color', 'w');
        save2pdf(['output/ghat_' num2str(j) '_' num2str(jx)]);
        close
  
    end

end 


% ------------------------------
% Robustness, cross-sectional
% ------------------------------
p_baseline=7;
nwlags=3;

userets={'wrds','br_bonds','br_fx','ff25','ind','op_inv','me_var','me_ni','me_beta','me_ac','me_mom'};

factorlist={'rmrf','ff3','mom','ipgrowth','liquidity','intermediary'};

Rall=[];
for i=1:length(userets)
    Rall=[Rall returns.(userets{i})];
end
Rall=Rall(:,~isnan(sum(Rall)));

rf=factors.rf;

NBOOT=10000;

for n=1:NBOOT
    
    if mod(n,100)==0
        n
    end
    
    % Select half
    FF=randsample(size(Rall,2),floor(size(Rall,2)/2),false);
    
    for j=1:length(factorlist)
        
        factorlist{j};
        
        G=factors.(factorlist{j});
        R=Rall(:,FF);
    
        % Remove missing
        nonmissingobs=find(~isnan(sum([R G],2)));
        R=R(nonmissingobs,:);
        G=G(nonmissingobs,:);
        
        Re=R-rf(nonmissingobs);

        res=PCAFM_fast(Re',G',p_baseline,nwlags);
        rp.(factorlist{j})(:,n)=res(2:end,2);
        
    end
      
end

figure    
hist(rp.rmrf',50)    
xlim([-0.002 0.01])
legend('RmRf')
set(gcf, 'units', 'inches', 'position', [2 2 12 7])
set(gcf, 'Color', 'w');
save2pdf('output/hist_cs_rmrf')
close 

figure    
hist(rp.ff3(2:3,:)',20)    
xlim([-0.002 0.01])
legend('SMB','HML')
set(gcf, 'units', 'inches', 'position', [2 2 12 7])
set(gcf, 'Color', 'w');
save2pdf('output/hist_cs_smbhml')
close 

figure    
hist(rp.mom',50)    
xlim([-0.002 0.01])
legend('MOM')
set(gcf, 'units', 'inches', 'position', [2 2 12 7])
set(gcf, 'Color', 'w');
save2pdf('output/hist_cs_mom')
close

    
figure    
hist(rp.ipgrowth',50)    
xlim([-0.002 0.01])
legend('ip growth')
set(gcf, 'units', 'inches', 'position', [2 2 12 7])
set(gcf, 'Color', 'w');
save2pdf('output/hist_cs_ip')
close 

figure    
hist(rp.liquidity',50)    
xlim([-0.002 0.01])
legend('Liquidity')
set(gcf, 'units', 'inches', 'position', [2 2 12 7])
set(gcf, 'Color', 'w');
save2pdf('output/hist_cs_liq')
close 

figure    
hist(rp.intermediary',50)    
xlim([-0.002 0.01])
legend('Intermediary (He et al.)','Intermediary (Adrian et al.)','location','northwest')
set(gcf, 'units', 'inches', 'position', [2 2 12 7])
set(gcf, 'Color', 'w');
save2pdf('output/hist_cs_interm')
close


% ------------------------------
% Robustness, time series
% ------------------------------

p_baseline=7;
nwlags=3;

userets={'wrds','br_bonds','br_fx','ff25','ind','op_inv','me_var','me_ni','me_beta','me_ac','me_mom'};

factorlist={'rmrf','ff3','mom','ipgrowth','liquidity','intermediary'};

Rall=[];
for i=1:length(userets)
    Rall=[Rall returns.(userets{i})];
end
Rall=Rall(:,~isnan(sum(Rall)));

rfall=factors.rf;

NBOOT=10000;

for n=1:NBOOT
    
    if mod(n,100)==0
         n
     end
    
    for j=1:length(factorlist)
        
        factorlist{j};
        
        G=factors.(factorlist{j});
        R=Rall;
        rf=rfall;
    
        % Remove missing
        nonmissingobs=find(~isnan(sum([R G],2)));
        R=R(nonmissingobs,:);
        G=G(nonmissingobs,:);
        rf=rf(nonmissingobs,:);
        
        % Resample
        select=randsample(size(R,1),floor(size(R,1)/2),false);
        
        R=R(select,:);
        G=G(select,:);
        rf=rf(select,:);
        
        Re=R-rf;
        
        res=PCAFM_fast(Re',G',p_baseline);
        rp.(factorlist{j})(:,n)=res(2:end,2);
        
    end
      
end

figure    
hist(rp.rmrf',50)    
xlim([-0.002 0.01])
legend('RmRf')
set(gcf, 'units', 'inches', 'position', [2 2 12 7])
set(gcf, 'Color', 'w');
save2pdf('output/hist_ts_rmrf')
close 

figure    
hist(rp.ff3(2:3,:)',20)    
xlim([-0.002 0.01])
legend('SMB','HML')
set(gcf, 'units', 'inches', 'position', [2 2 12 7])
set(gcf, 'Color', 'w');
save2pdf('output/hist_ts_smbhml')
close 

figure    
hist(rp.mom',50)    
xlim([-0.002 0.01])
legend('MOM')
set(gcf, 'units', 'inches', 'position', [2 2 12 7])
set(gcf, 'Color', 'w');
save2pdf('output/hist_ts_mom')
close

    
figure    
hist(rp.ipgrowth',50)    
xlim([-0.002 0.01])
legend('ip growth')
set(gcf, 'units', 'inches', 'position', [2 2 12 7])
set(gcf, 'Color', 'w');
save2pdf('output/hist_ts_ip')
close 

figure    
hist(rp.liquidity',50)    
xlim([-0.002 0.01])
legend('Liquidity')
set(gcf, 'units', 'inches', 'position', [2 2 12 7])
set(gcf, 'Color', 'w');
save2pdf('output/hist_ts_liq')
close 

figure    
hist(rp.intermediary',50)    
xlim([-0.002 0.01])
legend('Intermediary (He et al.)','Intermediary (Adrian et al.)','location','northwest')
set(gcf, 'units', 'inches', 'position', [2 2 12 7])
set(gcf, 'Color', 'w');
save2pdf('output/hist_ts_interm')
close


