clc
clear

c_fac=[];

factor_type=2;          % 2: cross-sectional regression: 
                        % factors are not excess returns   

% TwoPassGMM_JFQA.m

% Frontend for estimation of CDRM and other models by two-pass regression
% using two-pass-regression/ GMM a'la Cochrane produces results in paper
% "Creative Destruction and Asset Prices"

% Results are printed on screen

%% Initialization of globals

% 1. set here the analysis you want
model=1;        % relevant choices are:
                % 1=CAPM  
                % 2=Fama French model
                % 3=CDRM (referred to as PAG in the code)
                % 4=MP-CDRM with excess returns of mimicking portfolios
                % supplied in raw data (used for 45 degree plots and
                % covariances of excess returns with patenting growth
                % proper GMM inference is with model 6)
                
                % 5=GMM NOT IMPLEMENTED
                % inference for mimicking portfolio estimation, Wald
                % test for joint significance of mimicking portfolio
                % weights. Also performs regressions on SMB and HML to test
                % the significance of intercepts taking into account that
                % the portfolio weights are estimated
                
                % 6=MP-CDRM three-pass regression using GMM results,
                % including Lewellen/Wald test
                % 8=orthogonalized extended CDRM (orthogonalized HML and
                % SMB w.r.t. patenting grwoth)
                % 10=CDH portfolio model
                
% NOTE: model=4 uses a pre-estimated mimicking portfolio while model=6 
% estimates it and provides proper inference. The estimates are numeri-
% cally identical. Use model==4 to produce the goodness-of-fit plots

testassets=1;   % 1=25sb, value-weighted (VW)
                % 2=25sb, equally-weighted (EW)
                % 3=sb25+ind17ew, VW
                % 4=HML and SMB
                % 5=10 size single sorted (VW)
                % 6=10 B/M single sorted (VW)
                
freq=1;         % 1= annually
                % 2= quarterly
                
                
% Main test reports results on
% model - test asset - frequency

