sim1a

 function sim1
 performs deterministic simulations with lead or lag on one period

 INPUTS
   ...
 OUTPUTS
   ...
 ALGORITHM
   Laffargue, Boucekkine, Juillard (LBJ)
   see Juillard (1996) Dynare: A program for the resolution and
   simulation of dynamic models with forward variables through the use
   of a relaxation algorithm. CEPREMAP. Couverture Orange. 9602.

 SPECIAL REQUIREMENTS
   None.

0001 function sim1
0002 % function sim1
0003 % performs deterministic simulations with lead or lag on one period
0004 %
0005 % INPUTS
0006 %   ...
0007 % OUTPUTS
0008 %   ...
0009 % ALGORITHM
0010 %   Laffargue, Boucekkine, Juillard (LBJ)
0011 %   see Juillard (1996) Dynare: A program for the resolution and
0012 %   simulation of dynamic models with forward variables through the use
0013 %   of a relaxation algorithm. CEPREMAP. Couverture Orange. 9602.
0014 %
0015 % SPECIAL REQUIREMENTS
0016 %   None.
0017 
0018 % Copyright (C) 1996-2010 Dynare Team
0019 %
0020 % This file is part of Dynare.
0021 %
0022 % Dynare is free software: you can redistribute it and/or modify
0023 % it under the terms of the GNU General Public License as published by
0024 % the Free Software Foundation, either version 3 of the License, or
0025 % (at your option) any later version.
0026 %
0027 % Dynare is distributed in the hope that it will be useful,
0028 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0029 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0030 % GNU General Public License for more details.
0031 %
0032 % You should have received a copy of the GNU General Public License
0033 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0034 
0035 global M_ options_ oo_
0036 
0037 lead_lag_incidence = M_.lead_lag_incidence;
0038 
0039 ny = M_.endo_nbr;
0040 
0041 max_lag = M_.maximum_endo_lag;
0042 
0043 nyp = nnz(lead_lag_incidence(1,:)) ;
0044 iyp = find(lead_lag_incidence(1,:)>0) ;
0045 ny0 = nnz(lead_lag_incidence(2,:)) ;
0046 iy0 = find(lead_lag_incidence(2,:)>0) ;
0047 nyf = nnz(lead_lag_incidence(3,:)) ;
0048 iyf = find(lead_lag_incidence(3,:)>0) ;
0049 
0050 nd = nyp+ny0+nyf;
0051 nrc = nyf+1 ;
0052 isp = [1:nyp] ;
0053 is = [nyp+1:ny+nyp] ;
0054 isf = iyf+nyp ;
0055 isf1 = [nyp+ny+1:nyf+nyp+ny+1] ;
0056 stop = 0 ;
0057 iz = [1:ny+nyp+nyf];
0058 
0059 periods = options_.periods
0060 steady_state = oo_.steady_state;
0061 params = M_.params;
0062 endo_simul = oo_.endo_simul;
0063 exo_simul = oo_.exo_simul;
0064 i_cols_1 = nonzeros(lead_lag_incidence(2:3,:)');
0065 i_cols_A1 = find(lead_lag_incidence(2:3,:)');
0066 i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
0067 i_cols_j = 1:nd;
0068 i_upd = ny+(1:periods*ny);
0069 
0070 Y = endo_simul(:);
0071 
0072 disp (['-----------------------------------------------------']) ;
0073 disp (['MODEL SIMULATION :']) ;
0074 fprintf('\n') ;
0075 
0076 
0077 model_dynamic = str2func([M_.fname,'_dynamic']);
0078 z = Y(find(lead_lag_incidence'));
0079 [d1,jacobian] = model_dynamic(z,oo_.exo_simul, params, ...
0080                               steady_state,2);
0081 
0082 A = sparse([],[],[],periods*ny,periods*ny,periods*nnz(jacobian));
0083 res = zeros(periods*ny,1);
0084 
0085     
0086 h1 = clock ;
0087 for iter = 1:options_.maxit_
0088     h2 = clock ;
0089     
0090     i_rows = 1:ny;
0091     i_cols = find(lead_lag_incidence');
0092     i_cols_A = i_cols;
0093     
0094     for it = 2:(periods+1)
0095 
0096         [d1,jacobian] = model_dynamic(Y(i_cols),exo_simul, params, ...
0097                                       steady_state,it);
0098         if it == 2
0099             A(i_rows,i_cols_A1) = jacobian(:,i_cols_1);
0100         elseif it == periods+1
0101             A(i_rows,i_cols_A(i_cols_T)) = jacobian(:,i_cols_T);
0102         else
0103             A(i_rows,i_cols_A) = jacobian(:,i_cols_j);
0104         end
0105 
0106         res(i_rows) = d1;
0107         
0108         i_rows = i_rows + ny;
0109         i_cols = i_cols + ny;
0110         if it > 2
0111             i_cols_A = i_cols_A + ny;
0112         end
0113     end
0114         
0115     err = max(abs(res));
0116     
0117     if err < options_.dynatol.f
0118         stop = 1 ;
0119         fprintf('\n') ;
0120         disp([' Total time of simulation        :' num2str(etime(clock,h1))]) ;
0121         fprintf('\n') ;
0122         disp([' Convergency obtained.']) ;
0123         fprintf('\n') ;
0124         oo_.deterministic_simulation.status = 1;% Convergency obtained.
0125         oo_.deterministic_simulation.error = err;
0126         oo_.deterministic_simulation.iterations = iter;
0127         oo_.endo_simul = reshape(Y,ny,periods+2);
0128         break
0129     end
0130 
0131     dy = -A\res;
0132     
0133     Y(i_upd) =   Y(i_upd) + dy;
0134 
0135 end
0136 
0137 
0138 if ~stop
0139     fprintf('\n') ;
0140     disp(['     Total time of simulation        :' num2str(etime(clock,h1))]) ;
0141     fprintf('\n') ;
0142     disp(['WARNING : maximum number of iterations is reached (modify options_.maxit_).']) ;
0143     fprintf('\n') ;
0144     oo_.deterministic_simulation.status = 0;% more iterations are needed.
0145     oo_.deterministic_simulation.error = err;
0146     oo_.deterministic_simulation.errors = c/abs(err);    
0147     oo_.deterministic_simulation.iterations = options_.maxit_;
0148 end
0149 disp (['-----------------------------------------------------']) ;
0150