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