% The notation follows Simmet and Pohlmeier (2018)
% This version 02/2019

% The following code performes the Monte-Carlo simulation  
% of artifitial market portfolio excess returns, as well as 
% different proxy market returns with random and fixed
% measurement errors, and estimates the CAPM parameters
% for all the cases. 

clc;clear;

rng(5231); % random seed for replication of the Monte-Carlo results
%% Loading necessary data for the Monte Carlo simulation parameters
% Generation of excess returns
load sp500_edata.mat % price data for the constituents of the S&P500 from DataStream and 
                                 % market returns from French & Fama database 
%save('sp500price.txt',  '-ascii')                                 
sp500_stocks=table2array(sp500_stocks);            % price data for the constituents of the S&P500 from DataStream without columns with missing values
%save('sp500_stocks.txt',  '-ascii')
FactorData=table2array(FactorData);                     % market returns, 3 factors and risk free return from French & Fama database 
%save('FactorData.txt',  '-ascii')
RF=FactorData(end-495:end,5)/100;                      % risk free returns from French & Fama database 
returns_simple=(sp500_stocks(2:end,:)-sp500_stocks(1:end-1,:))./sp500_stocks(1:end-1,:); % simple returns
returns_excess=returns_simple-repmat(RF,1,205); % excess returns
r=returns_excess;
r=r*100;                                                               % returns in percentage
                     
% Generation of the vector of true weights for  artifitial market portfolio
load MV.mat                   % total market capitalization (01.01.1974-01.05.2015 monthly data) of 205 stocks of S&P500
%save('MV.txt',  '-ascii')
MVwoNA=MV(30:end,:);   % delete rows with missing data (01.06.1976-01.05.2015 monthly data)
meanbased=mean(MVwoNA,1)/sum(mean(MVwoNA,1)); %create sample supply vector (qs) that is the average share of the total market capitalization
qs=meanbased';                 



%% Parameters for the simulation

N = 205;                % total number of assets in the asset universe
T = 120;                % number of observations corresponding to 10 years of monthly data
CS=15;                 % dimension of the estimated CAPM system
Rep=10000;           % number of Monte Carlo replications

%% Additional parameters for the data generation process

bstar=qs;                      % true weighting scheme of the market portfolio b* 
gamma=0.4;                % absolute risk aversion parameter of the artifitial market
Omega=cov(r);             % true variance covariance matrix of the whole asset space
% randomly drawn assets for the CAPM regression estimation
INDC=datasample(1:N,CS,'Replace',false)'; 