% 1-1-1 (CAPM, value-weighted, annual)
% 2-1-1 (Fama-French, value-weighted, annual)
% 3-1-1 (CDRM, value-weighted, annual)
% 3-1-2 (CDRM, value-weighted, quarterly)
% 4-1-1 (MP-CDRM, value-weighted, annually) - just plots
% 5- -1 (HML and SMB as test assets, with MP-CDRM
% 6-1-1 (MP-CDRM, value-weighted, annually)
% 6-2-1 (MP-CDRM, value-weighted, quarterly)
% 6-4-1 (MP-CDRM with SMB and HML as test assets
% online appendix also has results for equally weighted portfolios

if freq==2 && testassets==3
    display('*************************************************');
    display('*  not implemented                              *');
    display('*************************************************');
    return;
end
if model==4 && freq==2
    display('*************************************************');
    display('*  not implemented                              *');
    display('*************************************************');
    return;  
end

display('*************************************************');
disp(strcat(' Analysis for model =  ',num2str(model)));
disp(strcat(' Frequency          =  ',num2str(freq)));
disp(strcat(' Test assets        =  ',num2str(testassets)));
display('*************************************************');

%% Manage data input

% 2. call function that loads in the data

alreadyin=0;

[ excess_returns, factors, mim_base_assets, mim_proj_var, pag, pagm ] = data_input( freq, alreadyin, testassets, model );

% define number of test assets, number of factors and length of time-series
N=size(excess_returns,2);
K=size(factors,2);
T=size(factors,1);

%% run the two-pass regression/ GMM for CAPM, Fama-French, CDRM or MP-CDRM

if model<=4 || model==8 || model==10
    % not for mimicking pf but the three base models
    
two_pass_regression_GMM;        % call this file: universal result for all 
                                % models 1-4
                                


%% Postprogressing of estimation results for model 1-4

% The two_pass_regression_GMM.m file does not prepare the betas and their
% standard errors in a convenient way, i.e in a way they are reported in
% the paper. This is done now: beta estimates come as a N*K vector (after
% the N alpha estimates, the betas for each asset are stacked and need to
% be reshaped to have them separate for each factor for reporting

% compute t-statistics for the estmiates
t=param_est_OLS./param_se_OLS;

betas=param_est_OLS(N+1:N*(K+1));
betas_se=param_se_OLS(N+1:N*(K+1));
betas_t=t(N+1:N*(K+1));

% The N*K vectors are reshaped as KxN matrices: the first column has the
% betas for the first test asset, and so forth, the first row holds the
% beta values belong to the first factor and so on: this is the lower right
% corner of the a_T matrix

betas_mat=reshape(betas,[K,N])';
betas_se_mat=reshape(betas_se,[K,N])';
betas_t_mat=reshape(betas_t,[K,N])';

% lambdas
lambdas_mat=param_est_OLS(N*(K+1)+1:size(param_est_OLS,1));
lambdas_se_mat=param_se_OLS(N*(K+1)+1:size(param_est_OLS,1));
lambdas_t_mat=t(N*(K+1)+1:size(param_est_OLS,1));

display('*************************************************');
display('* Estimated lambdas                             *');
disp(lambdas_mat);
display('* Lambdas t-stats                               *');
disp(lambdas_t_mat);
display('* Lambdas p-values in %                         *');
disp((1-normcdf(abs(lambdas_t_mat)))*100*2);
display('*************************************************');

% J-test results from the two_pass_regression_GMM.m file, we have the test
% statistic, the d.f. and the p-value

J_GRS_results=[J_last_N,N-K,p_value_J_last_N];

% R-squared along the lines of Maio Santa-Clara 2012
mexret_fit=model_implied_expected_return;
mexret=mean(excess_returns)';
alpha_hat=mexret-mexret_fit;
alpha_hat_demeaned=alpha_hat-mean(alpha_hat);
mexret_demeaned=mexret-mean(mexret);

% R^2
rsq=1-((sum(alpha_hat_demeaned.^2))/(sum(mexret_demeaned.^2)));
% uncentered R^2
rsq_uncentered=1-((sum(alpha_hat.^2))/(sum(mexret.^2)));

MAE=mean(abs(alpha_hat));

display('*************************************************');
display('*  Cross-sectional R_squared (Eqn(15))          *');
disp(rsq*100);
display('*************************************************');
display('*  Uncentered cross-sectional R_squared         *');
disp(rsq_uncentered*100);
display('*************************************************');
display('*  Mean absolute pricing error (MAE) (Eqn(14))  *');
disp(MAE*100);
display('*************************************************');

if model==3 || model==4     
    % 3: CDRM or 4: MP-CDRM
    
    % now reshape the beta estimates, s.e. and t in 5 by 5 matrices as
    % reported in the paper
    
    if testassets==1 || testassets==2
        CDR_annual_FF_beta_mkt=reshape(betas_mat(:,1),[5,5]);
        CDR_annual_FF_beta_PAG=reshape(betas_mat(:,2),[5,5]);
        
        display('*  Invention Betas  *');
        disp(CDR_annual_FF_beta_PAG');
        
        CDR_annual_FF_beta_mkt_se=reshape(betas_se_mat(:,1),[5,5]);
        CDR_annual_FF_beta_PAG_se=reshape(betas_se_mat(:,2),[5,5]);
        
        CDR_annual_FF_beta_mkt_t=reshape(betas_t_mat(:,1),[5,5]);
        CDR_annual_FF_beta_PAG_t=reshape(betas_t_mat(:,2),[5,5]);   
        
        display('*  Invention Betas t-stats  *');
        disp(CDR_annual_FF_beta_PAG_t');
        
        display('*  Invention Betas p_values (in %)  *');
        disp(((1-normcdf(abs(CDR_annual_FF_beta_PAG_t)))*100*2)');
        
        % Beta*lambda directly from two_pass_regression_GMM.m
        
        CDR_annual_FF_betalambda_mkt=beta_lambda_mat_factor_1;
        CDR_annual_FF_betalambda_PAG=beta_lambda_mat_factor_2;
        
        CDR_annual_FF_betalambda_mkt_se=se_beta_lambda_mat_factor_1;
        CDR_annual_FF_betalambda_PAG_se=se_beta_lambda_mat_factor_2;
        
        CDR_annual_FF_betalambda_mkt_t=t_beta_lambda_mat_factor_1;
        CDR_annual_FF_betalambda_PAG_t=t_beta_lambda_mat_factor_2;        
    else
        display('*  not implemented for test assets other than 25 B/M sorted portfolios  *');
    end

else
    display('*  not implemented  *');
    
end


% Time Series R-squared

r_sq_mat  = r_squared( param, excess_returns, factors , analysis, c_fac);
 
 
display('*************************************************');
display('*  Time Series R-squared                        *');
disp(r_sq_mat');
display('*************************************************');



if model==2 || model==3 || model==4 || model==8 || model==9 || model==10
   % ICAPM consistency based on beta-representation
   display('*  Plausibility of the implied risk aversion coefficients  *');   
   display('*  of the ICAPM (see Maio and Santa-Clara, 2012)           *');   
   
   gamma_vec=(repmat(1/T,[size(factors,2) size(factors,2)]).*(factors'*factors)-mean(factors)'*mean(factors))\lambdas_mat;
   
   % compare Eqn 12 in paper
   vc_lambda=par_var(size(par_var,1)-K+1:size(par_var,1),size(par_var,1)-K+1:size(par_var,1));
   
   vc_gamma=(repmat(1/T,[size(factors,2) size(factors,2)]).*(factors'*factors)-mean(factors)'*mean(factors))\(vc_lambda/(T))/((repmat(1/T,[size(factors,2) size(factors,2)]).*(factors'*factors)-mean(factors)'*mean(factors)));
   
   se_gamma=sqrt(diag(vc_gamma));
   
display('*************************************************');
display('*  implied gamma estimates and t-statistics (Eqn(12)) *');
disp([gamma_vec,gamma_vec./se_gamma]);
display('*************************************************');   
disp(strcat('* implied rra estimate  =  ',num2str(gamma_vec(1))));
disp([strcat('* with s.e.  =  ',num2str(se_gamma(1))),strcat(' with t-stat  =  ',num2str(gamma_vec(1)/se_gamma(1)))]);  

   if model==3 || model==4
       disp(strcat('* implied rra estimate  =  ',num2str(gamma_vec(2))));
       disp([strcat('* with s.e.  =  ',num2str(se_gamma(2))),strcat(' with t-stat  =  ',num2str(gamma_vec(2)/se_gamma(2)))]);  
       
       if model==4
          disp('Warning, inference on RRA and gamma_n does not take into account that mimicking weights are estmiated');
       end
       
   end
end

if model==3 || model==4
    % covariances of test assets excess returns with pg
    if model==3
        DataMat=[pag,excess_returns];
    elseif model==4
        DataMat=[pagm,excess_returns];
    end
    
    vcmatrix=repmat(1/T,[size(DataMat,2) size(DataMat,2)]).*(DataMat'*DataMat)-mean(DataMat)'*mean(DataMat);
    
    interesting_covs=vcmatrix(2:size(vcmatrix,1),1);
    
    if testassets==1 || testassets==2
        % only B/M sorted portfolios
        
       covmat_25=reshape(interesting_covs,[5,5]); 
       
       display('*************************************************');
       display('*  Covariances of patent growth (model=3) or    *');   
       display('*  mimicking portfolio excess return (model=4)  *');
       display('*  with excess returns arranged in 5x5 form     *'); 
       disp(covmat_25'*100);
       display('*************************************************');
       
    else
       display('*  not implemented    *');        
    end
    
end


elseif model==6     % mimicking portfolio CDRM with GMM
    % Postprogressing estimation results for model 6
    analysis=7;     % required in generic GMM procedures
    
    c_fac=K+size(mim_proj_var,2);
    
    param = get_parameters( mim_base_assets, mim_proj_var, analysis, excess_returns, factors );
    % estimates by OLS mimicking weights, betas and lambdas separately
    
    betas=param(N+1:N*(c_fac+1));
    betas_mat=reshape(betas,[c_fac,N])';
    
    % report estimates and t-stats as required for the paper
    betas_market=betas_mat(:,1);
    betas_pagm=betas_mat(:,2);
    
    if testassets==1 || testassets==2   %25 B/M sorted test assets
        betas_market_mat=reshape(betas_market,[5,5])';
        betas_pagm_mat=reshape(betas_pagm,[5,5])';
    else
       display(' not implemented '); 
    end
    
    S_est  = S_estimate( param, analysis, excess_returns, factors, factor_type, c_fac, mim_proj_var, mim_base_assets );
    [ ~,Avar_param ] = var_GMM( param, analysis, excess_returns, factors, factor_type, c_fac, mim_proj_var, mim_base_assets );
    
    [ a_T ] = get_a_T( param, excess_returns, analysis, factors, factor_type, c_fac, mim_base_assets);
    [ g_T ] = get_g_T( param, analysis, excess_returns, factors, factor_type, c_fac, mim_proj_var, mim_base_assets );
    
    display('*  check whether a_T * g_T = 0                  *');
    display('*************************************************');    
    display((a_T*g_T)');
    
    display('*************************************************');
    display('*  Parameters and t-stats:                      *');
    display('*  first alphas and betas, then lambdas         *');
    display('* (CS and TS), then mimicking                   *');
    display('*************************************************');
    t_stat=param./sqrt(diag(Avar_param)./size(mim_proj_var,1));
    disp([(1:size(param,1))',param,t_stat]);
    display('*************************************************');
    
    t_stat=param./sqrt(diag(Avar_param)./size(mim_proj_var,1));
    t_stat_betas=t_stat(N+1:N*(c_fac+1));
    t_stat_mat=reshape(t_stat_betas,[c_fac,N])';
    t_stat_market=t_stat_mat(:,1);
    t_stat_pagm=t_stat_mat(:,2);
    
    if testassets==1 || testassets==2
        t_stat_market_mat=reshape(t_stat_market,[5,5])';
        t_stat_pagm_mat=reshape(t_stat_pagm,[5,5])';
    else
        disp('not implemented for test assets other than 25 B/M');
    end
    
    rsq_vec=zeros(N,1);
    i=1;
    
    while i <= N
        % computes the time series R-squared
        mimportfolio=mim_base_assets*param(N*(c_fac+1)+2*c_fac+size(mim_proj_var,2)+1:N*(c_fac+1)+2*c_fac+size(mim_proj_var,2)+size(mim_base_assets,2));
        
        % construct Rē
        y=excess_returns(:,i);
        X=[ones(T,1),mimportfolio,factors];
        
        b=(X'*X)\(X'*y);        % 1. run regression of X on y
        u=y-X*b;                % 2. get residual of regression
        u_tilde=u-mean(u);
        y_tilde=y-mean(y);
        
        SS_u=u_tilde'*u_tilde;
        SS_y=y_tilde'*y_tilde;
        
        rsq=1-(SS_u/SS_y);      % 3. Rē
        
        
        rsq_vec(i)=rsq;
        i=i+1;
    end
    
    % reshape time series R-squared as required for table in the paper
    if testassets==1 || testassets==2
        rsq_mat=reshape(rsq_vec,[5,5])';
        display('*************************************************');
        display('*  Time Series R-squared                        *');
        disp(rsq_mat);
        display('*************************************************');
    else
        disp('not implemented for test assets other than 25 B/M sorted portfolios');
    end
    
    
    lambda_cs_market=param(N*(1+K+size(mim_proj_var,2))+1:N*(1+K+size(mim_proj_var,2))+1);
    lambda_cs_pagm=param(N*(1+K+size(mim_proj_var,2))+2:N*(1+K+size(mim_proj_var,2))+2);
    
    lambda_ts_market=param(N*(1+K+size(mim_proj_var,3))+1:N*(1+K+size(mim_proj_var,2))+3);
    lambda_ts_pagm=param(N*(1+K+size(mim_proj_var,2))+4:N*(1+K+size(mim_proj_var,2))+4);  
    
    lambda_cs_market_t_stat=t_stat(N*(1+K+size(mim_proj_var,2))+1:N*(1+K+size(mim_proj_var,2))+1);
    lambda_cs_pagm_t_stat=t_stat(N*(1+K+size(mim_proj_var,2))+2:N*(1+K+size(mim_proj_var,2))+2);
    
    lambda_ts_market_t_stat=t_stat(N*(1+K+size(mim_proj_var,2))+3:N*(1+K+size(mim_proj_var,2))+3);
    lambda_ts_pagm_t_stat=t_stat(N*(1+K+size(mim_proj_var,2))+4:N*(1+K+size(mim_proj_var,2))+4);
    
    % test for identity of lambdas
    % restriction matrix R
    R1=zeros(size(param,1),1);
    R2=zeros(size(param,1),1);
    R1(N*(1+K+size(mim_proj_var,2))+1,1)=1;
    R1(N*(1+K+size(mim_proj_var,2))+3,1)=-1;
    R2(N*(1+K+size(mim_proj_var,2))+2,1)=1;
    R2(N*(1+K+size(mim_proj_var,2))+4,1)=-1;
    R1=R1';
    R2=R2';
    % Wald test statistic
    Wald_market=(R1*param)'/(R1*(Avar_param/size(mim_base_assets,1))*R1')*(R1*param);
    Wald_pagm=(R2*param)'/(R2*(Avar_param/size(mim_base_assets,1))*R2')*(R2*param);
    % p value
    pval_market=1-chi2cdf(Wald_market,1);
    pval_pagm=1-chi2cdf(Wald_pagm,1);
    
    
    disp(strcat('* Market lambda estimate  =  ',num2str(lambda_cs_market)));
    disp(strcat('* t-stat                  =  ',num2str(lambda_cs_market_t_stat)));
    disp(strcat('* cross-sectional estimate of lambda_N  =  ',num2str(lambda_cs_pagm)));
    disp(strcat('* t-stat                  =  ',num2str(lambda_cs_pagm_t_stat)));
    p_val_cs=(1-normcdf(abs(lambda_cs_pagm_t_stat)))*100*2;
    disp(strcat('* p-value in %            =  ',num2str(p_val_cs)));
    disp(strcat('* time series estimate of lambda_N      =  ',num2str(lambda_ts_pagm)));
    disp(strcat('* t-stat                  =  ',num2str(lambda_ts_pagm_t_stat)));
    p_val_ts=(1-normcdf(abs(lambda_ts_pagm_t_stat)))*100*2;
    disp(strcat('* p-value in %            =  ',num2str(p_val_ts)));    
    
    display('*************************************************');
    display('*  Wald-statistic to test identity of           *');
    display('*  two estimates (Eqn(13)) - p-value in %       *');
    disp(pval_pagm*100);
    display('*************************************************');    
    
    [ g_T_var, par_var ] = var_GMM( param, analysis, excess_returns, factors, factor_type, c_fac, mim_proj_var, mim_base_assets );

    [ betalambda ] = get_beta_lambda( param, analysis, excess_returns, factors, c_fac );
    
    betalambda_market=betalambda(1,:);
    betalambda_pagm=betalambda(2,:);
    
    [ beta_lambda_vec_vc, se_beta_lambda ] = beta_lambda_vc( param, par_var, analysis, excess_returns, factors, c_fac);
    
    se_market=se_beta_lambda(1,:);
    se_pagm=se_beta_lambda(2,:);
    
    betalambda_t_stat_market=betalambda_market./se_market;
    betalambda_t_stat_pagm=betalambda_pagm./se_pagm;
    
    disp(' gamma estimates implied by lambda estimates (Eqn(12)) ');
    
    % extract lambda_W and lambda_N from parameter vector: lambdas_mat
    lambdas_mat=param(N*(1+K+size(mim_proj_var,2))+1:N*(1+K+size(mim_proj_var,2))+2);
    % extract the mimicking portfolio parameters
    mim_par=param(N*(c_fac+1)+2*c_fac+size(mim_proj_var,2)+1:N*(c_fac+1)+2*c_fac+size(mim_proj_var,2)+size(mim_base_assets,2));
    % construct the vector of factors: market return and mimicking pf
    mim_pag=mim_base_assets*mim_par;
    factors_mat=[factors,mim_pag];
    % extract the variance covariance sub-matrix of the lambdas
    vc_lambda=par_var(N*(1+K+size(mim_proj_var,2))+1:N*(1+K+size(mim_proj_var,2))+2,N*(1+K+size(mim_proj_var,2))+1:N*(1+K+size(mim_proj_var,2))+2);
    
    % compute implied gamma representation parameters
    gamma_vec=(repmat((1/T),[size(factors_mat,2),size(factors_mat,2)]).*(factors_mat'*factors_mat)-mean(factors_mat)'*mean(factors_mat))\lambdas_mat;
    % VC matrix as described in Appendix A.1
    vc_gamma=(repmat((1/T),[size(factors_mat,2),size(factors_mat,2)]).*(factors_mat'*factors_mat)-mean(factors_mat)'*mean(factors_mat))\(vc_lambda/size(factors_mat,1))/(repmat((1/T),[size(factors_mat,2),size(factors_mat,2)]).*(factors_mat'*factors_mat)-mean(factors_mat)'*mean(factors_mat));
    % standard errors of recovered parameters
    se_gamma=sqrt(diag(vc_gamma));
    
    display('*************************************************');
    display('*  implied gamma estimates and t-statistics     *');
    display('*  (Eqn(12))                                    *');
    disp([gamma_vec,gamma_vec./se_gamma]);
    display('*************************************************'); 
    
    
    if testassets==1 || testassets==2
        % Next, for test assets 1 and 2 (VW and EW) specifically for the 
        % 25 pf, output for beta*lambda table is generated
        
        % Test the significance of a 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 pf is 5.
        % Additionally, we pick the second row, that is the PAG factor
        % NOTE: for the CDR model, PG must be the second factor!
        
        disp('Differences in expected returns analysis for MP-CDRM');
        asset_1=5;          % small value
        asset_2=21;         % big growth
        
        [ beta_lambda_diff_r, beta_lambda_se, beta_lambda_pval ] = get_beta_lambda_diff( betalambda, beta_lambda_vec_vc, 2, asset_1, asset_2, excess_returns );
        
        disp('reported in paper');
        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('**************************************************');
        display('*  Anaylsis for asset 5, small value, against 21 *');
        display('*  big growth                                    *');   
    
        % reshape into matrices
        betalambda_market=reshape(betalambda_market,[5,5]);
        betalambda_pagm=reshape(betalambda_pagm,[5,5]);
        betalambda_t_stat_market=reshape(betalambda_t_stat_market,[5,5]);
        betalambda_t_stat_pagm=reshape(betalambda_t_stat_pagm,[5,5]);
        
        display('**************************************************');
        display('*  Beta*lambda for MP-CDRM                       *');
        display('*  horizontal: Book-to-Market, vertical: size    *');
        disp(betalambda_pagm');
        display('**************************************************');
        display('*  t-stats beta_N*lambda_N for MP-CDRM           *'); 
        disp(betalambda_t_stat_pagm');
        display('**************************************************');
        display('*  paper: large-growth and small-value (corners) *');
        display('*  p-value (%): beta_N*lambda_N for MP-CDRM      *');   
        disp(((1-normcdf(abs(betalambda_t_stat_pagm)))*100*2)');
        display('**************************************************');
        
        if freq==2
           disp('annualized compensation: (((beta_N*lambda_N+^)^4)-1)*100');
           disp((((betalambda_pagm+1)^4)-1)*100);
           disp('annualized beta_N*lambda_N*4*100');
           disp(betalambda_pagm*4*100);
        end
    
    else
       disp('not implemented for test assets other than 25 B/M sorted pf'); 
    end
    
    
    %% goodness of fit graphs:
    % compute mean realized excess returns and average predicted returns
    
    mexret=mean(excess_returns)';
    model_implied_expected_return=sum(betalambda);
    
    % Maio and Santa Clara type Rē
    mexret_fit=model_implied_expected_return';
    alpha_hat=mexret-mexret_fit;
    alpha_hat_demeaned=alpha_hat-mean(alpha_hat)';
    mexret_demeaned=mexret-mean(mexret)';
    
    rsq=1-((sum(alpha_hat_demeaned.^2))/(sum(mexret_demeaned.^2)));
    
    display('**************************************************');
    display('*  cross-sectional R-squared  (%) (Eqn(15))      *');
    disp(rsq*100);
    display('**************************************************');  
    MAE=mean(abs(alpha_hat));
    display('*  Mean absolute pricing error (%) (Eqn(14))     *');    
    disp(MAE*100);
    display('**************************************************');
    
    
    % 45 degree plot
    
    % get a 45 degree line: auxiliary variables

    max_deg=max([max(model_implied_expected_return),max(mean(excess_returns))])+0.005;


    h2=figure(2);
    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(h2,'fortyfive.eps');
    
    %% J-test: validity of asset pricing model (Eqn(16))
    
    [ J_last_N, p_value_J_last_N, g_T_last_N, t_stat_last_N, corr_g_T, var_g_T_last_N ] = Jstat_GRS( param, analysis, excess_returns, factors, factor_type, c_fac, mim_proj_var, mim_base_assets );
    
    display('**************************************************');
    disp(strcat('J-statistic (Eqn(16))  =  ',num2str(J_last_N)));
    disp(strcat('p-value in %           =  ',num2str(p_value_J_last_N)));
    display('**************************************************');
    display('*  Average pricing errors (g_T) and t-stats:     *');
    disp([g_T_last_N,t_stat_last_N]);
    display('**************************************************');
    
    
end

