Improve debugging options for perfect foresight

We occasionally have issues with singular Jacobians as in https://forum.dynare.org/t/blanchard-conditions-in-stochastic-and-deterministic-models/15980 It would be good to have a debugging options that would point to the culprit for e.g. a row or column of zeros. Maybe we could add this to model_diagnostics. Two test cases are attached, one for sim1 and one for sim1_linear.

ADS_NE.mod code_14_march.mod

added enhancement perfect foresight labels

assigned to @normann

A piece of code in model_diagnostics carries out this check for the static jacobian. It is reproduced below:

try
    rank_jacob = rank(jacob); %can sometimes fail
catch
    rank_jacob=size(jacob,1);
end
if rank_jacob < size(jacob,1)
    problem_dummy=1;
    singularity_problem = 1;
    disp(['MODEL_DIAGNOSTICS:  The Jacobian of the static model is ' ...
            'singular'])
    disp(['MODEL_DIAGNOSTICS:  there is ' num2str(endo_nbr-rank_jacob) ...
            ' colinear relationships between the variables and the equations'])
    ncol = null(jacob);
    n_rel = size(ncol,2);
    for i = 1:n_rel
        if n_rel  > 1
            disp(['Relation ' int2str(i)])
        end
        disp('Colinear variables:')
        for j=1:10
            k = find(abs(ncol(:,i)) > 10^-j);
            if max(abs(jacob(:,k)*ncol(k,i))) < 1e-6
                break
            end
        end
        fprintf('%s\n',endo_names{k})
    end
    neq = null(jacob');
    n_rel = size(neq,2);
    for i = 1:n_rel
        if n_rel  > 1
            disp(['Relation ' int2str(i)])
        end
        disp('Colinear equations')
        for j=1:10
            k = find(abs(neq(:,i)) > 10^-j);
            if max(abs(jacob(k,:)'*neq(k,i))) < 1e-6
                break
            end
        end
        disp(k')
    end
end

A few remarks, though:

1. I do not understand the point of writing

for j=1:10
    k = find(abs(ncol(:,i)) > 10^-j);
    if max(abs(jacob(:,k)*ncol(k,i))) < 1e-6
        break
    end
end

instead of

k = find(abs(ncol(:,i)) > 10^(-10))

2. The rank function and the null function both use a SVD. Calling the null function is enough as one may infer the rank from the dimension of the kernel, and saves one SVD.

The remarks above may stem from a failure in understanding subtle things. Any enlightenment would be appreciated in that regard.

I propose to create a function that returns the rank and the dependent colums of a matrix. It would be as follow:

function [rank, rel] = lin_dep_col(A)
    ker = null(A);
    d = size(ker, 2);
    n = size(A, 2);
    rank = n-d;
    rel = {};
    if d > 0 
        for i = 1:d
            rel{end+1} = find(abs(ker(:, i)) > 10^(-10));
        end
    end

It could be useful for identification issues as well (e.g. dependent columns in a information matrix to detect identification issues as Iskrev does)

Adapting the piece of code above for the perfect foresight problems is not straightforward, though. The involved jacobian in the perfect foresight problem is the stacked jacobian and its link with the basic dynamic jacobian .dynamic files compute is not theoretically clear.

It does not seem natural to check the stacked jacobian in model_diagnostics.m, as it considers the basic dynamic jacobian. A solution could be to check the stacked jacobian in each method, as the considered jacobian may differ in structure (e.g. perfect-foresight-models/sim1_linear), and return the incriminated variables and equations.

What do you think?

1. @MichelJuillard may know the answer.
2. model_diagnostics is not called repeatedly, so efficiency was not a major focus. That being said, factorization is always a good idea as it makes maintaining codes easier and creates efficient code for reuse. We should maybe benchmark whether the approach taken in identification_checks.m is more efficient.
3. It's a good question where to put such a check. model_diagnostics has no way of knowing the context, i.e. whether we are doing perfect foresight simulations. Checking the stacked Jacobian there does not seem to be a good idea.
4. I agree with @normann that it would better fit in the respective solvers. As a starting point, we could simply add such a check for the first iteration in an if options_.debug-clause in sim1 and sim1_linear

1. The idea was to display first the variables with the highest coefficients. It would be a good thing to separate computing the dependent columns and displaying the relevent information.
2. The stacked Jacobian is made with dynamic jacobians, but the determinacy of the stacked jacobian can't be deduced from the determinacy of the dynamic jacobian. Furthermore dynamic jacobian are rectangular matrices.
3. It is necessary to look at the null space of the sparse stacked jacobian. The most interesting information is with the equations or variables involved in the null space near the beginning or near the end of the stacked jacobian matrix.

changed the description

mentioned in commit JohannesPfeifer/dynare@b5a4c069

mentioned in merge request !1834 (merged)

mentioned in commit a42a9fff

closed with merge request !1834 (merged)

closed with commit b5a4c069

mentioned in commit JohannesPfeifer/dynare@7c589667