% Select factors following FGK(2020), the first-step Lasso.

clear; clc;
close all;

addpath('./FGX2020')
addpath('./FGX2020/glmnet_matlab')
addpath('./FGX2020/functions')


%% options

% sample period
st_date = 197401; 
et_date = 201612;

% oos period
L=360;
T_os = st_date+L/12*100;

% test asset
test_asset = 48; % ==48, will use FF 48portfolios as test assets; == 202, will use 202 portfolios as test assets.

% number of L1 grinds
NumLambda = 150;

% run DS Lasso
ds = 1;

% K fold
Kfld = 10;

% repetition under different random seeds
Jrep=200;

% alpha
alpha=1; % 1 means Lasso

% fix the random seed for any additional CV
seednum = 100;


%% import data


% riskfree rate
z = xlsread('Portfolio_ret.xlsx','FF5','A2:I653','basic'); 
tt = find(z(:,1)>=st_date & z(:,1)<=et_date);
rf = z(tt,7)/100; 
tdate = z(tt,1);

% anomalies
[z, factorname]= xlsread('Portfolio_ret.xlsx','Anomaly','A2:BS518'); 
tt = find(z(:,1)>=st_date & z(:,1)<=et_date);
anomalies = z(tt,2:end)/100;
anomalies = anomalies';

tt = find(z(:,1)>=st_date & z(:,1)<T_os);
anomalies_is = z(tt,2:end)/100;
anomalies_is = anomalies_is';

tt = find(z(:,1)>=T_os & z(:,1)<=et_date);
anomalies_os = z(tt,2:end)/100;
anomalies_os = anomalies_os';

P = size(anomalies,1)
factorname = factorname(2:end)';


% test assets
if test_asset==48
    z = xlsread('Portfolio_ret.xlsx','FF48vw','A3:AW530','basic');  % 48 Industries
    tt = find(z(:,1)>=st_date & z(:,1)<=et_date);
    N = size(z,2)-1;
    Ri = z(tt,2:end)/100 - rf*ones(1,N); %excess returns of FF48
    Ri = Ri'; % test asset, N-by-T

else
    z = xlsread('Portfolio_ret.xlsx','Dacheng202vw','A3:GU554','basic'); 
    tt = find(z(:,1)>=st_date & z(:,1)<=et_date);
    N = size(z,2)-1;
    Ri = z(tt,2:end)/100 - rf*ones(1,N); %excess returns of 202 chr portfolios
    Ri = Ri'; % test asset, N-by-T
end    

assert(size(Ri,2)==size(anomalies,2),"Check the sample periods of each return data.")


%% Standard Lasso (Step 1)
model_ss = DS_TSCV(Ri, anomalies(1,:), anomalies, alpha,NumLambda,seednum, Kfld, Jrep, 0);
num_L1factors =  sum(model_ss.coeff_path~=0); % number of non-zero coeffs
NF=[3,5,6];
step1_factors = cell(1,size(NF,2)+1); % save factors corresponding to [3,5,6,CV], in cell

for i=1:size(NF,2);   
    I =  num_L1factors==(NF(i)+1);   
    I = find(model_ss.cvm111 == min(model_ss.cvm111(I)), 1, 'last');
    L1coeff = model_ss.coeff_path(:,I);
    factorname(L1coeff(2:end)~=0) % convert the index to factornames, 
    step1_factors{i} = [tdate anomalies(L1coeff(2:end)~=0,:)'];
end   

factorname(model_ss.sel1) % convert the index to factornames, 
step1_factors{end} = [tdate anomalies(model_ss.sel1,:)']; %CV
if test_asset == 48
    save('FGX_factors.mat','step1_factors');
else 
    save('FGX_factors_IA.mat','step1_factors');
end


