adding new functions and new tests for risky steady state

matlab/AIM_first_order_solver.m

+function [dr,info]=AIM_first_order_solver(jacobia,M_,dr)
+    try
+        [dr,aimcode]=dynAIMsolver1(jacobia_,M_,dr);
+
+        % reuse some of the bypassed code and tests that may be needed 
+        
+        if aimcode ~=1
+            info(1) = aimcode;
+            info(2) = 1.0e+8;
+            return
+        end
+        [A,B] =transition_matrix(dr);
+        dr.eigval = eig(A);
+        nba = nd-sdim;
+
+        nyf = sum(kstate(:,2) > M_.maximum_endo_lag+1);
+
+        if nba ~= nyf
+            temp = sort(abs(dr.eigval));
+            if nba > nyf
+                temp = temp(nd-nba+1:nd-nyf)-1-options_.qz_criterium;
+                info(1) = 3;
+            elseif nba < nyf;
+                temp = temp(nd-nyf+1:nd-nba)-1-options_.qz_criterium;
+                info(1) = 4;
+            end
+            info(2) = temp'*temp;
+            return
+        end
+    catch
+        disp(lasterror.message)
+        error('Problem with AIM solver - Try to remove the "aim_solver" option')
+    end

matlab/dyn_first_order_solver.m

