1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
|
function [y, info] = solve_one_boundary(fname, y, x, params, steady_state, ...
y_index_eq, nze, periods, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, lambda, cutoff, stack_solve_algo, forward_backward, is_dynamic, verbose, M, options, oo)
% Computes the deterministic simulation of a block of equation containing
% lead or lag variables
%
% INPUTS
% fname [string] name of the file containing the block
% to simulate
% y [matrix] All the endogenous variables of the model
% x [matrix] All the exogenous variables of the model
% params [vector] All the parameters of the model
% steady_state [vector] steady state of the model
% y_index_eq [vector of int] The index of the endogenous variables of
% the block
% nze [integer] number of non-zero elements in the
% jacobian matrix
% periods [integer] number of simulation periods
% is_linear [integer] if is_linear=1 the block is linear
% if is_linear=0 the block is not linear
% Block_Num [integer] block number
% y_kmin [integer] maximum number of lag in the model
% maxit_ [integer] maximum number of iteration in Newton
% solve_tolf [double] convergence criteria
% lambda [double] initial value of step size in
% Newton
% cutoff [double] cutoff to correct the direction in Newton in case
% of singular jacobian matrix
% stack_solve_algo [integer] linear solver method used in the
% Newton algorithm :
% - 1 sparse LU
% - 2 GMRES
% - 3 BicGStab
% - 4 Optimal path length
% forward_backward [integer] The block has to be solve forward
% (1) or backward (0)
% is_dynamic [integer] (1) The block belong to the dynamic
% file and the oo_.deterministic_simulation field has to be uptated
% (0) The block belong to the static
% file and th oo_.deteerminstic
% field remains unchanged
% verbose [integer] (0) iterations are not printed
% (1) iterations are printed
% indirect_call [integer] (0) direct call to the fname
% (1) indirect call via the
% local_fname wrapper
% OUTPUTS
% y [matrix] All endogenous variables of the model
% info [integer] >=0 no error
% <0 error
%
% ALGORITHM
% Newton with LU or GMRES or BicGstab for dynamic block
%
% SPECIAL REQUIREMENTS
% none.
%
% Copyright (C) 1996-2017 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
Blck_size=size(y_index_eq,2);
g2 = [];
g3 = [];
correcting_factor=0.01;
ilu_setup.droptol=1e-10;
max_resa=1e100;
reduced = 0;
if forward_backward
incr = 1;
start = y_kmin+1;
finish = periods+y_kmin;
else
incr = -1;
start = periods+y_kmin;
finish = y_kmin+1;
end
for it_=start:incr:finish
cvg=0;
iter=0;
g1=spalloc( Blck_size, Blck_size, nze);
while ~(cvg==1 || iter>maxit_)
if is_dynamic
[r, y, g1, g2, g3] = feval(fname, y, x, params, steady_state, it_, 0);
else
[r, y, g1] = feval(fname, y, x, params);
end
if ~isreal(r)
max_res=(-(max(max(abs(r))))^2)^0.5;
else
max_res=max(max(abs(r)));
end
if verbose==1
disp(['iteration : ' int2str(iter+1) ' => ' num2str(max_res) ' time = ' int2str(it_)])
if is_dynamic
disp([M.endo_names(y_index_eq,:) num2str([y(it_,y_index_eq)' r g1])])
else
disp([M.endo_names(y_index_eq,:) num2str([y(y_index_eq) r g1])])
end
end
if ~isreal(max_res) || isnan(max_res)
cvg = 0;
elseif is_linear && iter>0
cvg = 1;
else
cvg=(max_res<solve_tolf);
end
if ~cvg
if iter>0
if ~isreal(max_res) || isnan(max_res) || (max_resa<max_res && iter>1)
if isnan(max_res) || (max_resa<max_res && iter>0)
detJ=det(g1a);
if(abs(detJ)<1e-7)
max_factor=max(max(abs(g1a)));
ze_elem=sum(diag(g1a)<cutoff);
if verbose
disp([num2str(full(ze_elem),'%d') ' elements on the Jacobian diagonal are below the cutoff (' num2str(cutoff,'%f') ')'])
end
if correcting_factor<max_factor
correcting_factor=correcting_factor*4;
if verbose
disp(['The Jacobain matrix is singular, det(Jacobian)=' num2str(detJ,'%f') '.'])
disp([' trying to correct the Jacobian matrix:'])
disp([' correcting_factor=' num2str(correcting_factor,'%f') ' max(Jacobian)=' num2str(full(max_factor),'%f')])
end
dx = - r/(g1+correcting_factor*speye(Blck_size));
y(it_,y_index_eq)=ya_save+lambda*dx;
continue
else
if verbose
disp('The singularity of the jacobian matrix could not be corrected')
end
info = -Block_Num*10;
return
end
end
elseif lambda>1e-8
lambda=lambda/2;
reduced = 1;
if verbose
disp(['reducing the path length: lambda=' num2str(lambda,'%f')])
end
if is_dynamic
y(it_,y_index_eq)=ya_save-lambda*dx;
else
y(y_index_eq)=ya_save-lambda*dx;
end
continue
else
if verbose
if cutoff==0
fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit".\n',Block_Num, it_, iter);
else
fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n',Block_Num, it_, iter);
end
end
if is_dynamic
oo_.deterministic_simulation.status = 0;
oo_.deterministic_simulation.error = max_res;
oo_.deterministic_simulation.iterations = iter;
oo_.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed.
oo_.deterministic_simulation.block(Block_Num).error = max_res;
oo_.deterministic_simulation.block(Block_Num).iterations = iter;
end
info = -Block_Num*10;
return
end
else
if lambda<1
lambda=max(lambda*2, 1);
end
end
end
if is_dynamic
ya = y(it_,y_index_eq)';
else
ya = y(y_index_eq);
end
ya_save=ya;
g1a=g1;
if ~is_dynamic && options.solve_algo==0
if verbose
disp('steady: fsolve')
end
if ~isoctave
if ~user_has_matlab_license('optimization_toolbox')
error('SOLVE_ONE_BOUNDARY: you can''t use solve_algo=0 since you don''t have MATLAB''s Optimization Toolbox')
end
end
options=optimset('fsolve');
options.MaxFunEvals = 50000;
options.MaxIter = 2000;
options.TolFun=1e-8;
options.Display = 'iter';
options.Jacobian = 'on';
if ~isoctave
[yn,fval,exitval,output] = fsolve(@local_fname, y(y_index_eq), ...
options, x, params, steady_state, y, y_index_eq, fname, 0);
else
% Under Octave, use a wrapper, since fsolve() does not have a 4th arg
func = @(z) local_fname(z, x, params, steady_state, y, y_index_eq, fname, 0);
% The Octave version of fsolve does not converge when it starts from the solution
fvec = feval(func,y(y_index_eq));
if max(abs(fvec)) >= options.solve_tolf
[yn,fval,exitval,output] = fsolve(func,y(y_index_eq),options);
else
yn = y(y_index_eq);
exitval = 3;
end
end
y(y_index_eq) = yn;
if exitval > 0
info = 0;
else
info = -Block_Num*10;
end
elseif (~is_dynamic && options.solve_algo==2) || (is_dynamic && stack_solve_algo==4)
if verbose==1 && ~is_dynamic
disp('steady: LU + lnsrch1')
end
lambda=1;
stpmx = 100 ;
if is_dynamic
stpmax = stpmx*max([sqrt(ya'*ya);size(y_index_eq,2)]);
else
stpmax = stpmx*max([sqrt(ya'*ya);size(y_index_eq,2)]);
end
nn=1:size(y_index_eq,2);
g = (r'*g1)';
f = 0.5*r'*r;
p = -g1\r ;
if is_dynamic
[ya,f,r,check]=lnsrch1(y(it_,:)',f,g,p,stpmax, ...
'lnsrch1_wrapper_one_boundary',nn, ...
y_index_eq, options.solve_tolx, y_index_eq, fname, y, x, params, steady_state, it_);
dx = ya' - y(it_, :);
else
[ya,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq, options.solve_tolx, x, ...
params, steady_state,0);
dx = ya - y(y_index_eq);
end
if is_dynamic
y(it_,:) = ya';
else
y = ya';
end
elseif ~is_dynamic && options.solve_algo==3
if verbose==1
disp('steady: csolve')
end
[yn,info] = csolve(@local_fname, y(y_index_eq),@ ...
local_fname,1e-6,500, x, params, steady_state, y, y_index_eq, fname, 1);
dx = ya - yn;
y(y_index_eq) = yn;
elseif (stack_solve_algo==1 && is_dynamic) || (stack_solve_algo==0 && is_dynamic) || (~is_dynamic && (options.solve_algo==1 || options.solve_algo==6))
if verbose==1 && ~is_dynamic
disp('steady: Sparse LU ')
end
dx = g1\r;
ya = ya - lambda*dx;
if is_dynamic
y(it_,y_index_eq) = ya';
else
y(y_index_eq) = ya;
end
elseif (stack_solve_algo==2 && is_dynamic) || (options.solve_algo==7 && ~is_dynamic)
flag1=1;
if isoctave
error('SOLVE_ONE_BOUNDARY: you can''t use solve_algo=7 since GMRES is not implemented in Octave')
end
if verbose == 1 && ~is_dynamic
disp('steady: GMRES ')
end
while flag1>0
[L1, U1]=ilu(g1,ilu_setup);
[dx,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1);
if flag1>0 || reduced
if verbose
if flag1==1
disp(['Error in simul: No convergence inside GMRES after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')])
elseif(flag1==2)
disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Num,'%3d')])
elseif(flag1==3)
disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')])
end
end
ilu_setup.droptol = ilu_setup.droptol/10;
reduced = 0;
else
ya = ya + lambda*dx;
if is_dynamic
y(it_,y_index_eq) = ya';
else
y(y_index_eq) = ya';
end
end
end
elseif (stack_solve_algo==3 && is_dynamic) || (options.solve_algo==8 && ~is_dynamic)
flag1=1;
if verbose == 1 && ~is_dynamic
disp('steady: BiCGStab')
end
while flag1>0
[L1, U1]=ilu(g1,ilu_setup);
phat = ya - U1 \ (L1 \ r);
if is_dynamic
y(it_,y_index_eq) = phat;
else
y(y_index_eq) = phat;
end
if is_dynamic
[r, y, g1, g2, g3] = feval(fname, y, x, params, ...
steady_state, it_, 0);
else
[r, y, g1] = feval(fname, y, x, params);
end
if max(abs(r))>=options.solve_tolf
[dx,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1);
else
flag1 = 0;
dx = phat - ya;
end
if flag1>0 || reduced
if verbose
if(flag1==1)
disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')])
elseif(flag1==2)
disp(['Error in simul: Preconditioner is ill-conditioned, in block' num2str(Block_Num,'%3d')])
elseif(flag1==3)
disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')])
end
end
ilu_setup.droptol = ilu_setup.droptol/10;
reduced = 0;
else
ya = ya + lambda*dx;
if is_dynamic
y(it_,y_index_eq) = ya';
else
y(y_index_eq) = ya';
end
end
end
else
if verbose
disp('unknown option : ')
if is_dynamic
disp(['options_.stack_solve_algo = ' num2str(stack_solve_algo) ' not implemented'])
else
disp(['options_.solve_algo = ' num2str(options.solve_algo) ' not implemented'])
end
end
info = -Block_Num*10;
return
end
iter=iter+1;
max_resa = max_res;
end
end
if cvg==0
if verbose
if cutoff == 0
fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit\".\n',Block_Num, it_,iter);
else
fprintf('Error in simul: Convergence not achieved in block %d, at time %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n',Block_Num, it_,iter);
end
end
if(is_dynamic)
oo_.deterministic_simulation.status = 0;
oo_.deterministic_simulation.error = max_res;
oo_.deterministic_simulation.iterations = iter;
oo_.deterministic_simulation.block(Block_Num).status = 0;% Convergency failed.
oo_.deterministic_simulation.block(Block_Num).error = max_res;
oo_.deterministic_simulation.block(Block_Num).iterations = iter;
end
info = -Block_Num*10;
return
end
end
if is_dynamic
info = 1;
oo_.deterministic_simulation.status = 1;
oo_.deterministic_simulation.error = max_res;
oo_.deterministic_simulation.iterations = iter;
oo_.deterministic_simulation.block(Block_Num).status = 1;
oo_.deterministic_simulation.block(Block_Num).error = max_res;
oo_.deterministic_simulation.block(Block_Num).iterations = iter;
else
info = 0;
end
function [err, G]=local_fname(yl, x, params, steady_state, y, y_index_eq, fname, is_csolve)
y(y_index_eq) = yl;
[err, y, G] = feval(fname, y, x, params, steady_state, 0);
if(is_csolve)
G = full(G);
end
|
