perfect_foresight_simulation

 Performs deterministic simulations with lead or lag on one period

 INPUTS
   endo_simul                  [double]     n*T matrix, where n is the number of endogenous variables.
   exo_simul                   [double]     q*T matrix, where q is the number of shocks.
   compute_linear_solution     [integer]    scalar equal to zero or one.

 OUTPUTS
   none

 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 info = perfect_foresight_simulation(compute_linear_solution,steady_state)
0002 % Performs deterministic simulations with lead or lag on one period
0003 %
0004 % INPUTS
0005 %   endo_simul                  [double]     n*T matrix, where n is the number of endogenous variables.
0006 %   exo_simul                   [double]     q*T matrix, where q is the number of shocks.
0007 %   compute_linear_solution     [integer]    scalar equal to zero or one.
0008 %
0009 % OUTPUTS
0010 %   none
0011 %
0012 % ALGORITHM
0013 %   Laffargue, Boucekkine, Juillard (LBJ)
0014 %   see Juillard (1996) Dynare: A program for the resolution and
0015 %   simulation of dynamic models with forward variables through the use
0016 %   of a relaxation algorithm. CEPREMAP. Couverture Orange. 9602.
0017 %
0018 % SPECIAL REQUIREMENTS
0019 %   None.
0020 
0021 % Copyright (C) 2009-2010 Dynare Team
0022 %
0023 % This file is part of Dynare.
0024 %
0025 % Dynare is free software: you can redistribute it and/or modify
0026 % it under the terms of the GNU General Public License as published by
0027 % the Free Software Foundation, either version 3 of the License, or
0028 % (at your option) any later version.
0029 %
0030 % Dynare is distributed in the hope that it will be useful,
0031 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0032 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0033 % GNU General Public License for more details.
0034 %
0035 % You should have received a copy of the GNU General Public License
0036 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0037 
0038 global M_ options_ it_ oo_
0039 
0040 persistent lead_lag_incidence dynamic_model ny nyp nyf nrs nrc iyf iyp isp is isf isf1 iz icf ghx iflag
0041 
0042 if ~nargin && isempty(iflag)% Initialization of the persistent variables.
0043     lead_lag_incidence = M_.lead_lag_incidence;
0044     dynamic_model = [M_.fname '_dynamic'];
0045     ny   = size(oo_.endo_simul,1);
0046     nyp  = nnz(lead_lag_incidence(1,:));% number of lagged variables.
0047     nyf  = nnz(lead_lag_incidence(3,:));% number of leaded variables.
0048     nrs  = ny+nyp+nyf+1;
0049     nrc  = nyf+1;
0050     iyf  = find(lead_lag_incidence(3,:)>0);% indices for leaded variables.
0051     iyp  = find(lead_lag_incidence(1,:)>0);% indices for lagged variables.
0052     isp  = 1:nyp;
0053     is   = (nyp+1):(nyp+ny); % Indices for contemporaneaous variables.
0054     isf  = iyf+nyp;
0055     isf1 = (nyp+ny+1):(nyf+nyp+ny+1);
0056     iz   = 1:(ny+nyp+nyf);
0057     icf  = 1:size(iyf,2);
0058     info = [];
0059     iflag = 1;
0060     return
0061 end
0062 
0063 initial_endo_simul = oo_.endo_simul;
0064 
0065 warning_old_state = warning;
0066 warning off all
0067 
0068 if nargin<1
0069     compute_linear_solution = 0;
0070 else
0071     if nargin<2
0072         error('The steady state (second input argument) is missing!');
0073     end
0074 end
0075 
0076 if ~isstruct(compute_linear_solution) && compute_linear_solution
0077     [dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
0078 elseif isstruct(compute_linear_solution)
0079     dr = compute_linear_solution;
0080     compute_linear_solution = 1;
0081 end
0082 
0083 if compute_linear_solution
0084     ghx(dr.order_var,:) = dr.ghx;
0085     ghx = ghx(iyf,:);
0086 end
0087 
0088 periods = options_.periods;
0089 
0090 stop    = 0 ;
0091 it_init = M_.maximum_lag+1;
0092 
0093 info.convergence = 1;
0094 info.time  = 0;
0095 info.error = 0;
0096 info.iterations.time  = zeros(options_.maxit_,1);
0097 info.iterations.error = info.iterations.time;
0098 
0099 last_line = options_.maxit_;
0100 error_growth = 0;
0101 
0102 h1 = clock;
0103 for iter = 1:options_.maxit_
0104     h2 = clock;
0105     if options_.terminal_condition
0106         c = zeros(ny*(periods+1),nrc);
0107     else
0108         c = zeros(ny*periods,nrc);
0109     end
0110     it_ = it_init;
0111     z = [ oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1) ];
0112     [d1,jacobian] = feval(dynamic_model,z,oo_.exo_simul, M_.params, it_);
0113     jacobian = [jacobian(:,iz) , -d1];
0114     ic = 1:ny;
0115     icp = iyp;
0116     c(ic,:) = jacobian(:,is)\jacobian(:,isf1) ;
0117     for it_ = it_init+(1:periods-1-(options_.terminal_condition==2))
0118         z = [ oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1)];
0119         [d1,jacobian] = feval(dynamic_model,z,oo_.exo_simul, M_.params, it_);
0120         jacobian = [jacobian(:,iz) , -d1];
0121         jacobian(:,[isf nrs]) = jacobian(:,[isf nrs])-jacobian(:,isp)*c(icp,:);
0122         ic = ic + ny;
0123         icp = icp + ny;
0124         c(ic,:) = jacobian(:,is)\jacobian(:,isf1);
0125     end
0126     if options_.terminal_condition
0127         if options_.terminal_condition==1% Terminal condition is Y_{T} = Y_{T+1}
0128             s = eye(ny);
0129             s(:,isf) = s(:,isf)+c(ic,1:nyf);
0130             ic = ic + ny;
0131             c(ic,nrc) = s\c(ic,nrc);
0132         else% Terminal condition is Y_{T+1}-Y^{\star} = TransitionMatrix*(Y_{T}-Y^{\star})
0133             it_ = it_+1;
0134             z = [ oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1) ] ;
0135             [d1,jacobian] = feval(dynamic_model,z,oo_.exo_simul, M_.params, it_);
0136             jacobian = [jacobian(:,iz) -d1];
0137             jacobian(:,[isf nrs]) = jacobian(:,[isf nrs])-jacobian(:,isp)*c(icp,:) ;
0138             ic = ic + ny;
0139             icp = icp + ny;
0140             s = jacobian(:,is);
0141             s(:,iyp) = s(:,iyp)+jacobian(:,isf)*ghx;
0142             c (ic,:) = s\jacobian(:,isf1);
0143         end
0144         c = bksup0(c,ny,nrc,iyf,icf,periods);
0145         c = reshape(c,ny,periods+1);
0146         oo_.endo_simul(:,it_init+(0:periods)) = oo_.endo_simul(:,it_init+(0:periods))+options_.slowc*c;
0147     else% Terminal condition is Y_{T}=Y^{\star}
0148         c = bksup0(c,ny,nrc,iyf,icf,periods);
0149         c = reshape(c,ny,periods);
0150         oo_.endo_simul(:,it_init+(0:periods-1)) = oo_.endo_simul(:,it_init+(0:periods-1))+options_.slowc*c;
0151     end
0152     err = max(max(abs(c)));
0153     info.iterations.time(iter)  = etime(clock,h2);
0154     info.iterations.error(iter) = err;
0155     if iter>1
0156         error_growth = error_growth + (info.iterations.error(iter)>info.iterations.error(iter-1));
0157     end
0158     if isnan(err) || error_growth>3
0159         last_line = iter;
0160         break
0161     end
0162     if err < options_.dynatol.f
0163         stop = 1;
0164         info.time  = etime(clock,h1);
0165         info.error = err;
0166         info.iterations.time = info.iterations.time(1:iter);
0167         info.iterations.error  = info.iterations.error(1:iter);
0168         break
0169     end
0170 end
0171 
0172 if stop && options_.terminal_condition==2
0173     % Compute the distance to the deterministic steady state (for the subset of endogenous variables with a non zero
0174     % steady state) at the last perdiod.
0175     idx = find(abs(oo_.steady_state)>0);
0176     distance_to_steady_state = abs(((oo_.endo_simul(idx,end)-oo_.steady_state(idx))./oo_.steady_state(idx)))*100;
0177     disp(['(max) Distance to steady state at the end is (in percentage):' num2str(max(distance_to_steady_state))])
0178 end
0179 
0180 if ~stop
0181     info.time  = etime(clock,h1);
0182     info.error = err;
0183     info.convergence = 0;
0184     info.iterations.time  = info.iterations.time(1:last_line);
0185     info.iterations.error = info.iterations.error(1:last_line);
0186     oo_.endo_simul = initial_endo_simul;
0187 end
0188 
0189 warning(warning_old_state);