clear;clc
%% For Table 2 and 3

st_date = 197401; et_date = 201612;
z = xlsread('Portfolio_ret.xlsx','FF5','A2:I653','basic'); 
tt = find(z(:,1)>=st_date & z(:,1)<=et_date);
rf = z(tt,7)/100;
mkt = z(tt,2)/100;
tdate = z(tt,1);
yy = floor(tdate/100); mm = floor(tdate - yy*100);

z = xlsread('Portfolio_ret.xlsx','Anomaly','A3:BS518','basic'); 
tt = find(z(:,1)>=st_date & z(:,1)<=et_date);
f5 = z(tt,2:7)/100; % benchmark: 1, 3, 5, 6 factors
G = z(tt,2:end)/100;
[T,L] = size(G);


%% 48 Industry 
z = xlsread('Portfolio_ret.xlsx','FF48vw','A3:AW530','basic'); 
tt = find(z(:,1)>=st_date & z(:,1)<=et_date);
N = size(z,2)-1;
R48 = z(tt,2:end)/100 - rf*ones(1,N);
 

%% 202 target
z = xlsread('Portfolio_ret.xlsx','Dacheng202vw','A3:GU554','basic'); 
tt = find(z(:,1)>=st_date & z(:,1)<=et_date);
N = size(z,2)-1;
R202 = z(tt,2:end)/100 - rf*ones(1,N);

%% all stocks
load stock_ret.mat;
ret = st_dataset; clear st_dataset;
ret = double(ret); % permno, year, month, xret, me, exchcd
ret_all = [ret(:,1), ret(:,2)*100+ret(:,3), ret(:,4)]; 

%% all-but-micro stocks
load stock_ret_nomicro.mat;
ret = st_dataset; clear st_dataset;
ret = double(ret); % permno, year, month, xret, me, exchcd
ret_but = [ret(:,1), ret(:,2)*100+ret(:,3), ret(:,4)]; 

%% Tests
NF = [1 3 5 6];
NoFF = 6;

MPE = nan(NoFF,size(NF,2),4);  TR2 = nan(NoFF,size(NF,2),4);  
SR = nan(NoFF,size(NF,2),4);   A2R = nan(NoFF,size(NF,2),4);   

for j =1:size(NF,2); % No. of factors
mu = mean(G);
y1 = G'*G/T-mu'*mu;
[E,v] = eig(y1);
[v,ind] = sort(diag(v),'descend');
E = E(:,ind(1:NF(j))); 
pcaf = G*E*inv(E'*E);                     %PCA

 
[~,~,~,~,~,~,~,stats] = plsregress(G,R48,NF(j));% PLS:R48 are basis assets
plsf48 = G*stats.W; 
 
rra48  = func_rraff(R48,G,NF(j)); % RRA: R48 are basis assets
 
[~,~,~,~,~,~,~,stats] = plsregress(G,R202,NF(j));% PlS: R202 are basis assets
plsf202 = G*stats.W; 
 
rra202  = func_rraff(R202,G,NF(j)); % RRA: R202 are basis assets
 
 FFs = nan(T,NF(j),NoFF);
 FFs(:,:,1)  = f5(:,1:NF(j));
 FFs(:,:,2)  = pcaf;  
 FFs(:,:,3)= plsf48;  
 FFs(:,:,4)= rra48;
 FFs(:,:,5)= plsf202;  
 FFs(:,:,6)= rra202;

 for i = 1:NoFF;  % No. of models
           
      ff = FFs(:,:,i);  
      
      res  = performance_portfolio(R48, ff); % OUTres = [MPE, TotalR2, SR, A2R];
      MPE(i,j,1) = res(1); TR2(i,j,1) = res(2); 
      SR(i,j,1)=res(3); A2R(i,j,1) = res(4);
      
      res = performance_portfolio(R202, ff);
      MPE(i,j,2)=res(1); TR2(i,j,2)= res(2); 
      SR(i,j,2) =res(3); A2R(i,j,2)= res(4);
      
      res = performance_stock(tdate, ret_all, ff);
      MPE(i,j,3)=res(1);  TR2(i,j,3)= res(2); 
      SR(i,j,3) =res(3); A2R(i,j,3)= res(4);
      
      res = performance_stock(tdate, ret_but, ff);
      MPE(i,j,4)=res(1);  TR2(i,j,4)= res(2); 
      SR(i,j,4) =res(3); A2R(i,j,4)= res(4);
   
    end
    clear ff FFs;
end

%% Table 2: Performance of different factor models: In-sample
Panel_A = [MPE(1:4,:,1),TR2(1:4,:,1)]
Panel_B = [MPE(1:4,:,2),TR2(1:4,:,2)]
Panel_C = [MPE(1:4,:,3),TR2(1:4,:,3)]
Panel_D = [MPE(1:4,:,4),TR2(1:4,:,4)]

Tab2 = [Panel_A;nan(1,8);...
    Panel_B;nan(1,8);...
    Panel_C;nan(1,8);...
    Panel_D];

%% Table 3 Extracting PLS and RRA factors using alternative basis assets
Panel_A = [MPE(5:6,:,1),TR2(5:6,:,1)]
Panel_B = [MPE(5:6,:,2),TR2(5:6,:,2)]
Panel_C = [MPE(5:6,:,3),TR2(5:6,:,3)]
Panel_D = [MPE(5:6,:,4),TR2(5:6,:,4)]

Tab3 = [Panel_A;nan(1,8);...
    Panel_B;nan(1,8);...
    Panel_C;nan(1,8);...
    Panel_D];

