function blnStability = check_stability(vecBeta, intN, Nlags)


% Initialise:
blnStability = true;

% Get companion form
matBetaTop = reshape(vecBeta, intN*Nlags+1, intN)';
% Get rid of constants (last col.):
matBetaTop = matBetaTop(:, 1:end-1);
matBetaBottom  = eye(intN*(Nlags-1), intN*Nlags);
matBetaC = [matBetaTop; matBetaBottom];

% NOTE: used to be matBetaBottom   = eye(intN, intN * Nlags);
% But this seems wrong and does not work for #lags ~= 2.


% Get modulus of eigenvalues:
vecModulus = abs(eig(matBetaC));

% Check stability:
if max(vecModulus) > 1
    blnStability = false;
end