function p = mbic(y, reg, x, maxp)
% (V)AR lag length selection by BIC
% INPUT: y is regressand nxm series
%        reg is an nxl matrix of variables that are included in the
%        regression, yet they are not lagged
%        x is nxk series
%        maxp is the max n. of lags
% OUTPUT: optimal lag length p
% based on Hayashi, page 398

[T m] = size(y);
[T k] = size(x); 
[T l] = size(reg);
T = T - maxp + 1;
ymat = y(maxp:end,:);
z = reg(maxp:end,:);

for i = 1:maxp
    z = [z x(maxp-i+1:end-i+1,:)];
end

i = 1;
while i <= maxp
    zmat = z(:,1:l+i*k);
    b = pinv(zmat'*zmat)*(zmat'*ymat);
    v = ymat - zmat*b;
    crit = log(det(v'*v/T))+(i*k + l)*m*log(T)/T; 
    if i <= 1; p = i; mincrit = crit; 
    elseif crit < mincrit; p = i; mincrit = crit; 
    end
    i=i+1;
end