function test = BP_test(result_ols,Z,P)
%  PURPOSE: Breusch-Pagan/Godfrey test 
% ---------------------------------------------------------------------
%  USAGE: test = BP_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     = 'Breusch-Pagan'
%         test.LM       = Lagrange multiplier statistics
%         test.chi2crit = critical value
%         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

%Breusch-Pagan test
test.meth = 'Breusch-Pagan';
    [nobs,nvar] = size(Z);
    g = (result_ols.resid.^2)./(result_ols.resid'*result_ols.resid/nobs)-1;
    
    test.LM = 1/2*(g'*Z*(Z'*Z)^(-1)*Z'*g);
    test.chi2crit = chi2inv(P,nvar-1);

    if test.LM < test.chi2crit
        test.result = 'hypothesis of homoscedasticity NOT REJECTED';
        test.homosc = 1; 
    else
        test.result = 'hypothesis of homoscedasticity REJECTED';
        test.homosc = 0;
    end