%% This Version: Nov 21, 2020. @copyright Stefano Giglio and Dacheng Xiu
%% If you use these codes, please cite the paper "Asset Pricing with Omitted Factors." (2020)
%% Main Simulation Codes for Many T and N Combinations

clear;
close all;
warning('off');
poolobj = gcp('nocreate');
 if isempty(poolobj)
     parpool('local',4)
 end
 


Tlist = [120,240,480];
Nlist = [50,100,200];

for tt = 1:length(Tlist)
    
    for nn = 1:length(Nlist)

M    = 1000;%5000;  % # of MCs
T    = Tlist(tt);   % # of time periods
n    = Nlist(nn);   % # of stocks
p    = 5;     % # of latent factors
d    = 4;     % # of factor proxies
pmax = 10;     % # 1:pmax is the list of number of factors
q    = floor(T^(1/4))+1;


load('Calibrated_Parameters.mat');

%%
mu           = zeros(p,1);     
alpha0       = 0;
Gammahat     = zeros(d+1,pmax,M);
Gammahat_nozero= zeros(d+1,pmax,M);
Gammahat_nozero_ridge = zeros(d,M);
GammahatLP     = zeros(d+1,pmax,M);
GammahatCK     = zeros(d+1,pmax,M);
GammaFM      = zeros(d+1,M);
GammaTrueFM  = zeros(d+1,M);
GammaTrueFM_nozero  = zeros(d+1,M);
gtMP = zeros(d,M);
gtavarMP= zeros(d,M);

GammahatOL   = zeros(d+1,M);
GammaFM_Crt  = zeros(d+1,M);
avarhat      = zeros(d+1,pmax,M);
avarhat_nozero      = zeros(d,pmax,M);
avarhatFM    = zeros(d+1,M);

gthat        = zeros(d,T,pmax,M);
R2F          = zeros(pmax,M);
R2G          = zeros(d,pmax,M);
phat         = zeros(3,M);
gtavar       = zeros(d,T,M);
gttrue       = zeros(d,T,M);
chi2test     = zeros(d,pmax,M);
beta         = repmat(beta0',n,1) + randn(n,p)*(Sigmabeta-beta0*beta0')^0.5; % factor loading

dlist        = divisors(n);
idx          = find(abs(dlist-floor(sqrt(n)))==min(abs(dlist-floor(sqrt(n)))),1);
Sigmau       = Sigmau(1:n,1:n);

xi           = zeros(d,1);
%gamma0       = 0; % UNCOMMENT THIS LINE IF NO ZERO BETA RATE IN THE TRUE DGP!

Gammatrue    = [gamma0;eta*gamma];       % zero-beta rate + true risk premiums of the proxies

R2FTrue      = gamma'*(Sigmabeta - beta0*beta0')*gamma/(gamma'*(Sigmabeta - beta0*beta0')*gamma + sigmaalpha^2);
R2GTrue      = diag(eta*Sigmav*eta')./diag(eta*Sigmav*eta'+Sigmaw); 


Sigma      = beta*Sigmav*beta'+Sigmau; 
lambdastar  = eta*Sigmav*beta'*inv(Sigma)*(beta*gamma);
parfor iMC = 1:M

  
    disp(iMC);
    alpha = (alpha0 +randn(n,1)*sigmaalpha);
    vt    = Sigmav^0.5*randn(p,T);          % factor innovations
    ut    = Sigmau^0.5*randn(n,T);          % residual innovations

    wt    = Sigmaw(1:d,1:d)^0.5*randn(d,T);          % proxies residual innovations
    rt    = repmat(ones(n,1)*gamma0+alpha,1,T) + repmat(beta*gamma,1,T) + beta*vt + ut; % returns
    gt    = repmat(xi,1,T) + eta*vt + wt;
    gttrue(:,:,iMC)= eta*vt; 
  
    %% Estimation
        
    res               = PCAFM(rt,gt,pmax,q);
    resFMvsMP         = FM(rt,vt,gt,q);
   

    GammaFM(:,iMC)    = res.GammaFM;    
    Gammahat(:,:,iMC) = res.Gammahat;
    Gammahat_nozero(:,:,iMC) = res.Gammahat_nozero;

    GammaTrueFM(:,iMC) = resFMvsMP.Gammahat;
    GammaTrueFM_nozero(:,iMC) = resFMvsMP.Gammahat_nozero;
    
    gtMP(:,iMC) = res.gtMP;
    gtavarMP(:,iMC) = res.avarMP;
    
    avarhat(:,:,iMC)  = res.avarhat;
    avarhat_nozero(:,:,iMC)  = res.avarhat_nozero;
    avarhatFM(:,iMC)  = res.avarhatFM;
    gthat(:,:,:,iMC)  = res.gthat;
    R2F(:,iMC)        = res.R2F;
    R2G(:,:,iMC)      = res.R2G;
    phat(:,iMC)       = res.phat;

    Gammahat_nozero_ridge(:,iMC)    = res.Gammahat_nozero_ridge; 

end


save(['Simulation_Results_T_',num2str(T),'_N_',num2str(n),'.mat'],'GammaFM','Gammahat','Gammahat_nozero','gtMP','gtavarMP','avarhat','avarhat_nozero','avarhatFM','gthat',...
    'R2F','R2G','phat','Gammatrue','d','p','M','T','Gammahat_nozero_ridge');
    end
end
%% 
