function [ J_last_N, p_value_J_last_N, g_T_last_N, t_stats_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 )
%  Auxiliary function for the two pass regression in a_T*g_T=0 setup:
%  Computes J-test for validity of the moment conditions for a two pass
%  regression test of the empirical asset pricing model.

%  The last N moment conditions are the moment conditions implied by the
%  asset pricing model. We read out the last N elements of g_T and its
%  respective variance-covariance matrix and construct a Wald-type test.
%  The vc-matrix is not of full rank per construction as usual for the
%  J-test. According to Cochrane (2005) this is a less restrictive version
%  of the more familiar GRS statistic (we do not have to assume e.g.
%  independence of factors and residuals, etc.)

% Under the null hypothesis that the asset pricing model is correctly
% specified, the test is chi-square with degrees of freedom N (number of
% test assets) - K (number of factors).

% This function returns the J-statistic, its p_value, the last N g_T
% (sample average moment conditions), the t-stats for a test that each of
% them is zero, the correlation matrix and the variance-covariance matrix
% of the last N g_T.

N=size(excess_returns,2);
K=size(factors,2);

% compute all sample moment conditions
g_T_full=get_g_T( param, analysis, excess_returns, factors, factor_type, c_fac, mim_proj_var, mim_base_assets );
% read out last N elements
g_T_last_N=g_T_full(size(g_T_full,1)-size(excess_returns,2)+1:size(g_T_full,1));

% variance-covariance matrix of the g_T
[Avar_g_T,~]=var_GMM( param, analysis, excess_returns, factors, factor_type, c_fac, mim_proj_var, mim_base_assets );
% read out lower right block of moment conditions (average pricing errors)
Avar_g_T_last_N=Avar_g_T(size(g_T_full,1)-size(excess_returns,2)+1:size(g_T_full,1),size(g_T_full,1)-size(excess_returns,2)+1:size(g_T_full,1));

% compute correlation of average pricing errors
corr_g_T=diag(ones(size(Avar_g_T_last_N,1),1)./sqrt(diag(Avar_g_T_last_N)))*Avar_g_T_last_N*diag(ones(size(Avar_g_T_last_N,1),1)./sqrt(diag(Avar_g_T_last_N)));

% compute J-statistic and its p-value
J_last_N=sqrt(size(excess_returns,1))*g_T_last_N'*pinv(Avar_g_T_last_N)*sqrt(size(excess_returns,1))*g_T_last_N;

p_value_J_last_N=(1-chi2cdf(J_last_N,N-K))*100;

% compute t-statistic for each of the sample average moment conditions
t_stats_last_N=sqrt(size(excess_returns,1))*g_T_last_N./sqrt(diag(Avar_g_T_last_N));


end

