% two_pass_regression_GMM.m

% you need to supply excess_returns, the test assets and factors
% is called from the two_pass_GMM_all.m

analysis=2;             % use 2 for regression-based estimation and test
                        % -ing of asset pricing models
                        
factor_type=2;          % 2: cross-sectional regression: 
                        % factors are not excess returns         
                        
GLS=0;                  % OLS

%% parameters, variance-covariance matrix and standard errors

% get parameters
param=two_pass(factors, excess_returns);

% compute asymptotic variance covariance matrix for average pricing errors
% and for parameter estimates: general GMM formulae!
[g_T_var,par_var]=var_GMM( param, analysis, excess_returns, factors, factor_type, c_fac, mim_proj_var, mim_base_assets );

% standard errors of the estimates
param_se_OLS=(sqrt(diag(par_var)))/(sqrt(size(excess_returns,1)));
param_est_OLS=param;

display('*************************************************');
display('*      Output of analysis                       *');
display('*************************************************');
display('*      Estimates and t-stats:                   *');
display('*      first N alpha, N*2 betas, K lambdas      *');
disp([param_est_OLS param_est_OLS./param_se_OLS]);
display('*************************************************');

%% beta*lambda analysis

% compute beta*lambda for each factor (in rows) assets (columns)
beta_lambda_r=get_beta_lambda( param, analysis, excess_returns, factors, c_fac );

% summing up the elements of each column gives the model-implied expected
% return per test asset (you can plot that on the 45 degree line plot)

