dr_block

 function [dr,info,M_,options_,oo_] = dr_block(dr,task,M_,options_,oo_)
 computes the reduced form solution of a rational expectation model (first
 approximation of the stochastic model around the deterministic steady state). 

 INPUTS
   dr         [matlab structure] Decision rules for stochastic simulations.
   task       [integer]          if task = 0 then dr1 computes decision rules.
                                 if task = 1 then dr1 computes eigenvalues.
   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
   first order block relaxation method applied to the model
    E[A Yt-1 + B Yt + C Yt+1 + ut] = 0
    
 SPECIAL REQUIREMENTS
   none.

0001 function [dr,info,M_,options_,oo_] = dr_block(dr,task,M_,options_,oo_)
0002 % function [dr,info,M_,options_,oo_] = dr_block(dr,task,M_,options_,oo_)
0003 % computes the reduced form solution of a rational expectation model (first
0004 % approximation of the stochastic model around the deterministic steady state).
0005 %
0006 % INPUTS
0007 %   dr         [matlab structure] Decision rules for stochastic simulations.
0008 %   task       [integer]          if task = 0 then dr1 computes decision rules.
0009 %                                 if task = 1 then dr1 computes eigenvalues.
0010 %   M_         [matlab structure] Definition of the model.
0011 %   options_   [matlab structure] Global options.
0012 %   oo_        [matlab structure] Results
0013 %
0014 % OUTPUTS
0015 %   dr         [matlab structure] Decision rules for stochastic simulations.
0016 %   info       [integer]          info=1: the model doesn't define current variables uniquely
0017 %                                 info=2: problem in mjdgges.dll info(2) contains error code.
0018 %                                 info=3: BK order condition not satisfied info(2) contains "distance"
0019 %                                         absence of stable trajectory.
0020 %                                 info=4: BK order condition not satisfied info(2) contains "distance"
0021 %                                         indeterminacy.
0022 %                                 info=5: BK rank condition not satisfied.
0023 %                                 info=6: The jacobian matrix evaluated at the steady state is complex.
0024 %   M_         [matlab structure]
0025 %   options_   [matlab structure]
0026 %   oo_        [matlab structure]
0027 %
0028 % ALGORITHM
0029 %   first order block relaxation method applied to the model
0030 %    E[A Yt-1 + B Yt + C Yt+1 + ut] = 0
0031 %
0032 % SPECIAL REQUIREMENTS
0033 %   none.
0034 %
0035 
0036 % Copyright (C) 2010-2011 Dynare Team
0037 %
0038 % This file is part of Dynare.
0039 %
0040 % Dynare is free software: you can redistribute it and/or modify
0041 % it under the terms of the GNU General Public License as published by
0042 % the Free Software Foundation, either version 3 of the License, or
0043 % (at your option) any later version.
0044 %
0045 % Dynare is distributed in the hope that it will be useful,
0046 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0047 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0048 % GNU General Public License for more details.
0049 %
0050 % You should have received a copy of the GNU General Public License
0051 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0052 
0053 info = 0;
0054 verbose = 0;
0055 if options_.order > 1
0056     error('2nd and 3rd order approximation not implemented with block option')
0057 end
0058 
0059 z = repmat(dr.ys,1,M_.maximum_lead + M_.maximum_lag + 1);
0060 zx = repmat([oo_.exo_simul oo_.exo_det_simul],M_.maximum_lead + M_.maximum_lag + 1, 1);
0061 if (isfield(M_,'block_structure'))
0062     data = M_.block_structure.block;
0063     Size = length(M_.block_structure.block);
0064 else
0065     data = M_;
0066     Size = 1;
0067 end;
0068 if (options_.bytecode)
0069     [chck, zz, data]= bytecode('dynamic','evaluate', z, zx, M_.params, dr.ys, 1, data);
0070 else
0071     [r, data] = feval([M_.fname '_dynamic'], z', zx, M_.params, dr.ys, M_.maximum_lag+1, data);
0072     chck = 0;
0073 end;
0074 mexErrCheck('bytecode', chck);
0075 dr.rank = 0;
0076 dr.eigval = [];
0077 dr.nstatic = 0;
0078 dr.nfwrd = 0;
0079 dr.npred = 0;
0080 dr.nboth = 0;
0081 dr.nd = 0;
0082 
0083 dr.ghx = [];
0084 dr.ghu = [];
0085 %Determine the global list of state variables:
0086 dr.state_var = M_.state_var;
0087 M_.block_structure.state_var = dr.state_var;
0088 n_sv = size(dr.state_var, 2);
0089 dr.ghx = zeros(M_.endo_nbr, length(dr.state_var));
0090 dr.exo_var = 1:M_.exo_nbr;
0091 dr.ghu = zeros(M_.endo_nbr, M_.exo_nbr);
0092 dr.nstatic = M_.nstatic;
0093 dr.nfwrd = M_.nfwrd;
0094 dr.npred = M_.npred;
0095 dr.nboth = M_.nboth;
0096 dr.nyf = 0;
0097 for i = 1:Size;
0098     ghx = [];
0099     indexi_0 = 0;
0100     if (verbose)
0101         disp('======================================================================');
0102         disp(['Block ' int2str(i)]);
0103         disp('-----------');
0104         data(i)
0105     end;
0106     n_pred = data(i).n_backward;
0107     n_fwrd = data(i).n_forward;
0108     n_both = data(i).n_mixed;
0109     n_static = data(i).n_static;
0110     dr.nyf = dr.nyf  + n_fwrd + n_both;
0111     nd = n_pred + n_fwrd + 2*n_both;
0112     dr.nd = dr.nd + nd;
0113     n_dynamic = n_pred + n_fwrd + n_both;
0114     exo_nbr = M_.block_structure.block(i).exo_nbr;
0115     exo_det_nbr = M_.block_structure.block(i).exo_det_nbr;
0116     other_endo_nbr = M_.block_structure.block(i).other_endo_nbr;
0117     jacob = full(data(i).g1);
0118     lead_lag_incidence = data(i).lead_lag_incidence;
0119     endo = data(i).variable;
0120     exo = data(i).exogenous;
0121     if (verbose)
0122         disp('jacob');
0123         disp(jacob);
0124         disp('lead_lag_incidence');
0125         disp(lead_lag_incidence);
0126     end;
0127     maximum_lag = data(i).maximum_endo_lag;
0128     maximum_lead = data(i).maximum_endo_lead;
0129     n = n_dynamic + n_static;
0130     
0131     block_type = M_.block_structure.block(i).Simulation_Type;
0132     if task ~= 1
0133         if block_type == 2 || block_type == 4 || block_type == 7 
0134             block_type = 8;
0135         end;
0136     end;
0137     if maximum_lag > 0 && (n_pred > 0  || n_both > 0) && block_type ~= 1 
0138         indexi_0 = min(lead_lag_incidence(2,:));
0139     end;
0140     switch block_type
0141       case 1
0142       %% ------------------------------------------------------------------
0143         %Evaluate Forward
0144         if maximum_lag > 0 && n_pred > 0
0145             indx_r = find(M_.block_structure.block(i).lead_lag_incidence(1,:));
0146             indx_c = M_.block_structure.block(i).lead_lag_incidence(1,indx_r);
0147             data(i).eigval = diag(jacob(indx_r, indx_c));
0148             data(i).rank = 0;
0149         else
0150             data(i).eigval = [];
0151             data(i).rank = 0;
0152         end
0153         dr.eigval = [dr.eigval ; data(i).eigval];
0154         dr.rank = dr.rank + data(i).rank;
0155         %First order approximation
0156         if task ~= 1
0157             [tmp1, tmp2, indx_c] = find(M_.block_structure.block(i).lead_lag_incidence(2,:));
0158             B = jacob(:,indx_c);
0159             if (maximum_lag > 0 && n_pred > 0)
0160                 [indx_r, tmp1, indx_r_v]  = find(M_.block_structure.block(i).lead_lag_incidence(1,:));
0161                 ghx = - B \ jacob(:,indx_r_v);
0162             end;
0163             if other_endo_nbr
0164                 fx = data(i).g1_o;
0165                 % retrieves the derivatives with respect to endogenous
0166                 % variable belonging to previous blocks
0167                 fx_tm1 = zeros(n,other_endo_nbr);
0168                 fx_t = zeros(n,other_endo_nbr);
0169                 fx_tp1 = zeros(n,other_endo_nbr);
0170                 % stores in fx_tm1 the lagged values of fx
0171                 [r, c, lag] = find(data(i).lead_lag_incidence_other(1,:));
0172                 fx_tm1(:,c) = fx(:,lag);
0173                 % stores in fx the current values of fx
0174                 [r, c, curr] = find(data(i).lead_lag_incidence_other(2,:));
0175                 fx_t(:,c) = fx(:,curr);
0176                 % stores in fx_tp1 the leaded values of fx
0177                 [r, c, lead] = find(data(i).lead_lag_incidence_other(3,:));
0178                 fx_tp1(:,c) = fx(:,lead);
0179 
0180                 l_x = dr.ghx(data(i).other_endogenous,:);
0181                 l_x_sv = dr.ghx(dr.state_var, 1:n_sv);
0182 
0183                 selector_tm1 = M_.block_structure.block(i).tm1;
0184                
0185                 ghx_other = - B \ (fx_t * l_x + (fx_tp1 * l_x * l_x_sv) + fx_tm1 * selector_tm1);
0186                 dr.ghx(endo, :) = dr.ghx(endo, :) + ghx_other;
0187             end;
0188             
0189             if exo_nbr
0190                 fu = data(i).g1_x;
0191                 exo = dr.exo_var;
0192                 if other_endo_nbr > 0
0193                     l_u_sv = dr.ghu(dr.state_var,:);
0194                     l_x = dr.ghx(data(i).other_endogenous,:);
0195                     l_u = dr.ghu(data(i).other_endogenous,:);
0196                     fu_complet = zeros(n, M_.exo_nbr);
0197                     fu_complet(:,data(i).exogenous) = fu;
0198                     ghu = - B \ (fu_complet + fx_tp1 * l_x * l_u_sv + (fx_t) * l_u );
0199                 else
0200                     fu_complet = zeros(n, M_.exo_nbr);
0201                     fu_complet(:,data(i).exogenous) = fu;
0202                     ghu = - B \ fu_complet;
0203                 end;
0204             else
0205                 exo = dr.exo_var;
0206                 if other_endo_nbr > 0
0207                      l_u_sv = dr.ghu(dr.state_var,:);
0208                      l_x = dr.ghx(data(i).other_endogenous,:);
0209                      l_u = dr.ghu(data(i).other_endogenous,:);
0210                      ghu = -B \ (fx_tp1 * l_x * l_u_sv + (fx_t) * l_u );
0211                 else
0212                     ghu = [];
0213                 end
0214             end
0215         end
0216       case 2
0217       %% ------------------------------------------------------------------
0218         %Evaluate Backward
0219         if maximum_lead > 0 && n_fwrd > 0
0220             indx_r = find(M_.block_structure.block(i).lead_lag_incidence(3,:));
0221             indx_c = M_.block_structure.block(i).lead_lag_incidence(3,indx_r);
0222             data(i).eigval = 1 ./ diag(jacob(indx_r, indx_c));
0223             data(i).rank = sum(abs(data(i).eigval) > 0);
0224         else
0225             data(i).eigval = [];
0226             data(i).rank = 0;
0227         end
0228         dr.eigval = [dr.eigval ; data(i).eigval];
0229         dr.rank = dr.rank + data(i).rank;
0230         %First order approximation
0231         if task ~= 1
0232             if (maximum_lag > 0)
0233                 indx_r = find(M_.block_structure.block(i).lead_lag_incidence(3,:));
0234                 indx_c = M_.block_structure.block(i).lead_lag_incidence(3,indx_r);
0235                 ghx = - inv(jacob(indx_r, indx_c));
0236             end;
0237             ghu =  - inv(jacob(indx_r, indx_c)) * data(i).g1_x;
0238         end
0239       case 3
0240       %% ------------------------------------------------------------------
0241         %Solve Forward single equation
0242         if maximum_lag > 0 && n_pred > 0
0243             data(i).eigval = - jacob(1 , 1 : n_pred) / jacob(1 , n_pred + n_static + 1 : n_pred + n_static + n_pred + n_both);
0244             data(i).rank = 0;
0245         else
0246             data(i).eigval = [];
0247             data(i).rank = 0;
0248         end;
0249         dr.eigval = [dr.eigval ; data(i).eigval];
0250         %First order approximation
0251         if task ~= 1
0252             if (maximum_lag > 0)
0253                  ghx = - jacob(1 , 1 : n_pred) / jacob(1 , n_pred + n_static + 1 : n_pred + n_static + n_pred + n_both);
0254             else
0255                  ghx = 0;
0256             end;
0257             if other_endo_nbr
0258                 fx = data(i).g1_o;
0259                 % retrieves the derivatives with respect to endogenous
0260                 % variable belonging to previous blocks
0261                 fx_tm1 = zeros(n,other_endo_nbr);
0262                 fx_t = zeros(n,other_endo_nbr);
0263                 fx_tp1 = zeros(n,other_endo_nbr);
0264                 % stores in fx_tm1 the lagged values of fx
0265                 [r, c, lag] = find(data(i).lead_lag_incidence_other(1,:));
0266                 fx_tm1(:,c) = fx(:,lag);
0267                 % stores in fx the current values of fx
0268                 [r, c, curr] = find(data(i).lead_lag_incidence_other(2,:));
0269                 fx_t(:,c) = fx(:,curr);
0270                 % stores in fx_tm1 the leaded values of fx
0271                 [r, c, lead] = find(data(i).lead_lag_incidence_other(3,:));
0272                 fx_tp1(:,c) = fx(:,lead);
0273 
0274                 l_x = dr.ghx(data(i).other_endogenous,:);
0275                 l_x_sv = dr.ghx(dr.state_var, 1:n_sv);
0276                 
0277                 selector_tm1 = M_.block_structure.block(i).tm1;
0278                 ghx_other = - (fx_t * l_x + (fx_tp1 * l_x * l_x_sv) + fx_tm1 * selector_tm1) / jacob(1 , n_pred + 1 : n_pred + n_static + n_pred + n_both);
0279                 dr.ghx(endo, :) = dr.ghx(endo, :) + ghx_other;
0280 
0281             end;
0282             if exo_nbr
0283                 fu = data(i).g1_x;
0284                 if other_endo_nbr > 0
0285                     l_u_sv = dr.ghu(dr.state_var,:);
0286                     l_x = dr.ghx(data(i).other_endogenous,:);
0287                     l_u = dr.ghu(data(i).other_endogenous,:);
0288                     fu_complet = zeros(n, M_.exo_nbr);
0289                     fu_complet(:,data(i).exogenous) = fu;
0290                     ghu = -(fu_complet + fx_tp1 * l_x * l_u_sv + (fx_t) * l_u ) / jacob(1 , n_pred + 1 : n_pred + n_static + n_pred + n_both);
0291                     exo = dr.exo_var;
0292                 else
0293                     ghu = - fu  / jacob(1 , n_pred + 1 : n_pred + n_static + n_pred + n_both);
0294                 end;
0295             else
0296                  if other_endo_nbr > 0
0297                      l_u_sv = dr.ghu(dr.state_var,:);
0298                      l_x = dr.ghx(data(i).other_endogenous,:);
0299                      l_u = dr.ghu(data(i).other_endogenous,:);
0300                      ghu = -(fx_tp1 * l_x * l_u_sv + (fx_t) * l_u ) / jacob(1 , n_pred + 1 : n_pred + n_static + n_pred + n_both);
0301                      exo = dr.exo_var;
0302                  else
0303                      ghu = [];
0304                  end
0305             end
0306         end
0307       case 4
0308       %% ------------------------------------------------------------------
0309         %Solve Backward single equation
0310         if maximum_lead > 0 && n_fwrd > 0
0311             data(i).eigval = - jacob(1 , n_pred + n - n_fwrd + 1 : n_pred + n) / jacob(1 , n_pred + n + 1 : n_pred + n + n_fwrd) ;
0312             data(i).rank = sum(abs(data(i).eigval) > 0);
0313         else
0314             data(i).eigval = [];
0315             data(i).rank = 0;
0316         end;
0317         dr.rank = dr.rank + data(i).rank;
0318         dr.eigval = [dr.eigval ; data(i).eigval];
0319       case 6
0320       %% ------------------------------------------------------------------
0321         %Solve Forward complete
0322         if maximum_lag > 0 && n_pred > 0
0323             data(i).eigval = eig(- jacob(: , 1 : n_pred) / ...
0324                                  jacob(: , (n_pred + n_static + 1 : n_pred + n_static + n_pred )));
0325             data(i).rank = 0;
0326         else
0327             data(i).eigval = [];
0328             data(i).rank = 0;
0329         end;
0330         dr.eigval = [dr.eigval ; data(i).eigval];
0331         if task ~= 1
0332             if (maximum_lag > 0)
0333                  ghx = - jacob(: , 1 : n_pred) / jacob(: , n_pred + n_static + 1 : n_pred + n_static + n_pred + n_both);
0334             else
0335                  ghx = 0;
0336             end;
0337             if other_endo_nbr
0338                 fx = data(i).g1_o;
0339                 % retrieves the derivatives with respect to endogenous
0340                 % variable belonging to previous blocks
0341                 fx_tm1 = zeros(n,other_endo_nbr);
0342                 fx_t = zeros(n,other_endo_nbr);
0343                 fx_tp1 = zeros(n,other_endo_nbr);
0344                 % stores in fx_tm1 the lagged values of fx
0345                 [r, c, lag] = find(data(i).lead_lag_incidence_other(1,:));
0346                 fx_tm1(:,c) = fx(:,lag);
0347                 % stores in fx the current values of fx
0348                 [r, c, curr] = find(data(i).lead_lag_incidence_other(2,:));
0349                 fx_t(:,c) = fx(:,curr);
0350                 % stores in fx_tm1 the leaded values of fx
0351                 [r, c, lead] = find(data(i).lead_lag_incidence_other(3,:));
0352                 fx_tp1(:,c) = fx(:,lead);
0353 
0354                 l_x = dr.ghx(data(i).other_endogenous,:);
0355                 l_x_sv = dr.ghx(dr.state_var, 1:n_sv);
0356                 
0357                 selector_tm1 = M_.block_structure.block(i).tm1;
0358                 ghx_other = - (fx_t * l_x + (fx_tp1 * l_x * l_x_sv) + fx_tm1 * selector_tm1) / jacob(: , n_pred + 1 : n_pred + n_static + n_pred + n_both);
0359                 dr.ghx(endo, :) = dr.ghx(endo, :) + ghx_other;
0360             end;
0361             if exo_nbr
0362                 fu = data(i).g1_x;
0363                 if other_endo_nbr > 0
0364                     l_u_sv = dr.ghu(dr.state_var,:);
0365                     l_x = dr.ghx(data(i).other_endogenous,:);
0366                     l_u = dr.ghu(data(i).other_endogenous,:);
0367                     fu_complet = zeros(n, M_.exo_nbr);
0368                     fu_complet(:,data(i).exogenous) = fu;
0369                     ghu = -(fu_complet + fx_tp1 * l_x * l_u_sv + (fx_t) * l_u ) / jacob(: , n_pred + 1 : n_pred + n_static + n_pred + n_both);
0370                     exo = dr.exo_var;
0371                 else
0372                     ghu = - fu  / jacob(: , n_pred + 1 : n_pred + n_static + n_pred + n_both);
0373                 end;
0374             else
0375                  if other_endo_nbr > 0
0376                      l_u_sv = dr.ghu(dr.state_var,:);
0377                      l_x = dr.ghx(data(i).other_endogenous,:);
0378                      l_u = dr.ghu(data(i).other_endogenous,:);
0379                      ghu = -(fx_tp1 * l_x * l_u_sv + (fx_t) * l_u ) / jacob(1 , n_pred + 1 : n_pred + n_static + n_pred + n_both);
0380                      exo = dr.exo_var;
0381                  else
0382                      ghu = [];
0383                  end
0384             end
0385         end
0386       case 7
0387       %% ------------------------------------------------------------------
0388         %Solve Backward complete
0389         if maximum_lead > 0 && n_fwrd > 0
0390             data(i).eigval = eig(- jacob(: , n_pred + n - n_fwrd + 1: n_pred + n))/ ...
0391                 jacob(: , n_pred + n + 1 : n_pred + n + n_fwrd);
0392             data(i).rank = sum(abs(data(i).eigval) > 0);
0393         else
0394             data(i).eigval = [];
0395             data(i).rank = 0;
0396         end;
0397         dr.rank = dr.rank + data(i).rank;
0398         dr.eigval = [dr.eigval ; data(i).eigval];
0399       case {5,8}
0400       %% ------------------------------------------------------------------
0401         %The lead_lag_incidence contains columns in the following order:
0402         %  static variables, backward variable, mixed variables and forward variables
0403         %
0404         %Proceeds to a QR decomposition on the jacobian matrix in order to reduce the problem size
0405         index_c = lead_lag_incidence(2,:);             % Index of all endogenous variables present at time=t
0406         index_s = lead_lag_incidence(2,1:n_static);    % Index of all static endogenous variables present at time=t
0407         if n_static > 0
0408             [Q, junk] = qr(jacob(:,index_s));
0409             aa = Q'*jacob;
0410         else
0411             aa = jacob;
0412         end;
0413         index_0m = (n_static+1:n_static+n_pred) + indexi_0 - 1;
0414         index_0p = (n_static+n_pred+1:n) + indexi_0 - 1;
0415         index_m = 1:(n_pred+n_both);
0416         indexi_p = max(lead_lag_incidence(2,:))+1;
0417         index_p = indexi_p:size(jacob, 2);
0418         nyf = n_fwrd + n_both;
0419         A = aa(:,index_m);  % Jacobain matrix for lagged endogeneous variables
0420         B = aa(:,index_c);  % Jacobian matrix for contemporaneous endogeneous variables
0421         C = aa(:,index_p);  % Jacobain matrix for led endogeneous variables
0422 
0423         row_indx = n_static+1:n;
0424         
0425         D = [[aa(row_indx,index_0m) zeros(n_dynamic,n_both) aa(row_indx,index_p)] ; [zeros(n_both, n_pred) eye(n_both) zeros(n_both, n_both + n_fwrd)]];
0426         E = [-aa(row_indx,[index_m index_0p])  ; [zeros(n_both, n_both + n_pred) eye(n_both, n_both + n_fwrd) ] ];
0427 
0428         [err, ss, tt, w, sdim, data(i).eigval, info1] = mjdgges(E,D,options_.qz_criterium);
0429 
0430         if (verbose)
0431             disp('eigval');
0432             disp(data(i).eigval);
0433         end;
0434         if info1
0435             info(1) = 2;
0436             info(2) = info1;
0437             return
0438         end
0439         nba = nd-sdim;
0440         if task == 1
0441             data(i).rank = rank(w(nd-nyf+1:end,nd-nyf+1:end));
0442             dr.rank = dr.rank + data(i).rank;
0443             if ~exist('OCTAVE_VERSION','builtin')
0444                 data(i).eigval = eig(E,D);
0445             end
0446             dr.eigval = [dr.eigval ; data(i).eigval];
0447         end
0448         if (verbose)
0449             disp(['sum eigval > 1 = ' int2str(sum(abs(data(i).eigval) > 1.)) ' nyf=' int2str(nyf) ' and dr.rank=' int2str(data(i).rank)]);
0450             disp(['data(' int2str(i) ').eigval']);
0451             disp(data(i).eigval);
0452         end;
0453         
0454         %First order approximation
0455         if task ~= 1
0456             if nba ~= nyf
0457                 sorted_roots = sort(abs(dr.eigval));
0458                 if isfield(options_,'indeterminacy_continuity')
0459                     if options_.indeterminacy_msv == 1
0460                         [ss,tt,w,q] = qz(e',d');
0461                         [ss,tt,w,junk] = reorder(ss,tt,w,q);
0462                         ss = ss';
0463                         tt = tt';
0464                         w  = w';
0465                         %nba = nyf;
0466                     end
0467                 else
0468                     if nba > nyf
0469                         temp = sorted_roots(nd-nba+1:nd-nyf)-1-options_.qz_criterium;
0470                         info(1) = 3;
0471                     elseif nba < nyf;
0472                         temp = sorted_roots(nd-nyf+1:nd-nba)-1-options_.qz_criterium;
0473                         info(1) = 4;
0474                     end
0475                     info(2) = temp'*temp;
0476                     return
0477                 end
0478             end
0479             indx_stable_root = 1: (nd - nyf);     %=> index of stable roots
0480             indx_explosive_root = n_pred + n_both + 1:nd;  %=> index of explosive roots
0481             % derivatives with respect to dynamic state variables
0482             % forward variables
0483             Z = w';
0484             Z11t = Z(indx_stable_root,    indx_stable_root)';
0485             Z21  = Z(indx_explosive_root, indx_stable_root);
0486             Z22  = Z(indx_explosive_root, indx_explosive_root);
0487             if ~isfloat(Z21) && (condest(Z21) > 1e9)
0488                 % condest() fails on a scalar under Octave
0489                 info(1) = 5;
0490                 info(2) = condest(Z21);
0491                 return;
0492             else
0493                 %gx = -inv(Z22) * Z21;
0494                 gx = - Z22 \ Z21;
0495             end
0496 
0497             % predetermined variables
0498             hx =  Z11t * inv(tt(indx_stable_root, indx_stable_root)) * ss(indx_stable_root, indx_stable_root) * inv(Z11t);
0499             
0500             k1 = 1:(n_pred+n_both);
0501             k2 = 1:(n_fwrd+n_both);
0502 
0503             ghx = [hx(k1,:); gx(k2(n_both+1:end),:)];
0504             
0505             %lead variables actually present in the model
0506             
0507             j4 = n_static+n_pred+1:n_static+n_pred+n_both+n_fwrd;   % Index on the forward and both variables
0508             j3 = nonzeros(lead_lag_incidence(2,j4)) - n_static - 2 * n_pred - n_both;  % Index on the non-zeros forward and both variables
0509             j4 = find(lead_lag_incidence(2,j4)); 
0510             
0511             if n_static > 0
0512                 B_static = B(:,1:n_static);  % submatrix containing the derivatives w.r. to static variables
0513             else
0514                 B_static = [];
0515             end;
0516             %static variables, backward variable, mixed variables and forward variables
0517             B_pred = B(:,n_static+1:n_static+n_pred+n_both);
0518             B_fyd = B(:,n_static+n_pred+n_both+1:end);
0519 
0520             % static variables
0521             if n_static > 0
0522                 temp = - C(1:n_static,j3)*gx(j4,:)*hx;
0523                 j5 = index_m;
0524                 b = aa(:,index_c);
0525                 b10 = b(1:n_static, 1:n_static);
0526                 b11 = b(1:n_static, n_static+1:n);
0527                 temp(:,j5) = temp(:,j5)-A(1:n_static,:);
0528                 temp = b10\(temp-b11*ghx);
0529                 ghx = [temp; ghx];
0530                 temp = [];
0531             end;
0532             
0533             A_ = real([B_static C(:,j3)*gx+B_pred B_fyd]); % The state_variable of the block are located at [B_pred B_both]
0534             
0535             if other_endo_nbr
0536                 if n_static > 0
0537                      fx = Q' * data(i).g1_o;
0538                 else
0539                     fx = data(i).g1_o;
0540                 end;
0541                 % retrieves the derivatives with respect to endogenous
0542                 % variable belonging to previous blocks
0543                 fx_tm1 = zeros(n,other_endo_nbr);
0544                 fx_t = zeros(n,other_endo_nbr);
0545                 fx_tp1 = zeros(n,other_endo_nbr);
0546                 % stores in fx_tm1 the lagged values of fx
0547                 [r, c, lag] = find(data(i).lead_lag_incidence_other(1,:));
0548                 fx_tm1(:,c) = fx(:,lag);
0549                 % stores in fx the current values of fx
0550                 [r, c, curr] = find(data(i).lead_lag_incidence_other(2,:));
0551                 fx_t(:,c) = fx(:,curr);
0552                 % stores in fx_tp1 the leaded values of fx
0553                 [r, c, lead] = find(data(i).lead_lag_incidence_other(3,:));
0554                 fx_tp1(:,c) = fx(:,lead);
0555 
0556                 l_x = dr.ghx(data(i).other_endogenous,:);
0557                 
0558                 l_x_sv = dr.ghx(dr.state_var, :);
0559                 
0560                 selector_tm1 = M_.block_structure.block(i).tm1; 
0561 
0562                 B_ = [zeros(size(B_static)) zeros(n,n_pred) C(:,j3) ];
0563                 C_ = l_x_sv;
0564                 D_ = (fx_t * l_x + fx_tp1 * l_x * l_x_sv + fx_tm1 * selector_tm1 );
0565                 % Solve the Sylvester equation:
0566                 % A_ * gx + B_ * gx * C_ + D_ = 0
0567                 if block_type == 5
0568                     vghx_other = - inv(kron(eye(size(D_,2)), A_) + kron(C_', B_)) * vec(D_);
0569                     ghx_other = reshape(vghx_other, size(D_,1), size(D_,2));
0570                 elseif options_.sylvester_fp == 1
0571                     ghx_other = gensylv_fp(A_, B_, C_, D_, i, options_.sylvester_fixed_point_tol);
0572                 else 
0573                     [err, ghx_other] = gensylv(1, A_, B_, C_, -D_);
0574                 end;
0575                 if options_.aim_solver ~= 1 && options_.use_qzdiv
0576                    % Necessary when using Sims' routines for QZ
0577                    ghx_other = real(ghx_other);
0578                 end
0579                 
0580                 dr.ghx(endo, :) = dr.ghx(endo, :) + ghx_other;
0581             end;
0582 
0583             if exo_nbr
0584                 if n_static > 0
0585                     fu = Q' * data(i).g1_x;
0586                 else
0587                     fu = data(i).g1_x;
0588                 end;
0589 
0590                 if other_endo_nbr > 0
0591                     l_u_sv = dr.ghu(dr.state_var,:);
0592                     l_x = dr.ghx(data(i).other_endogenous,:);
0593                     l_u = dr.ghu(data(i).other_endogenous,:);
0594                     fu_complet = zeros(n, M_.exo_nbr);
0595                     fu_complet(:,data(i).exogenous) = fu;
0596                     % Solve the equation in ghu:
0597                     % A_ * ghu + (fu_complet + fx_tp1 * l_x * l_u_sv + (fx_t + B_ * ghx_other) * l_u ) = 0
0598                     
0599                     ghu = -A_\ (fu_complet + fx_tp1 * l_x * l_u_sv + fx_t * l_u + B_ * ghx_other  * l_u_sv  );
0600                     exo = dr.exo_var;
0601                 else
0602                     ghu = - A_ \ fu;
0603                 end;
0604             else
0605                 if other_endo_nbr > 0
0606                     l_u_sv = dr.ghu(dr.state_var,:);
0607                     l_x = dr.ghx(data(i).other_endogenous,:);
0608                     l_u = dr.ghu(data(i).other_endogenous,:);
0609                     % Solve the equation in ghu:
0610                     % A_ * ghu + (fx_tp1 * l_x * l_u_sv + (fx_t + B_ * ghx_other) * l_u ) = 0
0611                     ghu = -real(A_)\ (fx_tp1 * l_x * l_u_sv + (fx_t * l_u + B_ * ghx_other * l_u_sv) );
0612                     exo = dr.exo_var;
0613                 else
0614                     ghu = [];
0615                 end;
0616             end
0617 
0618 
0619             
0620             if options_.loglinear == 1
0621                 error('log linear option is for the moment not supported in first order approximation for a block decomposed mode');
0622 %                 k = find(dr.kstate(:,2) <= M_.maximum_endo_lag+1);
0623 %                 klag = dr.kstate(k,[1 2]);
0624 %                 k1 = dr.order_var;
0625 %
0626 %                 ghx = repmat(1./dr.ys(k1),1,size(ghx,2)).*ghx.* ...
0627 %                       repmat(dr.ys(k1(klag(:,1)))',size(ghx,1),1);
0628 %                 ghu = repmat(1./dr.ys(k1),1,size(ghu,2)).*ghu;
0629             end
0630 
0631 
0632             if options_.aim_solver ~= 1 && options_.use_qzdiv
0633                 % Necessary when using Sims' routines for QZ
0634                 ghx = real(ghx);
0635                 ghu = real(ghu);
0636             end
0637 
0638             %exogenous deterministic variables
0639             if exo_det_nbr > 0
0640                 error('deterministic exogenous are not yet implemented in first order approximation for a block decomposed model');
0641 %                 f1 = sparse(jacobia_(:,nonzeros(M_.lead_lag_incidence(M_.maximum_endo_lag+2:end,order_var))));
0642 %                 f0 = sparse(jacobia_(:,nonzeros(M_.lead_lag_incidence(M_.maximum_endo_lag+1,order_var))));
0643 %                 fudet = data(i).g1_xd;
0644 %                 M1 = inv(f0+[zeros(n,n_static) f1*gx zeros(n,nyf-n_both)]);
0645 %                 M2 = M1*f1;
0646 %                 dr.ghud = cell(M_.exo_det_length,1);
0647 %                 dr.ghud{1} = -M1*fudet;
0648 %                 for i = 2:M_.exo_det_length
0649 %                     dr.ghud{i} = -M2*dr.ghud{i-1}(end-nyf+1:end,:);
0650 %                 end
0651             end
0652         end
0653     end;
0654     if task ~=1
0655         if (maximum_lag > 0 && (n_pred > 0 || n_both > 0))
0656             sorted_col_dr_ghx = M_.block_structure.block(i).sorted_col_dr_ghx;
0657             dr.ghx(endo, sorted_col_dr_ghx) = dr.ghx(endo, sorted_col_dr_ghx) + ghx;
0658             data(i).ghx = ghx;
0659             data(i).pol.i_ghx = sorted_col_dr_ghx;
0660         else
0661             data(i).pol.i_ghx = [];
0662         end;
0663         data(i).ghu = ghu;
0664         dr.ghu(endo, exo) = ghu;
0665         data(i).pol.i_ghu = exo;
0666     end;
0667     
0668    if (verbose)
0669         disp('dr.ghx');
0670         dr.ghx
0671         disp('dr.ghu');
0672         dr.ghu
0673    end; 
0674    
0675 end;
0676 M_.block_structure.block = data ;
0677 if (verbose)
0678         disp('dr.ghx');
0679         disp(real(dr.ghx));
0680         disp('dr.ghu');
0681         disp(real(dr.ghu));
0682 end; 
0683 if (task == 1)
0684     return;
0685 end;