+function [dr,info] = dyn_first_order_solver(jacobia,b,M_,dr,options,task)
+    
+    info = 0;
+    
+    dr.ghx = [];
+    dr.ghu = [];
+    
+    klen = M_.maximum_endo_lag+M_.maximum_endo_lead+1;
+    kstate = dr.kstate;
+    kad = dr.kad;
+    kae = dr.kae;
+    nstatic = dr.nstatic;
+    nfwrd = dr.nfwrd;
+    npred = dr.npred;
+    nboth = dr.nboth;
+    order_var = dr.order_var;
+    nd = size(kstate,1);
+    lead_lag_incidence = M_.lead_lag_incidence;
+    nz = nnz(lead_lag_incidence);
+
+    sdyn = M_.endo_nbr - nstatic;
+
+    [junk,cols_b,cols_j] = find(lead_lag_incidence(M_.maximum_endo_lag+1,...
+                                                   order_var));
+    
+    if nstatic > 0
+        [Q,R] = qr(b(:,1:nstatic));
+        aa = Q'*jacobia;
+    else
+        aa = jacobia;
+    end
+    k1 = find([1:klen] ~= M_.maximum_endo_lag+1);
+    a = aa(:,nonzeros(lead_lag_incidence(k1,:)'));
+    b(:,cols_b) = aa(:,cols_j);
+    b10 = b(1:nstatic,1:nstatic);
+    b11 = b(1:nstatic,nstatic+1:end);
+    b2 = b(nstatic+1:end,nstatic+1:end);
+    if any(isinf(a(:)))
+        info = 1;
+        return
+    end
+
+    % buildind D and E
+    d = zeros(nd,nd) ;
+    e = d ;
+
+    k = find(kstate(:,2) >= M_.maximum_endo_lag+2 & kstate(:,3));
+    d(1:sdyn,k) = a(nstatic+1:end,kstate(k,3)) ;
+    k1 = find(kstate(:,2) == M_.maximum_endo_lag+2);
+    e(1:sdyn,k1) =  -b2(:,kstate(k1,1)-nstatic);
+    k = find(kstate(:,2) <= M_.maximum_endo_lag+1 & kstate(:,4));
+    e(1:sdyn,k) = -a(nstatic+1:end,kstate(k,4)) ;
+    k2 = find(kstate(:,2) == M_.maximum_endo_lag+1);
+    k2 = k2(~ismember(kstate(k2,1),kstate(k1,1)));
+    d(1:sdyn,k2) = b2(:,kstate(k2,1)-nstatic);
+
+    if ~isempty(kad)
+        for j = 1:size(kad,1)
+            d(sdyn+j,kad(j)) = 1 ;
+            e(sdyn+j,kae(j)) = 1 ;
+        end
+    end
+
+    % 1) if mjdgges.dll (or .mexw32 or ....) doesn't exit, 
+    % matlab/qz is added to the path. There exists now qz/mjdgges.m that 
+    % contains the calls to the old Sims code 
+    % 2) In  global_initialization.m, if mjdgges.m is visible exist(...)==2, 
+    % this means that the DLL isn't avaiable and use_qzdiv is set to 1
+    
+    [err,ss,tt,w,sdim,dr.eigval,info1] = mjdgges(e,d,options.qz_criterium);
+    mexErrCheck('mjdgges', err);
+
+    if info1
+        if info1 == -30
+            info(1) = 7;
+        else
+            info(1) = 2;
+            info(2) = info1;
+            info(3) = size(e,2);
+        end
+        return
+    end
+
+    nba = nd-sdim;
+
+    nyf = sum(kstate(:,2) > M_.maximum_endo_lag+1);
+
+    if task == 1
+        dr.rank = rank(w(1:nyf,nd-nyf+1:end));
+        % Under Octave, eig(A,B) doesn't exist, and
+        % lambda = qz(A,B) won't return infinite eigenvalues
+        if ~exist('OCTAVE_VERSION')
+            dr.eigval = eig(e,d);
+        end
+        return
+    end
+
+    if nba ~= nyf
+        temp = sort(abs(dr.eigval));
+        if nba > nyf
+            temp = temp(nd-nba+1:nd-nyf)-1-options.qz_criterium;
+            info(1) = 3;
+        elseif nba < nyf;
+            temp = temp(nd-nyf+1:nd-nba)-1-options.qz_criterium;
+            info(1) = 4;
+        end
+        info(2) = temp'*temp;
+        return
+    end
+
+    np = nd - nyf;
+    n2 = np + 1;
+    n3 = nyf;
+    n4 = n3 + 1;
+    % derivatives with respect to dynamic state variables
+    % forward variables
+    w1 =w(1:n3,n2:nd);
+    if ~isscalar(w1) && (condest(w1) > 1e9);
+        % condest() fails on a scalar under Octave
+        info(1) = 5;
+        info(2) = condest(w1);
+        return;
+    else
+        gx = -w1'\w(n4:nd,n2:nd)';
+    end  
+
+    % predetermined variables
+    hx = w(1:n3,1:np)'*gx+w(n4:nd,1:np)';
+    hx = (tt(1:np,1:np)*hx)\(ss(1:np,1:np)*hx);
+
+    k1 = find(kstate(n4:nd,2) == M_.maximum_endo_lag+1);
+    k2 = find(kstate(1:n3,2) == M_.maximum_endo_lag+2);
+    dr.ghx = [hx(k1,:); gx(k2(nboth+1:end),:)];
+
+    %lead variables actually present in the model
+    j3 = nonzeros(kstate(:,3));
+    j4  = find(kstate(:,3));
+    % derivatives with respect to exogenous variables
+    if M_.exo_nbr
+        fu = aa(:,nz+(1:M_.exo_nbr));
+        a1 = b;
+        aa1 = [];
+        if nstatic > 0
+            aa1 = a1(:,1:nstatic);
+        end
+        dr.ghu = -[aa1 a(:,j3)*gx(j4,1:npred)+a1(:,nstatic+1:nstatic+ ...
+                                                 npred) a1(:,nstatic+npred+1:end)]\fu;
+    else
+        dr.ghu = [];
+    end
+
+    % static variables
+    if nstatic > 0
+        temp = -a(1:nstatic,j3)*gx(j4,:)*hx;
+        j5 = find(kstate(n4:nd,4));
+        temp(:,j5) = temp(:,j5)-a(1:nstatic,nonzeros(kstate(:,4)));
+        temp = b10\(temp-b11*dr.ghx);
+        dr.ghx = [temp; dr.ghx];
+        temp = [];
+    end
+
+    if options.use_qzdiv
+        %% Necessary when using Sims' routines for QZ
+        gx = real(gx);
+        hx = real(hx);
+        dr.ghx = real(dr.ghx);
+        dr.ghu = real(dr.ghu);
+    end
+
+    dr.Gy = hx;
\ No newline at end of file

matlab/dyn_ramsey_linearized_foc.m