if size(beta_lambda_r',2)==1
model_implied_expected_return=beta_lambda_r';           %% Achtung hier abgeändert!!! von sum zu so, damit die dimensionen passen
else
model_implied_expected_return=(sum(beta_lambda_r))';
end

% compare realized and predicted expected return
display('*  observed excess return vs. model implied     *');
disp([mean(excess_returns)' model_implied_expected_return]);
display('*************************************************');

%% 45 degree plots

% get a 45 degree line: auxiliary variables
max_deg=max([max(model_implied_expected_return),max(mean(excess_returns))])+0.005;


h1=figure(1);
scatter(model_implied_expected_return,mean(excess_returns)','x','black')
title('predicted vs. model-implied pricing errors')
xlabel('predicted excess returns')
ylabel('observed excess returns')
xlim([0,max_deg])
ylim([0,max_deg])
refline(1,0);

% if freq==2
% xtick=[0,0.05,0.01,0];    
% ytick=[0,0.05,0.01,0];
% elseif freq==1
% xtick=[0,0.2,0.05,0];    
% ytick=[0,0.2,0.05,0];    
% end
saveas(h1,'fortyfive.eps');

%% J-test: statistic for validity of the asset pricing model

[ J_last_N, p_value_J_last_N, g_T_last_N, t_stat_last_N, corr_g_T, Avar_g_T_last_N ] = Jstat_GRS( param, analysis, excess_returns, factors, factor_type, c_fac, mim_proj_var, mim_base_assets );

display('*************************************************');
display('* J-statistic from GMM '); disp(J_last_N);
display('* p-value (%) '); disp(p_value_J_last_N);
display('*************************************************');
display('* Avg. pricing errors (g_T) and t-stats         *'); 
disp([g_T_last_N t_stat_last_N]);
display('*************************************************');

%% compute variance covariance matric of vectorized beta*lambda and the standard errors of the vectorized beta*lambda

[ beta_lambda_vec_vc, se_beta_lambda ] = beta_lambda_vc( param, par_var, analysis, excess_returns, factors, c_fac);

display('*************************************************');
display('* each row is a factor for which we have        *');
display('* beta * lambda for the stocks                  *'); 
display('* Beta * lambda                                 *'); 
disp(beta_lambda_r);
display('*************************************************');
display('* standard error beta * lambda                  *');
disp(se_beta_lambda);
display('*************************************************');
display('* t-statistics beta * lambda                    *');
disp(beta_lambda_r./se_beta_lambda);
display('*************************************************');

if testassets==3                % 17 industry portfolios
    if model==3 || model==4     % for CDRM variants
        industry_beta_lambda=beta_lambda_r(2,size(beta_lambda_r,2)-17+1:size(beta_lambda_r,2));
        industry_beta_lambda_se=se_beta_lambda(2,size(beta_lambda_r,2)-17+1:size(beta_lambda_r,2));
        industry_excess_returns=excess_returns(:,size(beta_lambda_r,2)-17+1:size(beta_lambda_r,2));
        industry_pricing_errors=g_T_last_N(size(g_T_last_N,1)-17+1:size(g_T_last_N,1),:);
        industry_pricing_errors_t=t_stat_last_N(size(g_T_last_N,1)-17+1:size(g_T_last_N,1),:);
 
        display('* mean ex. ret. || Beta_N x Lambda_N x 100 17 Industries || t-stats beta * lambda || mean pricing error || t-stat ');               
        disp([mean(industry_excess_returns)'*100,industry_beta_lambda'*100,(industry_beta_lambda./industry_beta_lambda_se)',industry_pricing_errors*100,industry_pricing_errors_t]);
    end
end

% now format beta*lambda results into matrices used in paper

display('* Beta * lambda result by factor                *');
display('*************************************************');

if testassets==1 || testassets==2
    factor=1;
    % selects the first factor (first row) - market return for all three
    % models
    
    % reshape beta*lambda to write in paper format
    display('* Beta*lambda result for first factor (market)  *');
    display('*************************************************');
    
    beta_lambda_mat_factor_1=reshape(beta_lambda_r(factor,:),[5,5]);
    se_beta_lambda_mat_factor_1=reshape(se_beta_lambda(factor,:),[5,5]);
    t_beta_lambda_mat_factor_1=reshape(beta_lambda_r(factor,:)./se_beta_lambda(factor,:),[5,5]);
    
    disp([strcat('* Beta*lambda for selected factor ( ',num2str(factor)),strcat(' ) ')]);
    display('* horizonal: Book-to-market - vertical: size    *');
    disp(beta_lambda_mat_factor_1');
    display('*************************************************');  
    disp([strcat('* standard error matrix beta*lambda for selected factor ( ',num2str(factor)),strcat(' ) ')]);
    disp(se_beta_lambda_mat_factor_1');
    display('*************************************************');  
    disp([strcat('* t-stats matrix beta*lambda for selected factor ( ',num2str(factor)),strcat(' ) ')]);
    disp(t_beta_lambda_mat_factor_1');
    display('*************************************************');
    disp([strcat('* p-value (%) matrix beta*lambda for selected factor ( ',num2str(factor)),strcat(' ) ')]);
    disp(((1-normcdf(abs(t_beta_lambda_mat_factor_1)))*100*2)');
else
    display('Not implemented for test assets other than B/M portfolios');
end


if model==2 || model==3 || model==4     % models with a least 2 factors
   if testassets==1 || testassets==2
       % select factor number (row) 
       factor=2;    
       % select the 2 PAG factor, 1 market
       
       % reshape beta*lambda to write in paper format
       display('* Beta*lambda result for second factor (PAG for CDR model) *');
       display('************************************************************');  
       
       beta_lambda_mat_factor_2=reshape(beta_lambda_r(factor,:),[5,5]);
       se_beta_lambda_mat_factor_2=reshape(se_beta_lambda(factor,:),[5,5]);
       t_beta_lambda_mat_factor_2=reshape(beta_lambda_r(factor,:)./se_beta_lambda(factor,:),[5,5]);
       
       disp([strcat('* Beta*lambda for selected factor ( ',num2str(factor)),strcat(' ) ')]);
       display('* horizonal: Book-to-market - vertical: size    *');
       disp(beta_lambda_mat_factor_2');
       display('*************************************************');  
       disp([strcat('* standard error matrix beta*lambda for selected factor ( ',num2str(factor)),strcat(' ) ')]);
       disp(se_beta_lambda_mat_factor_2');
       display('*************************************************');  
       disp([strcat('* t-stats matrix beta*lambda for selected factor ( ',num2str(factor)),strcat(' ) ')]);
       disp(t_beta_lambda_mat_factor_2');
       display('*************************************************');
       disp([strcat('* p-value (%) matrix beta*lambda for selected factor ( ',num2str(factor)),strcat(' ) ')]);
       disp(((1-normcdf(abs(t_beta_lambda_mat_factor_2)))*100*2)');
       
       if freq==2 && model==3;          % for CDRM: annualized compensations
          display('annualized compensation (((beta_N x lambda_N +1)^4)-1)*100'); 
          disp((((beta_lambda_mat_factor_2+1)^4)-1)*100);
          display('annualized beta_N x lambda_N x 4 x 100');  
          disp(beta_lambda_mat_factor_2*4*100);
       end
       
       % Difference in expected return analysis for CDR model
       if model==3 || model==4  
           % This is only done for CDR model (possibly with mimicking PF)
           
           % Test the significance of beta*lambda difference of two stocks.
           
           % Of the last three arguments, the first is the row that
           % identifies the factor (differences can only be tested holding
           % the factor constant, the remaining two are the assets. The big
           % growth portfolio is number 21, the small value portfolio is 5.
           % We pick the second row, that is the PAG factor
           
           % NOTE that fit the CDR model, PAG must be the second factor !!
           
           display('*********************************************************');  
           display('* Differences in expected return analysis for CDR model *'); 
           
           asset_1=5;           % first asset to be compared (small value)
           asset_2=21;          % second asset (big growth)
           
           [ beta_lambda_diff_r, beta_lambda_se, beta_lambda_pval ] = get_beta_lambda_diff( beta_lambda_r, beta_lambda_vec_vc, factor, asset_1, asset_2, excess_returns );
           display('Difference beta*lambda for factor'); disp(factor);
           disp([strcat('Assets  =  ',num2str(asset_1)),strcat(' - ',num2str(asset_2))]);
           disp('diff.  -  s.e.  -  t-stat  -  p_value in % ');
           disp([beta_lambda_diff_r,beta_lambda_se,beta_lambda_diff_r./beta_lambda_se,beta_lambda_pval*100]);
           display('Analysis should be for asset 5, small value, against 21, big growth');
       end
       
   else
       display('Not implemented for test assets other than 1 or 2');
   end
end


if model==2                             % models with three factors: 
                                        % Fama-French model
   if testassets==1 || testassets==2
       % select factor number (row) 
       factor=3;    
       % select the 2 PAG factor, 1 market
       
       % reshape beta*lambda to write in paper format
       display('* Beta*lambda result for third factor (only FF model) *');
       display('*******************************************************');    
       
       beta_lambda_mat_factor_3=reshape(beta_lambda_r(factor,:),[5,5]);
       se_beta_lambda_mat_factor_3=reshape(se_beta_lambda(factor,:),[5,5]);
       t_beta_lambda_mat_factor_3=reshape(beta_lambda_r(factor,:)./se_beta_lambda(factor,:),[5,5]);
       
       disp([strcat('* Beta*lambda for selected factor ( ',num2str(factor)),strcat(' ) ')]);
       display('* horizonal: Book-to-market - vertical: size    *');
       disp(beta_lambda_mat_factor_3');
       display('*************************************************');  
       disp([strcat('* standard error matrix beta*lambda for selected factor ( ',num2str(factor)),strcat(' ) ')]);
       disp(se_beta_lambda_mat_factor_3');
       display('*************************************************');  
       disp([strcat('* t-stats matrix beta*lambda for selected factor ( ',num2str(factor)),strcat(' ) ')]);
       disp(t_beta_lambda_mat_factor_3');
       display('*************************************************');
       disp([strcat('* p-value (%) matrix beta*lambda for selected factor ( ',num2str(factor)),strcat(' ) ')]);
       disp(((1-normcdf(abs(t_beta_lambda_mat_factor_3)))*100*2)');
       display('*************************************************');
   else
       display('Not implemented for test assets other than 1 or 2');
   end
    
end
