function test = glejser_test(result_ols,Z,P)
%  PURPOSE: Breusch-Pagan/Godfrey test 
% ---------------------------------------------------------------------
%  USAGE: test = glejser_test(result_ols,Z,P)
%  where: result_ols = a structure returned by function OLS
%         Z = independent variables matrix Z including vector of ones
%                                           (nobs x nvar)   
%         P = level of significance (optional) by default P=0.95
% ---------------------------------------------------------------------
%  RETURNS: a structure
%         test.meth     = 'Glejser'
%         test.Fstat_e2 = F statistic for regression e^2
%         test.Fstat_abse = F statistic for regression abse
%         test.Fstat_logabse = F statistic for regression logabse
%         test.Fcrit    = critical value
%         test.result_e2      = 'hypothesis of homoscedasticity (not)
%                                rejected'
%         test.result_abse    = 'hypothesis of homoscedasticity (not) rejected'
%         test.result_logabse = 'hypothesis of homoscedasticity (not) rejected'
%         test.homosc_e2      = 1 (not rejected)or 0
%         test.homosc_abse    = 1 (not rejected)or 0
%         test.homosc_logabse = 1 (not rejected)or 0
%         test.result   = 'hypothesis of homoscedasticity (not) rejected'
%         test.homosc   = 1 (not rejected)or 0
% ---------------------------------------------------------------------

nIn = nargin;

if nIn < 2
  error('Not enough input arguments.');
elseif nIn > 3
  error('Too many input arguments.');
end

if nIn == 2
  P = 0.95;
end

%Glejser's test
test.meth = 'Glejser';
    [nobs,nvar] = size(Z);
    
    J = nvar-1;
    R = [zeros(J,1) eye(J)];
    q = zeros(J,1);
    
    
    e2 = result_ols.resid.^2;
    abse = abs(result_ols.resid);
    logabse = log(abse);
    
    res_e2 = ols(e2,Z);
    res_abse = ols(abse,Z);
    res_logabse = ols(logabse,Z);

        b=res_e2.beta;
    test.Fstat_e2 = (R*b-q)'*(R*(Z'*Z)^(-1)*R')^(-1)*(R*b-q)/J;
        test.Fstat_e2 = test.Fstat_e2/(res_e2.resid'*res_e2.resid/(nobs-nvar));    
    
        b=res_abse.beta;
    test.Fstat_abse = (R*b-q)'*(R*(Z'*Z)^(-1)*R')^(-1)*(R*b-q)/J;
        test.Fstat_abse = test.Fstat_abse/(res_abse.resid'*res_abse.resid/(nobs-nvar));    
    
        b=res_logabse.beta;
    test.Fstat_logabse = (R*b-q)'*(R*(Z'*Z)^(-1)*R')^(-1)*(R*b-q)/J;
        test.Fstat_logabse = test.Fstat_logabse/(res_logabse.resid'*res_logabse.resid/(nobs-nvar));   
    
    test.Fcrit = finv(P,J,nobs-nvar);

    if test.Fstat_e2 < test.Fcrit
        test.result_e2 = 'hypothesis of homoscedasticity NOT REJECTED';
        test.homosc_e2 = 1; 
    else
        test.result_e2 = 'hypothesis of homoscedasticity REJECTED';
        test.homosc_e2 = 0;
    end
    
    if test.Fstat_abse < test.Fcrit
        test.result_abse = 'hypothesis of homoscedasticity NOT REJECTED';
        test.homosc_abse = 1; 
    else
        test.result_abse = 'hypothesis of homoscedasticity REJECTED';
        test.homosc_abse = 0;
    end
    
    if test.Fstat_logabse < test.Fcrit
        test.result_logabse = 'hypothesis of homoscedasticity NOT REJECTED';
        test.homosc_logabse = 1; 
    else
        test.result_logabse = 'hypothesis of homoscedasticity REJECTED';
        test.homosc_logabse = 0;
    end

%SUMMARY
    if ((test.Fstat_e2 < test.Fcrit) + (test.Fstat_abse < test.Fcrit) ...
            + (test.Fstat_logabse < test.Fcrit)) > 1
        test.result = 'hypothesis of homoscedasticity NOT REJECTED';
        test.homosc = 1; 
    else
        test.result = 'hypothesis of homoscedasticity REJECTED';
        test.homosc = 0;
    end