% mean values for true market return and individual asset returns
 mu=repmat(gamma*Omega*bstar,1,T);                  % expected value of the asset returns from the equilibrium formula
 mu_m_star=bstar'*mu;                                         % expected value of the true market return 
 beta_formula=Omega*bstar/(bstar'*Omega*bstar);  % the vector of CAPM betas according to the formula (9)

%% Misspecified weighting schemes
% Weights for equally weighted market portfolio 
b.ew=ones(N,1)/N; 
% Standard deviaion of the measurement error for the weights with random error
sigma_nu=[0.01 0.025 0.05];
% Obtaining weights with fixed error
% First, we sort the assets in the asset universe 
% according to their true weight
[sortedValue_X , X_Ranked] = sort(bstar,'descend'); % returns sorted weights and the corresponding index of the asset
cumsumqs=cumsum(sortedValue_X);                     % returns the vector of the cummulative weight of the first largest assets
M=[1:1:204];                                                         % number of assets included in the proxy index

b.I2=zeros(N,N-1);                                                 % Space holder for standardized weights for a subset of assets
b.I3=zeros(N,N-1);                                                 % Space holder for equally weighted scheme for a subset of assets
 
for i=1:204
  ind=X_Ranked(1:M(i),1);
  b.I2(ind,i)=bstar(ind,1)/sum(bstar(ind,1));               % Weights for the index consisting of i assets; can be changed to b.I2(ind,i)=bstar(ind,1) to get results for index without scaling
  b.I3(ind,i)=1/M(i);                                                % Weights for the equally-weighted index consisting of i assets
end

%% Memory allocation 
error=zeros(Rep,N,T);      % simulated sampling error 

% coefficients' estimates from the CAPM models with different market proxies
coeff.true=zeros(Rep,CS,2);     coeff.ew=zeros(Rep,CS,2);    coeff.re01=zeros(Rep,CS,2);
coeff.re025=zeros(Rep,CS,2);   coeff.re05=zeros(Rep,CS,2); 
coeff.I2=zeros(204,Rep,CS,2);  coeff.I3=zeros(204,Rep,CS,2);

% standard deviations of the estimated coefficients
std.true=zeros(Rep,CS,2);     std.ew=zeros(Rep,CS,2); std.re01=zeros(Rep,CS,2);
std.re025=zeros(Rep,CS,2); std.re05=zeros(Rep,CS,2);
std.I2=zeros(204,Rep,CS,2); std.I3=zeros(204,Rep,CS,2);

% t-statistics for the test of the significance of the regression coefficient
tstat.true=zeros(Rep,CS,2);    tstat.ew=zeros(Rep,CS,2); tstat.re01=zeros(Rep,CS,2); 
tstat.re025=zeros(Rep,CS,2); tstat.re05=zeros(Rep,CS,2); 
tstat.I2=zeros(204,Rep,CS,2); tstat.I3=zeros(204,Rep,CS,2);

% p-value  for the test of the significance of the regression coefficient
pval.true=zeros(Rep,CS,2);    pval.ew=zeros(Rep,CS,2); pval.re01=zeros(Rep,CS,2);
pval.re025=zeros(Rep,CS,2); pval.re05=zeros(Rep,CS,2);
pval.I2=zeros(204,Rep,CS,2); pval.I3=zeros(204,Rep,CS,2);

% Wald-test rejection for a joint hypothesis of absence of abnormal returns with the true VC matrix,
%  F-test rejection for a joint hypothesis for the absence of abnormal returns with the estimated VC matrix
% and the value of F-test statistic
wald_F_stat.true=zeros(Rep,3);  wald_F_stat.ew=zeros(Rep,3); wald_F_stat.re01=zeros(Rep,3);
wald_F_stat.re025=zeros(Rep,3); wald_F_stat.re05=zeros(Rep,3);
wald_F_stat.I2=zeros(204,Rep,3); wald_F_stat.I3=zeros(204,Rep,3);

% MSE of the estimated coefficients
MSE.true=zeros(Rep,CS,2);   MSE.ew=zeros(Rep,CS,2); MSE.re01=zeros(Rep,CS,2); 
MSE.re025=zeros(Rep,CS,2); MSE.re05=zeros(Rep,CS,2);
MSE.I2=zeros(204,Rep,CS,2); MSE.I3=zeros(204,Rep,CS,2); 

% Correlation coefficient between the proxy return and 
% the true market return
cor.ew=zeros(Rep,1); cor.re01=zeros(Rep,1);
cor.re025=zeros(Rep,1); cor.re05=zeros(Rep,1);
cor.I2=zeros(204,Rep,1); cor.I3=zeros(204,Rep,1);

% Slope coefficient, lambda_1, from the linear projection of the true market return on the
% proxy return
lambda.ew=zeros(Rep,1); lambda.re01=zeros(Rep,1);
lambda.re025=zeros(Rep,1); lambda.re05=zeros(Rep,1);
lambda.I2=zeros(204,Rep,1); lambda.I3=zeros(204,Rep,1);

% Reliability ratio
RR.ew=zeros(Rep,1); RR.re01=zeros(Rep,1);
RR.re025=zeros(Rep,1); RR.re05=zeros(Rep,1);
RR.I2=zeros(204,Rep,1); RR.I3=zeros(204,Rep,1);

% Performance measures
% Sharpe ratios
performance.SR_true=zeros(Rep,CS); performance.SR_ew=zeros(Rep,CS);performance.SR_re01=zeros(Rep,CS);
performance.SR_re025=zeros(Rep,CS); performance.SR_re05=zeros(Rep,CS);
performance.SR_I2=zeros(204,Rep,CS); performance.SR_I3=zeros(204,Rep,CS);

% Treynor's ratios
performance.TR_true=zeros(Rep,CS);  performance.TR_ew=zeros(Rep,CS); performance.TR_re01=zeros(Rep,CS);
performance.TR_re025=zeros(Rep,CS); performance.TR_re05=zeros(Rep,CS);
performance.TR_I2=zeros(204,Rep,CS); performance.TR_I3=zeros(204,Rep,CS);

% Estimated slope of the Sequrity Market Line
SML.true=zeros(Rep,1); SML.ew=zeros(Rep,1); SML.re01=zeros(Rep,1);
SML.re025=zeros(Rep,1); SML.re05=zeros(Rep,1);
SML.I2=zeros(204,Rep,1); SML.I3=zeros(204,Rep,1);

%% Simulation and estimation block
tic;
 
 for r=1:Rep
 % Generate normal random errors
 epsilon= mvnrnd(zeros(T,N),Omega)'; 
 error(r,:,:)=epsilon;
 % Generate assets' returns and the true market return
 simret=mu+epsilon;                      % individual assets' returns
 Rmstar=bstar'*simret;                   % true market return
 SML.true(r,:)=mean(Rmstar);         % estimated slope of the SML
  
 %% Estimation of the CAPM parameters for 15 random assets by SUR method with the true markt return
 % see function sure.m
 [tmp1_p,VC_U,tmp1_cv,~,~,tmp1_t,tmp1_w1,tmp1_w2]=sure(simret(INDC,:)',Rmstar',1,Omega(INDC,INDC)); % const=1 to add a column of 1s
 
 % Store the estimates
 coeff.true(r,:,:) = tmp1_p';                         % estimated coefficients
 std.true(r,:,:)=reshape(diag(tmp1_cv),2,[])'; % standard deviations of the estimates
 tstat.true(r,:,:)=reshape(tmp1_t,2,[])';         % t-statistics
 wald_F_stat.true(r,2)=(tmp1_w2>finv(0.95,CS,T-CS-1)); % rejection of the null hypothesis for the F test
 wald_F_stat.true(r,1)=(tmp1_w1>chi2inv(0.95,CS-1));    % rejection of the null hypothesis for the Wald test
 wald_F_stat.true(r,3)=(tmp1_w2); % F test statistic
 MSE.true(r,:,1)=(tmp1_p(1,:)'-zeros(CS,1)).^2;              %MSE alpha
 MSE.true(r,:,2)=(tmp1_p(2,:)'-beta_formula(INDC)).^2;   %MSE beta
 % performance measures
 performance.TR_true(r,:)=(mean(simret(INDC,:),2)./tmp1_p(2,:)')'; 
 performance.SR_true(r,:)=(mean(Rmstar)*tmp1_p(2,:)'./diag(cov(simret(INDC,:)')))';
 
 %% equally weighted
 R.ew=b.ew'*simret;
 SML.ew(r,:)=mean(R.ew);
 % reliability ratio, correlation and lambda
 RR.ew(r,:)=cov(Rmstar)/cov(R.ew);
 corrm=corrcoef(Rmstar,R.ew);
 cor.ew(r,:)=corrm(2,1);
 lambdas = inv([ones(T,1) R.ew']'*[ones(T,1) R.ew'])*[ones(T,1) R.ew']'*Rmstar';     
 lambda.ew(r,:) = lambdas(2,1);
 
  % CAPM regressions estimation
 [tmp1_p,VC_U,tmp1_cv,~,~,tmp1_t,tmp1_w1,tmp1_w2]=sure(simret(INDC,:)',R.ew',1,Omega(INDC,INDC)); % const=1 to add a column of 1s
 
 % Store data
 coeff.ew(r,:,:) = tmp1_p';
 std.ew(r,:,:)=reshape(diag(tmp1_cv),2,[])'; 
 tstat.ew(r,:,:)=reshape(tmp1_t,2,[])';  
 wald_F_stat.ew(r,2)=(tmp1_w2>finv(0.95,CS,T-CS-1));
 wald_F_stat.ew(r,1)=(tmp1_w1>chi2inv(0.95,CS-1)); 
 wald_F_stat.ew(r,3)=(tmp1_w2);
 MSE.ew(r,:,1)=(tmp1_p(1,:)'-zeros(CS,1)).^2;
 MSE.ew(r,:,2)=(tmp1_p(2,:)'-beta_formula(INDC)).^2;
 % performance measures
 performance.TR_ew(r,:)=(mean(simret(INDC,:),2)./tmp1_p(2,:)')'; 
 performance.SR_ew(r,:)=(mean(R.ew)*tmp1_p(2,:)'./diag(cov(simret(INDC,:)')))';
 
 %% random error with standard deviation 0.01
 b.re01=bstar+normrnd(0,sigma_nu(1),205,1);
 R.re01=b.re01'*simret;
 SML.re01(r,:)=mean(R.re01);
 % reliability ratio, correlation and lambda
 RR.re01(r,:)=cov(Rmstar)/cov(R.re01);
 corrm=corrcoef(Rmstar,R.re01);
 cor.re01(r,:)=corrm(2,1);
 lambdas = inv([ones(T,1) R.re01']'*[ones(T,1) R.re01'])*[ones(T,1) R.re01']'*Rmstar';     
 lambda.re01(r,:) = lambdas(2,1);
 
 
 % CAPM regressions estimation
 [tmp1_p,VC_U,tmp1_cv,~,~,tmp1_t,tmp1_w1,tmp1_w2]=sure(simret(INDC,:)',R.re01',1,Omega(INDC,INDC)); % const=1 to add a column of 1s

 % Store data
 coeff.re01(r,:,:) = tmp1_p';
 std.re01(r,:,:)=reshape(diag(tmp1_cv),2,[])'; 
 tstat.re01(r,:,:)=reshape(tmp1_t,2,[])';  
 wald_F_stat.re01(r,2)=(tmp1_w2> finv(0.95,CS,T-CS-1));
 wald_F_stat.re01(r,1)=(tmp1_w1>chi2inv(0.95,CS-1)); 
 wald_F_stat.re01(r,3)=(tmp1_w2);
 MSE.re01(r,:,1)=(tmp1_p(1,:)'-zeros(CS,1)).^2;
 MSE.re01(r,:,2)=(tmp1_p(2,:)'-beta_formula(INDC)).^2;
 % performance measures
 performance.TR_re01(r,:)=(mean(simret(INDC,:),2)./tmp1_p(2,:)')'; 
 performance.SR_re01(r,:)=(mean(R.re01)*tmp1_p(2,:)'./diag(cov(simret(INDC,:)')))';
 
 %% random error with standard deviation 0.025
 b.re025=bstar+normrnd(0,sigma_nu(2),205,1);
 R.re025=b.re025'*simret;
 SML.re025(r,:)=mean(R.re025);
 % reliability ratio, correlation and lambda
 RR.re025(r,:)=cov(Rmstar)/cov(R.re025);
 corrm=corrcoef(Rmstar,R.re025);
 cor.re025(r,:)=corrm(2,1);
 lambdas = inv([ones(T,1) R.re025']'*[ones(T,1) R.re025'])*[ones(T,1) R.re025']'*Rmstar';     
 lambda.re025(r,:) = lambdas(2,1);
 
 
 % CAPM regressions estimation
 [tmp1_p,VC_U,tmp1_cv,~,~,tmp1_t,tmp1_w1,tmp1_w2]=sure(simret(INDC,:)',R.re025',1,Omega(INDC,INDC)); % const=1 to add a column of 1s
  
 % Store data
 coeff.re025(r,:,:) = tmp1_p';
 std.re025(r,:,:)=reshape(diag(tmp1_cv),2,[])'; 
 tstat.re025(r,:,:)=reshape(tmp1_t,2,[])';  
 wald_F_stat.re025(r,2)=(tmp1_w2>finv(0.95,CS,T-CS-1));
 wald_F_stat.re025(r,1)=(tmp1_w1>chi2inv(0.95,CS-1)); 
 wald_F_stat.re025(r,3)=(tmp1_w2);
 MSE.re025(r,:,1)=(tmp1_p(1,:)'-zeros(CS,1)).^2;
 MSE.re025(r,:,2)=(tmp1_p(2,:)'-beta_formula(INDC)).^2;
 % performance measures
 performance.TR_re025(r,:)=(mean(simret(INDC,:),2)./tmp1_p(2,:)')'; 
 performance.SR_re025(r,:)=(mean(R.re025)*tmp1_p(2,:)'./diag(cov(simret(INDC,:)')))';
 
 %% random error  with standard deviation 0.05
 b.re05=bstar+normrnd(0,sigma_nu(3),205,1);
 R.re05=b.re05'*simret;
 SML.re05(r,:)=mean(R.re05);
 % reliability ratio, correlation and lambda
 RR.re05(r,:)=cov(Rmstar)/cov(R.re05);
 corrm=corrcoef(Rmstar,R.re05);
 cor.re05(r,:)=corrm(2,1);
 lambdas = inv([ones(T,1) R.re05']'*[ones(T,1) R.re05'])*[ones(T,1) R.re05']'*Rmstar';     
 lambda.re05(r,:) = lambdas(2,1);
 
 
 % CAPM regressions estimation
 [tmp1_p,VC_U,tmp1_cv,~,~,tmp1_t,tmp1_w1,tmp1_w2]=sure(simret(INDC,:)',R.re05',1,Omega(INDC,INDC)); % const=1 to add a column of 1s
 
 % Store data
 coeff.re05(r,:,:) = tmp1_p';
 std.re05(r,:,:)=reshape(diag(tmp1_cv),2,[])'; 
 tstat.re05(r,:,:)=reshape(tmp1_t,2,[])';  
 wald_F_stat.re05(r,2)=(tmp1_w2>finv(0.95,CS,T-CS-1));
 wald_F_stat.re05(r,1)=(tmp1_w1>chi2inv(0.95,CS-1)); 
  wald_F_stat.re05(r,3)=(tmp1_w2);
 MSE.re05(r,:,1)=(tmp1_p(1,:)'-zeros(CS,1)).^2;
 MSE.re05(r,:,2)=(tmp1_p(2,:)'-beta_formula(INDC)).^2;
 % performance measures
 performance.TR_re05(r,:)=(mean(simret(INDC,:),2)./tmp1_p(2,:)')'; 
 performance.SR_re05(r,:)=(mean(R.re05)*tmp1_p(2,:)'./diag(cov(simret(INDC,:)')))';
 
 % Estimation of the CAPM with the proxies subject to fixed measurement error
 for m=1:204
   %% Proxy based on the weighting scheme b_I2
 R.I2=b.I2(:,m)'*simret;
 SML.I2(m,r,:)=mean(R.I2);
 % reliability ratio, correlation and lambda
 RR.I2(m,r,:)=cov(Rmstar)/cov(R.I2);
 corrm=corrcoef(Rmstar,R.I2);
 cor.I2(m,r,:)=corrm(2,1);
 lambdas = inv([ones(T,1) R.I2']'*[ones(T,1) R.I2'])*[ones(T,1) R.I2']'*Rmstar';     
 lambda.I2(m,r,:) = lambdas(2,1);
 
 
 % CAPM regressions estimation
 [tmp1_p,VC_U,tmp1_cv,~,~,tmp1_t,tmp1_w1,tmp1_w2]=sure(simret(INDC,:)',R.I2',1,Omega(INDC,INDC)); % const=1 to add a column of 1s
 
 % Store data
 coeff.I2(m,r,:,:) = tmp1_p';
 std.I2(m,r,:,:)=reshape(diag(tmp1_cv),2,[])'; 
 tstat.I2(m,r,:,:)=reshape(tmp1_t,2,[])';  
 wald_F_stat.I2(m,r,2)=(tmp1_w2>finv(0.95,CS,T-CS-1));
 wald_F_stat.I2(m,r,1)=(tmp1_w1>chi2inv(0.95,CS-1)); 
 wald_F_stat.I2(m,r,3)=(tmp1_w2);
 MSE.I2(m,r,:,1)=(tmp1_p(1,:)'-zeros(CS,1)).^2;
 MSE.I2(m,r,:,2)=(tmp1_p(2,:)'-beta_formula(INDC)).^2;
 % performance measures
 performance.TR_I2(m,r,:)=(mean(simret(INDC,:),2)./tmp1_p(2,:)')'; 
 performance.SR_I2(m,r,:)=(mean(R.I2)*tmp1_p(2,:)'./diag(cov(simret(INDC,:)')))';
 
 %% Proxy based on the weighting scheme b_I3
 R.I3=b.I3(:,m)'*simret;
 SML.I3(m,r,:)=mean(R.I3);
 % reliability ratio, correlation and lambda
 RR.I3(m,r,:)=cov(Rmstar)/cov(R.I3);
 corrm=corrcoef(Rmstar,R.I3);
 cor.I3(m,r,:)=corrm(2,1);
 lambdas = inv([ones(T,1) R.I3']'*[ones(T,1) R.I3'])*[ones(T,1) R.I3']'*Rmstar';     
 lambda.I3(m,r,:) = lambdas(2,1);
 
 
 % CAPM regressions estimation
 [tmp1_p,VC_U,tmp1_cv,~,~,tmp1_t,tmp1_w1,tmp1_w2]=sure(simret(INDC,:)',R.I3',1,Omega(INDC,INDC)); % const=1 to add a column of 1s

 % Store data
 coeff.I3(m,r,:,:) = tmp1_p';
 std.I3(m,r,:,:)=reshape(diag(tmp1_cv),2,[])'; 
 tstat.I3(m,r,:,:)=reshape(tmp1_t,2,[])';  
 wald_F_stat.I3(m,r,2)=(tmp1_w2> finv(0.95,CS,T-CS-1));
 wald_F_stat.I3(m,r,1)=(tmp1_w1>chi2inv(0.95,CS-1));
 wald_F_stat.I3(m,r,3)=(tmp1_w2);
 MSE.I3(m,r,:,1)=(tmp1_p(1,:)'-zeros(CS,1)).^2;
 MSE.I3(m,r,:,2)=(tmp1_p(2,:)'-beta_formula(INDC)).^2;
 % performance measures
 performance.TR_I3(m,r,:)=(mean(simret(INDC,:),2)./tmp1_p(2,:)')'; 
 performance.SR_I3(m,r,:)=(mean(R.I3)*tmp1_p(2,:)'./diag(cov(simret(INDC,:)')))';
 end
 
 
 end;
 
 toc;

% clear the data we don't need
clear qs r ku A N MV sortedValue_X X_Ranked Rep tmp corr lambdas
        
%% Calculate empirical rejection rate for the t-test of the significance of the coefficients
rej_alpha.true=mean(abs(tstat.true(:,:,1))> norminv(0.975),1)';
rej_beta.true=mean(tstat.true(:,:,2)>norminv(0.95),1)';

rej_alpha.ew=mean(abs(tstat.ew(:,:,1))> norminv(0.975),1)';
rej_beta.ew=mean(tstat.ew(:,:,2)>norminv(0.95),1)';

rej_alpha.re01=mean(abs(tstat.re01(:,:,1))> norminv(0.975),1)';
rej_beta.re01=mean(tstat.re01(:,:,2)>norminv(0.95),1)';

rej_alpha.re025=mean(abs(tstat.re025(:,:,1))> norminv(0.975),1)';
rej_beta.re025=mean(tstat.re025(:,:,2)>norminv(0.95),1)'; 

rej_alpha.re05=mean(abs(tstat.re05(:,:,1))> norminv(0.975),1)';
rej_beta.re05=mean(tstat.re05(:,:,2)>norminv(0.95),1)';

% consider the cases of the fixed measurement error with market proxies based on
% 4, 7 and 20 assets
rej_alpha.I2_4=mean(abs(tstat.I2(4,:,:,1))> norminv(0.975),2);
rej_beta.I2_4=mean(tstat.I2(4,:,:,2)>norminv(0.95),2); 
rej_alpha.I2_7=mean(abs(tstat.I2(7,:,:,1))> norminv(0.975),2);
rej_beta.I2_7=mean(tstat.I2(7,:,:,2)>norminv(0.95),2); 
rej_alpha.I2_20=mean(abs(tstat.I2(20,:,:,1))> norminv(0.975),2);
rej_beta.I2_20=mean(tstat.I2(20,:,:,2)>norminv(0.95),2); 

rej_alpha.I3_4=mean(abs(tstat.I3(4,:,:,1))> norminv(0.975),2);
rej_beta.I3_4=mean(tstat.I3(4,:,:,2)>norminv(0.95),2); 
rej_alpha.I3_7=mean(abs(tstat.I3(7,:,:,1))> norminv(0.975),2);
rej_beta.I3_7=mean(tstat.I3(7,:,:,2)>norminv(0.95),2); 
rej_alpha.I3_20=mean(abs(tstat.I3(20,:,:,1))> norminv(0.975),2);
rej_beta.I3_20=mean(tstat.I3(20,:,:,2)>norminv(0.95),2); 

%% OUTPUT TABLES

disp('############################################################################################')
disp('##                                             Output for random error                                                                            ##')
disp('############################################################################################')
OUTPUT=[ 1 mean(RR.re05)  mean(RR.re025) mean(RR.re01) ; % estimated RR 
 1  mean(cor.re05) mean(cor.re025) mean(cor.re01);
  1 mean(lambda.re05) mean(lambda.re025) mean(lambda.re01) ;
  mean(mean(coeff.true(:,:,1),1)) mean(mean(coeff.re05(:,:,1),1)) mean(mean(coeff.re025(:,:,1),1)) mean(mean(coeff.re01(:,:,1),1)) ;
  mean(sqrt(mean(MSE.true(:,:,1),1))) mean(sqrt(mean(MSE.re05(:,:,1),1))) mean(sqrt(mean(MSE.re025(:,:,1),1))) mean(sqrt(mean(MSE.re01(:,:,1),1)));
  mean(rej_alpha.true) mean(rej_alpha.re05) mean(rej_alpha.re025) mean(rej_alpha.re01);
  mean(wald_F_stat.true(:,2)) mean(wald_F_stat.re05(:,2)) mean(wald_F_stat.re025(:,2)) mean(wald_F_stat.re01(:,2));
  mean(wald_F_stat.true(:,1)) mean(wald_F_stat.re05(:,1)) mean(wald_F_stat.re025(:,1)) mean(wald_F_stat.re01(:,1));
  mean(mean(coeff.true(:,:,2),1)) mean(mean(coeff.re05(:,:,2),1)) mean(mean(coeff.re025(:,:,2),1)) mean(mean(coeff.re01(:,:,2),1));
  mean(sqrt(mean(MSE.true(:,:,2),1))) mean(sqrt(mean(MSE.re05(:,:,2),1))) mean(sqrt(mean(MSE.re025(:,:,2),1))) mean(sqrt(mean(MSE.re01(:,:,2),1)));
  sum(mean(coeff.true(:,:,2),1)>1) sum(mean(coeff.re05(:,:,2),1)>1) sum(mean(coeff.re025(:,:,2),1)>1) sum(mean(coeff.re01(:,:,2),1)>1);
  mean(mean(performance.SR_true)) mean(mean(performance.SR_re05)) mean(mean(performance.SR_re025)) mean(mean(performance.SR_re01));
  mean(mean(performance.TR_true)) mean(mean(performance.TR_re05)) mean(mean(performance.TR_re025)) mean(mean(performance.TR_re01))]

toc;

disp('############################################################################################')
disp('##                                                       Output for fixed error                                                                     ##')
disp('############################################################################################')

OUTPUT2=[mean(RR.I2(4,:)) mean(RR.I2(7,:)) mean(RR.I2(20,:)) mean(RR.I3(4,:)) mean(RR.I3(7,:)) mean(RR.I3(20,:)) mean(RR.ew);
mean(cor.I2(4,:)) mean(cor.I2(7,:)) mean(cor.I2(20,:)) mean(cor.I3(4,:)) mean(cor.I3(7,:)) mean(cor.I3(20,:)) mean(cor.ew);
 mean(lambda.I2(4,:)) mean(lambda.I2(7,:)) mean(lambda.I2(20,:)) mean(lambda.I3(4,:)) mean(lambda.I3(7,:)) mean(lambda.I3(20,:)) mean(lambda.ew);
 
 mean(mean(coeff.I2(4,:,:,1),2)) mean(mean(coeff.I2(7,:,:,1),2)) mean(mean(coeff.I2(20,:,:,1),2)) mean(mean(coeff.I3(4,:,:,1),2)) mean(mean(coeff.I3(7,:,:,1),2)) mean(mean(coeff.I3(20,:,:,1),2)) mean(mean(coeff.ew(:,:,1),1));
 mean(sqrt(mean(MSE.I2(4,:,:,1),2))) mean(sqrt(mean(MSE.I2(7,:,:,1),2))) mean(sqrt(mean(MSE.I2(20,:,:,1),2))) mean(sqrt(mean(MSE.I3(4,:,:,1),2))) mean(sqrt(mean(MSE.I3(7,:,:,1),2))) mean(sqrt(mean(MSE.I3(20,:,:,1),2))) mean(sqrt(mean(MSE.ew(:,:,1),1)));
 mean(rej_alpha.I2_4) mean(rej_alpha.I2_7) mean(rej_alpha.I2_20) mean(rej_alpha.I3_4) mean(rej_alpha.I3_7) mean(rej_alpha.I3_20) mean(rej_alpha.ew);
 mean(wald_F_stat.I2(4,:,2)) mean(wald_F_stat.I2(7,:,2)) mean(wald_F_stat.I2(20,:,2)) mean(wald_F_stat.I3(4,:,2)) mean(wald_F_stat.I3(7,:,2)) mean(wald_F_stat.I3(20,:,2)) mean(wald_F_stat.ew(:,2));
 mean(wald_F_stat.I2(4,:,1)) mean(wald_F_stat.I2(7,:,1))  mean(wald_F_stat.I2(20,:,1)) mean(wald_F_stat.I3(4,:,1)) mean(wald_F_stat.I3(7,:,1)) mean(wald_F_stat.I3(20,:,1)) mean(wald_F_stat.ew(:,1));

 mean(mean(coeff.I2(4,:,:,2),2)) mean(mean(coeff.I2(7,:,:,2),2)) mean(mean(coeff.I2(20,:,:,2),2)) mean(mean(coeff.I3(4,:,:,2),2)) mean(mean(coeff.I2(7,:,:,2),2)) mean(mean(coeff.I3(20,:,:,2),2)) mean(mean(coeff.ew(:,:,2),1));
 mean(sqrt(mean(MSE.I2(4,:,:,2),2))) mean(sqrt(mean(MSE.I2(7,:,:,2),2))) mean(sqrt(mean(MSE.I2(20,:,:,2),2))) mean(sqrt(mean(MSE.I3(4,:,:,2),2))) mean(sqrt(mean(MSE.I3(7,:,:,2),2))) mean(sqrt(mean(MSE.I3(20,:,:,2),2))) mean(sqrt(mean(MSE.ew(:,:,2),1)));
 sum(mean(coeff.I2(4,:,:,2),2)>1) sum(mean(coeff.I2(7,:,:,2),2)>1) sum(mean(coeff.I2(20,:,:,2),2)>1) sum(mean(coeff.I3(4,:,:,2),2)>1) sum(mean(coeff.I3(7,:,:,2),2)>1) sum(mean(coeff.I3(20,:,:,2),2)>1) sum(mean(coeff.ew(:,:,2),1)>1);

 mean(mean(performance.SR_I2(4,:,:))) mean(mean(performance.SR_I2(7,:,:))) mean(mean(performance.SR_I2(20,:,:))) mean(mean(performance.SR_I3(4,:,:))) mean(mean(performance.SR_I3(7,:,:))) mean(mean(performance.SR_I3(20,:,:))) mean(mean(performance.SR_ew));
 mean(mean(performance.TR_I2(4,:,:))) mean(mean(performance.TR_I2(7,:,:))) mean(mean(performance.TR_I2(20,:,:))) mean(mean(performance.TR_I3(4,:,:))) mean(mean(performance.TR_I3(7,:,:))) mean(mean(performance.TR_I3(20,:,:))) mean(mean(performance.TR_ew))]

%% OUTPUT FIGURES
%% Figure 1. Correlation between estimated alphas and betas for equally
% weighted index

% Correlation in the case of equally weighted index with 205 stocks
alpha_est=mean(coeff.ew(:,:,1),1);
beta_est=mean(coeff.ew(:,:,2),1);
matrix=corrcoef(alpha_est',beta_est'); % correlation matrix 

% Correlation in the case of equally weighted index based on the subset of
% stocks
correlation=zeros(1,205); % space holder for correlation coefficient
correlation(1,205)=matrix(2,1);
for g=1:204
      % space holder
alpha_est=zeros(1,15);
beta_est=zeros(1,15);
for k=1:15
alpha_est(1,k)=mean(coeff.I3(g,:,k,1),2);
beta_est(1,k)=mean(coeff.I3(g,:,k,2),2);
end
ccc=corrcoef(alpha_est',beta_est');
correlation(1,g)=ccc(2,1);    % store correlation coefficient
end 
figure(1)
plot(correlation,'LineWidth',5)
xlim([0 205])
title('Correlation coefficient between estimated alphas and betas')
xlabel('Number of the largest stocks included in the index','FontSize',20)
ylabel('Correlation coefficient','FontSize',20)
hold on
zeroline = plot([0,205],[0,0],'r-','LineWidth',5); % maye first argument xlim
hold off

%% Figure 2. Empirical rejection rate of the absence of abnormal returns in the CAPM with the 
% index based on the subset of largest stocks
errT.I2=mean(mean(abs(tstat.I2(:,:,:,1))> norminv(0.975),2),3);   % empirical rejection rate for t-test
errT.I3=mean(mean(abs(tstat.I3(:,:,:,1))> norminv(0.975),2),3);   % empirical rejection rate for t-test
errT.I2(205,1)=mean(mean(abs(tstat.true(:,:,1))> norminv(0.975),1),2); 
errT.I3(205,1)=mean(mean(abs(tstat.ew(:,:,1))> norminv(0.975),1),2); 

errF.I2(1:204,1)=mean(wald_F_stat.I2(:,:,2),2); % empirical rejection rate for F-test in case of the index based on the subset of stocks
errF.I3(1:204,1)=mean(wald_F_stat.I3(:,:,2),2); % empirical rejection rate for F-test in case of the equally-weighted index based on the subset of stocks
errF.I2(205,1)=mean(wald_F_stat.true(:,2));      % empirical rejection rate for F-test in case of the true market index
errF.I3(205,1)=mean(wald_F_stat.ew(:,2));        % empirical rejection rate for F-test in case of the equally-weighted index based on 205 stocks

%space holders
power_true=zeros(1,101); power_ew=zeros(1,101); power_I2_4=zeros(1,101); power_I2_7=zeros(1,101);
power_I2_20=zeros(1,101); power_I3_4=zeros(1,101); power_I3_7=zeros(1,101); power_I3_20=zeros(1,101);

signlevels=[0.0:0.01:1];

%the ROC curve of the F-test
for i=1:101
power_true(1,i)=mean(mean(abs(tstat.true(:,:,1))> norminv(1-signlevels(1,i)/2),1)'); 
power_ew(1,i)=mean(wald_F_stat.true(:,1)>finv(1-signlevels(1,i),CS,T-CS-1));
power_I2_4(1,i)=mean(wald_F_stat.I2(4,:,3)>finv(1-signlevels(1,i),CS,T-CS-1),2);
power_I2_7(1,i)=mean(wald_F_stat.I2(7,:,3)>finv(1-signlevels(1,i),CS,T-CS-1),2);
power_I2_20(1,i)=mean(wald_F_stat.I2(20,:,3)>finv(1-signlevels(1,i),CS,T-CS-1),2);

power_I3_4(1,i)=mean(wald_F_stat.I3(4,:,3)>finv(1-signlevels(1,i),CS,T-CS-1),2);
power_I3_7(1,i)=mean(wald_F_stat.I3(7,:,3)>finv(1-signlevels(1,i),CS,T-CS-1),2);
power_I3_20(1,i)=mean(wald_F_stat.I3(20,:,3)>finv(1-signlevels(1,i),CS,T-CS-1),2);
end
% for the ROC curve of the t-test (comment if not needed)
for i=1:101
power_true(1,i)=mean(mean(abs(tstat.true(:,:,1))> norminv(1-signlevels(1,i)/2),1)');  
power_ew(1,i)=mean(mean(abs(tstat.ew(:,:,1))> norminv(1-signlevels(1,i)/2),1)');
power_I2_4(1,i)=mean(mean(abs(tstat.I2(4,:,:,1))> norminv(1-signlevels(1,i)/2),2));
power_I2_7(1,i)=mean(mean(abs(tstat.I2(7,:,:,1))> norminv(1-signlevels(1,i)/2),2));
power_I2_20(1,i)=mean(mean(abs(tstat.I2(20,:,:,1))> norminv(1-signlevels(1,i)/2),2));

power_I3_4(1,i)=mean(mean(abs(tstat.I3(4,:,:,1))> norminv(1-signlevels(1,i)/2),2));
power_I3_7(1,i)=mean(mean(abs(tstat.I3(7,:,:,1))> norminv(1-signlevels(1,i)/2),2));
power_I3_20(1,i)=mean(mean(abs(tstat.I3(20,:,:,1))> norminv(1-signlevels(1,i)/2),2));
end


% Grid for the graph
numberstocks=[1:1:205];                       % number of largest stocks chosen for construction of the proxy
numberstocks=numberstocks*100/205;   % percentage of largest stocks chosen for construction of the proxy
gr=[1:4:205];


figure(2)
subplot(2,2,1)
p1 = plot(numberstocks(gr),errT.I2(gr),'b+');
xlabel('Percentage of the largest stocks included','FontSize',14)
ylabel('ERR(.05)','FontSize',14)
hold on
p2 = plot(xlim,[0.05,0.05],'r-');
set(p1,{'LineWidth'},{3})
hold off
title('ERR(.05) for the index with normalized weights','FontSize',14)

subplot(2,2,2)
p1 = plot(numberstocks(gr),errT.I3(gr),'b+');
xlabel('Percentage of the largest stocks included','FontSize',14)
ylabel('ERR(.05)','FontSize',14)
hold on
p2 = plot(xlim,[0.05,0.05],'r-');
set(p1,{'LineWidth'},{3})
hold off
title('ERR(.05) for the index with equal weights','FontSize',14)

subplot(2,2,3)
plot(signlevels,power_true)
xlabel('Size','FontSize',14,'Interpreter','Latex')
ylabel('Power','FontSize',14,'Interpreter','Latex')
hold on
f1=plot(signlevels,power_I2_4,'-.')
f2=plot(signlevels,power_I2_7,':k')
f3=plot(signlevels,power_I2_20,'--')
hold off
title('ROC curve for the index with normalized weights','FontSize',14)
lgt = legend([f1, f2, f3], 'M=4', 'M=7', 'M=20');
set(f1,{'LineWidth'},{3})
set(f2,{'LineWidth'},{3})
set(f3,{'LineWidth'},{3})
lgt.Location='southeast';
lgt.FontSize=14;

subplot(2,2,4)
plot(signlevels,power_true)
xlabel('Size','FontSize',14,'Interpreter','Latex')
ylabel('Power','FontSize',14,'Interpreter','Latex')
hold on
l1=plot(signlevels,power_I3_4,'-.')
l2=plot(signlevels,power_I3_7,':k')
l3=plot(signlevels,power_I3_20,'--')
hold off
title('ROC curve for the index with equal weights','FontSize',14)
lgt = legend([l1, l2, l3], 'M=4', 'M=7', 'M=20');
set(l1,{'LineWidth'},{3})
set(l2,{'LineWidth'},{3})
set(l3,{'LineWidth'},{3})
lgt.Location='southeast';
lgt.FontSize=14;


%% Figure 3. Security market line and CAPM estimates with and without measurement error
slope_estimates=[mu(INDC) mean(coeff.true(:,:,2),1)' mean(coeff.ew(:,:,2),1)']

figure(3)
p1 = plot(slope_estimates(:,2),slope_estimates(:,1),'b+');
xlabel('\beta','FontSize',14)
ylabel('E[r_j] - r_f','FontSize',14)
hold on
p2 = plot([0,max(slope_estimates(:,2))],[0,mean(mu_m_star)*max(slope_estimates(:,2))],'r-');
p4 = plot([0,max(slope_estimates(:,2))],[0,mean(SML.ew)*max(slope_estimates(:,2))],'b-');
p3 = plot(slope_estimates(:,3),slope_estimates(:,1),'go');
set(p2,{'LineWidth'},{2})
set(p1,{'LineWidth'},{2})
set(p3,{'LineWidth'},{2})
lgt = legend([p1, p3, p2, p4], 'true \beta', '\beta^{EW}', 'true SML','SML^{EW}' );
lgt.Location='southeast';
lgt.FontSize=14;
hold off

%% Plots of the Treynor ratios
figure(4)
p1=boxplot(performance.TR_true(:,:))
set(p1,{'LineWidth'},{2})
title('Treynor ratio estimates in case of known r^*_{m,t}')
ylim([0 20])
ylabel('TR')
h=findobj(gca,'tag','Outliers');
delete(h)
hold on
p2=plot(xlim,[mu_m_star(1,1) mu_m_star(1,1)])
set(p2,{'LineWidth'},{2})
hold off

figure(5)
p1=boxplot(performance.TR_ew(:,:))
set(p1,{'LineWidth'},{2})
title('Treynor ratio estimates in case of equally-weighted index')
ylim([0 20])
ylabel('TR')
h=findobj(gca,'tag','Outliers');
delete(h)
hold on
p2=plot(xlim,[mu_m_star(1,1) mu_m_star(1,1)])
set(p2,{'LineWidth'},{2})
hold off

figure(6)
p1=boxplot(performance.TR_re05(:,:))
set(p1,{'LineWidth'},{2})
title('Treynor ratio estimates in case of random measurement error with \sigma_{\nu}=0.05')
ylim([-20 80])
ylabel('TR')
h=findobj(gca,'tag','Outliers');
delete(h)
hold on
p2=plot(xlim,[mu_m_star(1,1) mu_m_star(1,1)])
set(p2,{'LineWidth'},{2})
hold off
%% Figure 5. Sharpe ratios
SRrefer=(mean(mu(INDC,:),2)./diag(cov(simret(INDC,:)')))'; % true Sharpe ratio
SRrefer=repmat(SRrefer,10,1);

figure(7)
p1=boxplot(performance.SR_re05(:,:))
set(p1,{'LineWidth'},{2})
title('Sharpe ratio estimates in case of random measurement error with \sigma_{\nu}=0.05')
ylabel('SR')
h=findobj(gca,'tag','Outliers');
delete(h)
hold on
p2=boxplot(SRrefer(:,:), 'Color', 'g')
set(p2,{'LineWidth'},{2})
ylim([-0.1 0.3])
hold off

figure(8)
p1=boxplot(performance.SR_re025(:,:))
set(p1,{'LineWidth'},{2})
title('Sharpe ratio estimates in case of random measurement error with \sigma_{\nu}=0.025')
ylabel('SR')
h=findobj(gca,'tag','Outliers');
delete(h)
hold on
p2=boxplot(SRrefer(:,:), 'Color', 'g')
set(p2,{'LineWidth'},{2})
ylim([-0.1 0.3])
hold off

figure(9)
p1=boxplot(performance.SR_re01(:,:))
set(p1,{'LineWidth'},{2})
title('Sharpe ratio estimates in case of random measurement error with \sigma_{\nu}=0.01')
ylabel('SR')
h=findobj(gca,'tag','Outliers');
delete(h)
hold on
p2=boxplot(SRrefer(:,:), 'Color', 'g')
set(p2,{'LineWidth'},{2})
ylim([0 0.3])
hold off
%% Check of the autocorrelations of the returns in the simulation study
subplot(2,1,1)
autocorr(Rmstar)
subplot(2,1,2)
parcorr(Rmstar)

subplot(2,1,1)
autocorr(R.I3)
subplot(2,1,2)
parcorr(R.I3)

subplot(2,1,1)
autocorr(simret(1,:))
subplot(2,1,2)
parcorr(simret(1,:))