+function [jacobia_,dr,info,M_,oo_] = dyn_ramsey_linearized_foc(dr,M_,options_,oo_)
+% function [jacobia_,dr,info,M_,oo_] = dyn_ramsey_linearized_foc(dr,M_,options_,oo_)
+% computes the Jacobian of the linear approximation of the F.O.C of a
+% Ramsey problem
+%
+% INPUTS
+%   dr         [matlab structure] Decision rules for stochastic simulations.
+%   M_         [matlab structure] Definition of the model.           
+%   options_   [matlab structure] Global options.
+%   oo_        [matlab structure] Results 
+%    
+% OUTPUTS
+%   dr         [matlab structure] Decision rules for stochastic simulations.
+%   info       [integer]          info=1: the model doesn't define current variables uniquely
+%                                 info=2: problem in mjdgges.dll info(2) contains error code. 
+%                                 info=3: BK order condition not satisfied info(2) contains "distance"
+%                                         absence of stable trajectory.
+%                                 info=4: BK order condition not satisfied info(2) contains "distance"
+%                                         indeterminacy.
+%                                 info=5: BK rank condition not satisfied.
+%                                 info=6: The jacobian matrix evaluated at the steady state is complex.        
+%   M_         [matlab structure]            
+%   options_   [matlab structure]
+%   oo_        [matlab structure]
+%  
+% ALGORITHM
+%   ...
+%    
+% SPECIAL REQUIREMENTS
+%   none.
+%  
+
+% Copyright (C) 1996-2009 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+info = 0;
+
+if isfield(M_,'orig_model')
+        orig_model = M_.orig_model;
+        M_.endo_nbr = orig_model.endo_nbr;
+        M_.orig_endo_nbr = orig_model.orig_endo_nbr;
+        M_.aux_vars = orig_model.aux_vars;
+        M_.endo_names = orig_model.endo_names;
+        M_.lead_lag_incidence = orig_model.lead_lag_incidence;
+        M_.maximum_lead = orig_model.maximum_lead;
+        M_.maximum_endo_lead = orig_model.maximum_endo_lead;
+        M_.maximum_lag = orig_model.maximum_lag;
+        M_.maximum_endo_lag = orig_model.maximum_endo_lag;
+    end
+
+    if options_.steadystate_flag
+        k_inst = [];
+        instruments = options_.instruments;
+        for i = 1:size(instruments,1)
+            k_inst = [k_inst; strmatch(options_.instruments(i,:), ...
+                                       M_.endo_names,'exact')];
+        end
+        ys = oo_.steady_state;
+        [inst_val,info1] = dynare_solve('dyn_ramsey_static_', ...
+                                oo_.steady_state(k_inst),0, ...
+                                M_,options_,oo_,it_);
+        M_.params = evalin('base','M_.params;');
+        ys(k_inst) = inst_val;
+        [x,check] = feval([M_.fname '_steadystate'],...
+                          ys,[oo_.exo_steady_state; ...
+                            oo_.exo_det_steady_state]);
+        if size(x,1) < M_.endo_nbr 
+            if length(M_.aux_vars) > 0
+                x = add_auxiliary_variables_to_steadystate(x,M_.aux_vars,...
+                                                           M_.fname,...
+                                                           oo_.exo_steady_state,...
+                                                           oo_.exo_det_steady_state,...
+                                                           M_.params);
+            else
+                error([M_.fname '_steadystate.m doesn''t match the model']);
+            end
+        end
+        oo_.steady_state = x;
+        [junk,junk,multbar] = dyn_ramsey_static_(oo_.steady_state(k_inst),M_,options_,oo_,it_);
+    else
+        [oo_.steady_state,info1] = dynare_solve('dyn_ramsey_static_', ...
+                                        oo_.steady_state,0,M_,options_,oo_,it_);
+        [junk,junk,multbar] = dyn_ramsey_static_(oo_.steady_state,M_,options_,oo_,it_);
+    end
+        
+    check1 = max(abs(feval([M_.fname '_static'],...
+                           oo_.steady_state,...
+                           [oo_.exo_steady_state; ...
+                        oo_.exo_det_steady_state], M_.params))) > options_.dynatol ;
+    if check1
+        info(1) = 20;
+        info(2) = check1'*check1;
+        return
+    end
+    
+    [jacobia_,M_] = dyn_ramsey_dynamic_(oo_.steady_state,multbar,M_,options_,oo_,it_);
+    klen = M_.maximum_lag + M_.maximum_lead + 1;
+    dr.ys = [oo_.steady_state;zeros(M_.exo_nbr,1);multbar];
+    oo_.steady_state = dr.ys;
+    
+    if options_.noprint == 0
+        disp_steady_state(M_,oo_)
+        for i=M_.orig_endo_nbr:M_.endo_nbr
+            if strmatch('mult_',M_.endo_names(i,:))
+                disp(sprintf('%s \t\t %g',M_.endo_names(i,:), ...
+                             dr.ys(i)));
+            end
+        end
+    end
+

