function test = BPKB_test(result_ols,Z,P)
%  PURPOSE: Breusch-Pagan test (Koenker and Basset variant)
% ---------------------------------------------------------------------
%  USAGE: test = BPKB_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 (Koenker-Basset)'
%         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 (Koenker-Basset)
test.meth = 'Breusch-Pagan (Koenker-Basset)';
    [nobs,nvar] = size(Z);
    u = result_ols.resid.^2;
    u_ = result_ols.resid'*result_ols.resid/nobs;
    V = mean((u-u_).^2);
        
    test.LM = 1/V*((u-u_)'*Z*(Z'*Z)^(-1)*Z'*(u-u_));
    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