1function info = perfect_foresight_simulation(compute_linear_solution,steady_state) 2% Performs deterministic simulations with lead or lag on one period 3% 4% INPUTS 5% endo_simul [double] n*T matrix, where n is the number of endogenous variables. 6% exo_simul [double] q*T matrix, where q is the number of shocks. 7% compute_linear_solution [integer] scalar equal to zero or one. 8% 9% OUTPUTS 10% none 11% 12% ALGORITHM 13% Laffargue, Boucekkine, Juillard (LBJ) 14% see Juillard (1996) Dynare: A program for the resolution and 15% simulation of dynamic models with forward variables through the use 16% of a relaxation algorithm. CEPREMAP. Couverture Orange. 9602. 17% 18% SPECIAL REQUIREMENTS 19% None. 20 21% Copyright (C) 2009-2015 Dynare Team 22% 23% This file is part of Dynare. 24% 25% Dynare is free software: you can redistribute it and/or modify 26% it under the terms of the GNU General Public License as published by 27% the Free Software Foundation, either version 3 of the License, or 28% (at your option) any later version. 29% 30% Dynare is distributed in the hope that it will be useful, 31% but WITHOUT ANY WARRANTY; without even the implied warranty of 32% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 33% GNU General Public License for more details. 34% 35% You should have received a copy of the GNU General Public License 36% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 37 38global M_ options_ it_ oo_ 39 40persistent lead_lag_incidence dynamic_model ny nyp nyf nrs nrc iyf iyp isp is isf isf1 iz icf ghx iflag 41 42if ~nargin && isempty(iflag)% Initialization of the persistent variables. 43 lead_lag_incidence = M_.lead_lag_incidence; 44 dynamic_model = [M_.fname '.dynamic']; 45 ny = size(oo_.endo_simul,1); 46 nyp = nnz(lead_lag_incidence(1,:));% number of lagged variables. 47 nyf = nnz(lead_lag_incidence(3,:));% number of leaded variables. 48 nrs = ny+nyp+nyf+1; 49 nrc = nyf+1; 50 iyf = find(lead_lag_incidence(3,:)>0);% indices for leaded variables. 51 iyp = find(lead_lag_incidence(1,:)>0);% indices for lagged variables. 52 isp = 1:nyp; 53 is = (nyp+1):(nyp+ny); % Indices for contemporaneaous variables. 54 isf = iyf+nyp; 55 isf1 = (nyp+ny+1):(nyf+nyp+ny+1); 56 iz = 1:(ny+nyp+nyf); 57 icf = 1:size(iyf,2); 58 info = []; 59 iflag = 1; 60 return 61end 62 63initial_endo_simul = oo_.endo_simul; 64 65warning_old_state = warning; 66warning off all 67 68if nargin<1 69 compute_linear_solution = 0; 70else 71 if nargin<2 72 error('The steady state (second input argument) is missing!'); 73 end 74end 75 76if ~isstruct(compute_linear_solution) && compute_linear_solution 77 [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_); 78elseif isstruct(compute_linear_solution) 79 dr = compute_linear_solution; 80 compute_linear_solution = 1; 81end 82 83if compute_linear_solution 84 ghx(dr.order_var,:) = dr.ghx; 85 ghx = ghx(iyf,:); 86end 87 88periods = options_.periods; 89 90stop = 0 ; 91it_init = M_.maximum_lag+1; 92 93info.convergence = 1; 94info.time = 0; 95info.error = 0; 96info.iterations.time = zeros(options_.simul.maxit,1); 97info.iterations.error = info.iterations.time; 98 99last_line = options_.simul.maxit; 100error_growth = 0; 101 102h1 = clock; 103for iter = 1:options_.simul.maxit 104 h2 = clock; 105 if options_.terminal_condition 106 c = zeros(ny*(periods+1),nrc); 107 else 108 c = zeros(ny*periods,nrc); 109 end 110 it_ = it_init; 111 z = [ oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1) ]; 112 [d1,jacobian] = feval(dynamic_model,z,oo_.exo_simul, M_.params, it_); 113 jacobian = [jacobian(:,iz) , -d1]; 114 ic = 1:ny; 115 icp = iyp; 116 c(ic,:) = jacobian(:,is)\jacobian(:,isf1) ; 117 for it_ = it_init+(1:periods-1-(options_.terminal_condition==2)) 118 z = [ oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1)]; 119 [d1,jacobian] = feval(dynamic_model,z,oo_.exo_simul, M_.params, it_); 120 jacobian = [jacobian(:,iz) , -d1]; 121 jacobian(:,[isf nrs]) = jacobian(:,[isf nrs])-jacobian(:,isp)*c(icp,:); 122 ic = ic + ny; 123 icp = icp + ny; 124 c(ic,:) = jacobian(:,is)\jacobian(:,isf1); 125 end 126 if options_.terminal_condition 127 if options_.terminal_condition==1% Terminal condition is Y_{T} = Y_{T+1} 128 s = eye(ny); 129 s(:,isf) = s(:,isf)+c(ic,1:nyf); 130 ic = ic + ny; 131 c(ic,nrc) = s\c(ic,nrc); 132 else% Terminal condition is Y_{T+1}-Y^{\star} = TransitionMatrix*(Y_{T}-Y^{\star}) 133 it_ = it_+1; 134 z = [ oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1) ] ; 135 [d1,jacobian] = feval(dynamic_model,z,oo_.exo_simul, M_.params, it_); 136 jacobian = [jacobian(:,iz) -d1]; 137 jacobian(:,[isf nrs]) = jacobian(:,[isf nrs])-jacobian(:,isp)*c(icp,:) ; 138 ic = ic + ny; 139 icp = icp + ny; 140 s = jacobian(:,is); 141 s(:,iyp) = s(:,iyp)+jacobian(:,isf)*ghx; 142 c (ic,:) = s\jacobian(:,isf1); 143 end 144 c = bksup0(c,ny,nrc,iyf,icf,periods); 145 c = reshape(c,ny,periods+1); 146 oo_.endo_simul(:,it_init+(0:periods)) = oo_.endo_simul(:,it_init+(0:periods))+options_.slowc*c; 147 else% Terminal condition is Y_{T}=Y^{\star} 148 c = bksup0(c,ny,nrc,iyf,icf,periods); 149 c = reshape(c,ny,periods); 150 oo_.endo_simul(:,it_init+(0:periods-1)) = oo_.endo_simul(:,it_init+(0:periods-1))+options_.slowc*c; 151 end 152 err = max(max(abs(c))); 153 info.iterations.time(iter) = etime(clock,h2); 154 info.iterations.error(iter) = err; 155 if iter>1 156 error_growth = error_growth + (info.iterations.error(iter)>info.iterations.error(iter-1)); 157 end 158 if isnan(err) || error_growth>3 159 last_line = iter; 160 break 161 end 162 if err < options_.dynatol.f 163 stop = 1; 164 info.time = etime(clock,h1); 165 info.error = err; 166 info.iterations.time = info.iterations.time(1:iter); 167 info.iterations.error = info.iterations.error(1:iter); 168 break 169 end 170end 171 172if stop && options_.terminal_condition==2 173 % Compute the distance to the deterministic steady state (for the subset of endogenous variables with a non zero 174 % steady state) at the last perdiod. 175 idx = find(abs(oo_.steady_state)>0); 176 distance_to_steady_state = abs(((oo_.endo_simul(idx,end)-oo_.steady_state(idx))./oo_.steady_state(idx)))*100; 177 disp(['(max) Distance to steady state at the end is (in percentage):' num2str(max(distance_to_steady_state))]) 178end 179 180if ~stop 181 info.time = etime(clock,h1); 182 info.error = err; 183 info.convergence = 0; 184 info.iterations.time = info.iterations.time(1:last_line); 185 info.iterations.error = info.iterations.error(1:last_line); 186 oo_.endo_simul = initial_endo_simul; 187end 188 189warning(warning_old_state); 190