matlab/dyn_risky_steadystate_solver.m

+function [dr,info] = dyn_risky_steadystate_solver(ys0,M, ...
+                                                  dr,options,oo)
+    
+    info = 0;
+    lead_lag_incidence = M.lead_lag_incidence;
+    order_var = dr.order_var;
+    exo_nbr = M.exo_nbr;
+    
+    M.var_order_endo_names = M.endo_names(dr.order_var,:);
+ 
+    [junk,dr.i_fwrd_g,i_fwrd_f] = find(lead_lag_incidence(3,order_var));
+    dr.i_fwrd_f = i_fwrd_f;
+    nd = nnz(lead_lag_incidence) + M.exo_nbr;
+    dr.nd = nd;
+    kk = reshape(1:nd^2,nd,nd);
+    kkk = reshape(1:nd^3,nd^2,nd);
+    dr.i_fwrd2_f = kk(i_fwrd_f,i_fwrd_f);
+    dr.i_fwrd2a_f = kk(i_fwrd_f,:);
+    dr.i_fwrd3_f = kkk(dr.i_fwrd2_f,:);
+    dr.i_uu = kk(end-exo_nbr+1:end,end-exo_nbr+1:end);
+    if options.k_order_solver
+        func = @risky_residuals_k_order;
+    else
+        func = @risky_residuals;
+    end
+    
+    if isfield(options,'portfolio') && options.portfolio == 1
+        eq_tags = M.equations_tags;
+        n_tags = size(eq_tags,1);
+        l_var = zeros(n_tags,1);
+        for i=1:n_tags
+            l_var(i) = find(strncmp(eq_tags(i,3),M.endo_names, ...
+                                    length(cell2mat(eq_tags(i,3)))));
+        end
+        dr.ys = ys0;
+        x0 = ys0(l_var);
+        %        dr = first_step_ds(x0,M,dr,options,oo);
+        n = size(ys0);
+        %x0 = ys0;
+        [x, info] = solve1(@risky_residuals_ds,x0,1:n_tags,1:n_tags,0,1,M,dr,options,oo);
+        %[x, info] = solve1(@risky_residuals,x0,1:n,1:n,0,1,M,dr,options,oo);
+        %        ys0(l_var) = x;
+        ys0(l_var) = x;
+        dr.ys = ys0;
+        oo.dr = dr;
+        oo.steady_state = ys0;
+        disp_steady_state(M,oo);
+    end
+        
+    [ys, info] = csolve(func,ys0,[],1e-10,100,M,dr,options,oo);
+    
+    if options.k_order_solver
+        [resid,dr] = risky_residuals_k_order(ys,M,dr,options,oo);
+    else
+        [resid,dr] = risky_residuals(ys,M,dr,options,oo);
+    end
+    
+    dr.ys = ys;
+    for i=1:M.endo_nbr
+        disp(sprintf('%16s %12.6f %12.6f',M.endo_names(i,:),ys0(i), ys(i)))
+    end
+    
+    dr.ghs2 = zeros(size(dr.ghs2));
+
+    k_var = setdiff(1:M.endo_nbr,l_var);
+    dr.ghx(k_var,:) = dr.ghx;
+    dr.ghu(k_var,:) = dr.ghu;
+    dr.ghxx(k_var,:) = dr.ghxx;
+    dr.ghxu(k_var,:) = dr.ghxu;
+    dr.ghuu(k_var,:) = dr.ghuu;
+    dr.ghs2(k_var,:) = dr.ghs2;
+end
+
+function [resid,dr] = risky_residuals(ys,M,dr,options,oo)
+    persistent old_ys old_resid
+    
+    lead_lag_incidence = M.lead_lag_incidence;
+    iyv = lead_lag_incidence';
+    iyv = iyv(:);
+    iyr0 = find(iyv) ;
+    
+    if M.exo_nbr == 0
+        oo.exo_steady_state = [] ;
+    end
+    
+    z = repmat(ys,1,3);
+    z = z(iyr0) ;
+    [resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
+                                     [oo.exo_simul ...
+                        oo.exo_det_simul], M.params, dr.ys, 2);
+    if ~isreal(d1) || ~isreal(d2)
+        pause
+    end
+    
+    if options.use_dll
+        % In USE_DLL mode, the hessian is in the 3-column sparse representation
+        d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
+                         size(d1, 1), size(d1, 2)*size(d1, 2));
+    end
+
+    if isfield(options,'portfolio') && options.portfolio == 1
+        eq_tags = M.equations_tags;
+        n_tags = size(eq_tags,1);
+        portfolios_eq = cell2mat(eq_tags(strcmp(eq_tags(:,2), ...
+                                                'portfolio'),1));
+        eq = setdiff(1:M.endo_nbr,portfolios_eq);
+        l_var = zeros(n_tags,1);
+        for i=1:n_tags
+            l_var(i) = find(strncmp(eq_tags(i,3),M.endo_names, ...
+                                    length(cell2mat(eq_tags(i,3)))));
+        end
+        k_var = setdiff(1:M.endo_nbr,l_var);
+        lli1 = lead_lag_incidence(:,k_var);
+        lead_incidence = lli1(3,:)';
+        k = find(lli1');
+        lli2 = lli1';
+        lli2(k) = 1:nnz(lli1);
+        lead_lag_incidence = lli2';
+        x = ys(l_var);
+        dr = first_step_ds(x,M,dr,options,oo);
+
+        
+        M.lead_lag_incidence = lead_lag_incidence;
+        lli1a = [nonzeros(lli1'); size(d1,2)+(-M.exo_nbr+1:0)'];
+        d1a = d1(eq,lli1a);
+        ih = 1:size(d2,2);
+        ih = reshape(ih,size(d1,2),size(d1,2));
+        ih1 = ih(lli1a,lli1a);
+        d2a = d2(eq,ih1);
+        
+        M.endo_nbr = M.endo_nbr-n_tags;
+        dr = set_state_space(dr,M);
+    
+        [junk,dr.i_fwrd_g] = find(lead_lag_incidence(3,dr.order_var));
+        i_fwrd_f = nonzeros(lead_incidence(dr.order_var));
+        i_fwrd2_f = ih(i_fwrd_f,i_fwrd_f);
+        dr.i_fwrd_f = i_fwrd_f;
+        dr.i_fwrd2_f = i_fwrd2_f;
+        nd = nnz(lead_lag_incidence) + M.exo_nbr;
+        dr.nd = nd;
+        kk = reshape(1:nd^2,nd,nd);
+        kkk = reshape(1:nd^3,nd^2,nd);
+        dr.i_fwrd2a_f = kk(i_fwrd_f,:);
+        %        dr.i_fwrd3_f = kkk(i_fwrd2_f,:);
+        dr.i_uu = kk(end-M.exo_nbr+1:end,end-M.exo_nbr+1:end);
+    else
+        d1a = d1;
+        d2a = d2;
+    end
+    
+% $$$     [junk,cols_b,cols_j] = find(lead_lag_incidence(2,dr.order_var));
+% $$$     b = zeros(M.endo_nbr,M.endo_nbr);
+% $$$     b(:,cols_b) = d1a(:,cols_j);
+% $$$ 
+% $$$     [dr,info] = dyn_first_order_solver(d1a,b,M,dr,options,0);
+% $$$     if info
+% $$$         [m1,m2]=max(abs(ys-old_ys));
+% $$$         disp([m1 m2])
+% $$$         %        print_info(info,options.noprint);
+% $$$         resid = old_resid+info(2)/40;
+% $$$         return
+% $$$     end
+% $$$     
+% $$$     dr = dyn_second_order_solver(d1a,d2a,dr,M);
+    
+    gu1 = dr.ghu(dr.i_fwrd_g,:);
+
+    resid = resid1+0.5*(d1(:,dr.i_fwrd_f)*dr.ghuu(dr.i_fwrd_g,:)+ ...
+                        d2(:,dr.i_fwrd2_f)*kron(gu1,gu1))*vec(M.Sigma_e);
+    disp(d1(:,dr.i_fwrd_f)*dr.ghuu(dr.i_fwrd_g,:)*vec(M.Sigma_e));
+    old_ys = ys;
+    disp(max(abs(resid)))
+    old_resid = resid;
+end
+
+function [resid,dr] = risky_residuals_ds(x,M,dr,options,oo)
+    persistent old_ys old_resid old_resid1 old_d1 old_d2
+    
+    dr = first_step_ds(x,M,dr,options,oo);
+
+    lead_lag_incidence = M.lead_lag_incidence;
+    iyv = lead_lag_incidence';
+    iyv = iyv(:);
+    iyr0 = find(iyv) ;
+    
+    if M.exo_nbr == 0
+        oo.exo_steady_state = [] ;
+    end
+    
+    eq_tags = M.equations_tags;
+    n_tags = size(eq_tags,1);
+    portfolios_eq = cell2mat(eq_tags(strcmp(eq_tags(:,2), ...
+                                            'portfolio'),1));
+    eq = setdiff(1:M.endo_nbr,portfolios_eq);
+    l_var = zeros(n_tags,1);
+    for i=1:n_tags
+        l_var(i) = find(strncmp(eq_tags(i,3),M.endo_names, ...
+                                length(cell2mat(eq_tags(i,3)))));
+    end
+    k_var = setdiff(1:M.endo_nbr,l_var);
+    lli1 = lead_lag_incidence(:,k_var);
+    k = find(lli1');
+    lli2 = lli1';
+    lli2(k) = 1:nnz(lli1);
+    lead_lag_incidence = lli2';
+
+    ys = dr.ys;
+    ys(l_var) = x;
+    
+    z = repmat(ys,1,3);
+    z = z(iyr0) ;
+    [resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
+                                     [oo.exo_simul ...
+                        oo.exo_det_simul], M.params, dr.ys, 2);
+% $$$     if isempty(old_resid)
+% $$$         old_resid1 = resid1;
+% $$$         old_d1 = d1;
+% $$$         old_d2 = d2;
+% $$$         old_ys = ys;
+% $$$     else
+% $$$         if ~isequal(resid1,old_resid)
+% $$$             disp('ys')
+% $$$             disp((ys-old_ys)');
+% $$$             disp('resids1')
+% $$$             disp((resid1-old_resid1)')
+% $$$             old_resid1 = resid1;
+% $$$             pause
+% $$$         end
+% $$$         if ~isequal(d1,old_d1)
+% $$$             disp('d1')
+% $$$             disp(d1-old_d1);
+% $$$             old_d1 = d1;
+% $$$             pause
+% $$$         end
+% $$$         if ~isequal(d2,old_d2)
+% $$$             disp('d2')
+% $$$             disp(d2-old_d2);
+% $$$             old_d2 = d2;
+% $$$             pause
+% $$$         end
+% $$$     end
+    if ~isreal(d1) || ~isreal(d2)
+        pause
+    end
+    
+    if options.use_dll
+        % In USE_DLL mode, the hessian is in the 3-column sparse representation
+        d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
+                         size(d1, 1), size(d1, 2)*size(d1, 2));
+    end
+
+% $$$     if isfield(options,'portfolio') && options.portfolio == 1
+% $$$         lli1a = [nonzeros(lli1'); size(d1,2)+(-M.exo_nbr+1:0)'];
+% $$$         d1a = d1(eq,lli1a);
+% $$$         ih = 1:size(d2,2);
+% $$$         ih = reshape(ih,size(d1,2),size(d1,2));
+% $$$         ih1 = ih(lli1a,lli1a);
+% $$$         d2a = d2(eq,ih1);
+% $$$         
+% $$$         M.endo_nbr = M.endo_nbr-n_tags;