Inconsistency in treatment of qz_criterium in dyn_first_order_solver.m / mjdgges.c
I was under the impression that qz_criterium
was intended to be a threshold on the absolute value of dr.eigval
. See e.g. line 85 of AIM_first_order_solver.m
:
nba = nd-sum( abs(dr.eigval) < qz_criterium );
However, it appears that in dyn_first_order_solver.m
it is instead a threshold on the squared absolute value of dr.eigval
. See line 190 of dyn_first_order_solver.m
:
[err, ss, tt, w, sdim, dr.eigval, info1] = mjdgges(E, D, DynareOptions.qz_criterium, DynareOptions.qz_zero_threshold);
line 31 of mjdgges.c
:
return ((*alphar **alphar + *alphai **alphai) < criterium **beta **beta);
and line 107 of mjdgges.c
:
criterium = *mxGetPr(prhs[2]);
The eigenvalues are ( alphar[j] + alphai[j] * 1i ) / beta[j]
according to the MKL manual. Thus, it seems that this code is comparing the squared absolute value of the eigenvalue to qz_criterium
.
The easiest fix would be to replace line 190 of dyn_first_order_solver.m
with:
[err, ss, tt, w, sdim, dr.eigval, info1] = mjdgges(E, D, DynareOptions.qz_criterium .^ 2, DynareOptions.qz_zero_threshold);
so you're comparing squared to squared.