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
|
function [y, oo]= solve_two_boundaries(fname, y, x, params, steady_state, y_index, nze, periods, y_kmin_l, y_kmax_l, is_linear, Block_Num, y_kmin, maxit_, solve_tolf, lambda, cutoff, stack_solve_algo,options,M, oo)
% Computes the deterministic simulation of a block of equation containing
% both lead and lag variables using relaxation methods
%
% 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 [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
% y_kmin_l [integer] maximum number of lag in the block
% y_kmax_l [integer] maximum number of lead in the block
% 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 sprse LU
% - 2 GMRES
% - 3 BicGStab
% - 4 Optimal path length
% M [structure] Model description
% oo [structure] Results
%
% OUTPUTS
% y [matrix] All endogenous variables of the model
% oo [structure] Results
%
% ALGORITHM
% Newton with LU or GMRES or BicGstab
%
% SPECIAL REQUIREMENTS
% none.
%
% Copyright (C) 1996-2018 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/>.
verbose = options.verbosity;
cvg=0;
iter=0;
Per_u_=0;
g2 = [];
g3 = [];
Blck_size=size(y_index,2);
correcting_factor=0.01;
ilu_setup.droptol=1e-10;
ilu_setup.type = 'ilutp';
%ilu_setup.milu = 'col';
ilu_setup.milu = 'off';
ilu_setup.thresh = 1;
ilu_setup.udiag = 0;
max_resa=1e100;
Jacobian_Size=Blck_size*(y_kmin+y_kmax_l +periods);
g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods);
reduced = 0;
while ~(cvg==1 || iter>maxit_)
[r, y, g1, g2, g3, b]=feval(fname, y, x, params, steady_state, periods, 0, y_kmin, Blck_size,options.periods);
preconditioner = 2;
g1a=g1(:, y_kmin*Blck_size+1:(periods+y_kmin)*Blck_size);
term1 = g1(:, 1:y_kmin_l*Blck_size)*reshape(y(1+y_kmin-y_kmin_l:y_kmin,y_index)',1,y_kmin_l*Blck_size)';
term2 = g1(:, (periods+y_kmin_l)*Blck_size+1:(periods+y_kmin_l+y_kmax_l)*Blck_size)*reshape(y(periods+y_kmin+1:periods+y_kmin+y_kmax_l,y_index)',1,y_kmax_l*Blck_size)';
b = b - term1 - term2;
[max_res, max_indx]=max(max(abs(r')));
if ~isreal(r)
max_res = (-max_res^2)^0.5;
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 verbose && ~isreal(max_res)
disp(['Variable ' M.endo_names{max_indx} ' (' int2str(max_indx) ') returns an undefined value']);
end
if isnan(max_res)
detJ=det(g1aa);
if abs(detJ)<1e-7
max_factor=max(max(abs(g1aa)));
ze_elem=sum(diag(g1aa)<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 = (g1aa+correcting_factor*speye(periods*Blck_size))\ba- ya;
y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)';
continue
else
disp('The singularity of the jacobian matrix could not be corrected');
return
end
end
elseif lambda>1e-8
lambda=lambda/2;
reduced = 1;
if verbose
disp(['reducing the path length: lambda=' num2str(lambda,'%f')]);
end
y(1+y_kmin:periods+y_kmin,y_index)=reshape((ya_save+lambda*dx)',length(y_index),periods)';
continue
else
if verbose
if cutoff==0
fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.simul.maxit".\n',Block_Num, iter);
else
fprintf('Error in simul: Convergence not achieved in block %d, after %d iterations.\n Increase "options_.simul.maxit" or set "cutoff=0" in model options.\n',Block_Num, iter);
end
end
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;
return
end
else
if lambda<1
lambda=max(lambda*2, 1);
end
end
end
ya = reshape(y(y_kmin+1:y_kmin+periods,y_index)',1,periods*Blck_size)';
ya_save=ya;
g1aa=g1a;
ba=b;
max_resa=max_res;
if stack_solve_algo==0
dx = g1a\b- ya;
ya = ya + lambda*dx;
y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';
elseif stack_solve_algo==1
for t=1:periods
first_elem = (t-1)*Blck_size+1;
last_elem = t*Blck_size;
next_elem = (t+1)*Blck_size;
Elem = first_elem:last_elem;
Elem_1 = last_elem+1:next_elem;
B1_inv = inv(g1a(Elem, Elem));
if (t < periods)
S1 = B1_inv * g1a(Elem, Elem_1);
end
g1a(Elem, Elem_1) = S1;
b(Elem) = B1_inv * b(Elem);
g1a(Elem, Elem) = ones(Blck_size, Blck_size);
if t<periods
g1a(Elem_1, Elem_1) = g1a(Elem_1, Elem_1) - g1a(Elem_1, Elem) * S1;
b(Elem_1) = b(Elem_1) - g1a(Elem_1, Elem) * b(Elem);
g1a(Elem_1, Elem) = zeros(Blck_size, Blck_size);
end
end
za = b(Elem);
zaa = za;
y_Elem = (periods - 1) * Blck_size + 1:(periods) * Blck_size;
dx = ya;
dx(y_Elem) = za - ya(y_Elem);
ya(y_Elem) = ya(y_Elem) + lambda*dx(y_Elem);
for t=periods-1:-1:1
first_elem = (t-1)*Blck_size+1;
last_elem = t*Blck_size;
next_elem = (t+1)*Blck_size;
Elem_1 = last_elem+1:next_elem;
Elem = first_elem:last_elem;
za = b(Elem) - g1a(Elem, Elem_1) * zaa;
zaa = za;
y_Elem = Blck_size * (t-1)+1:Blck_size * (t);
dx(y_Elem) = za - ya(y_Elem);
ya(y_Elem) = ya(y_Elem) + lambda*dx(y_Elem);
y(y_kmin + t, y_index) = ya(y_Elem);
end
elseif stack_solve_algo==2
flag1=1;
while flag1>0
if preconditioner==2
[L1, U1]=ilu(g1a,ilu_setup);
elseif preconditioner==3
Size = Blck_size;
gss1 = g1a(Size + 1: 2*Size,Size + 1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
[L1, U1]=lu(gss1);
L(1:Size,1:Size) = L1;
U(1:Size,1:Size) = U1;
gss2 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
[L2, U2]=lu(gss2);
L(Size+1:(periods-1)*Size,Size+1:(periods-1)*Size) = kron(eye(periods-2), L2);
U(Size+1:(periods-1)*Size,Size+1:(periods-1)*Size) = kron(eye(periods-2), U2);
gss2 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size);
[L3, U3]=lu(gss2);
L((periods-1)*Size+1:periods*Size,(periods-1)*Size+1:periods*Size) = L3;
U((periods-1)*Size+1:periods*Size,(periods-1)*Size+1:periods*Size) = U3;
L1 = L;
U1 = U;
elseif preconditioner==4
Size = Blck_size;
gss1 = g1a(1: 3*Size, 1: 3*Size);
[L, U] = lu(gss1);
L1 = kron(eye(ceil(periods/3)),L);
U1 = kron(eye(ceil(periods/3)),U);
L1 = L1(1:periods * Size, 1:periods * Size);
U1 = U1(1:periods * Size, 1:periods * Size);
end
[za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);
if (flag1>0 || reduced)
if verbose
if flag1==1
disp(['Error in simul: No convergence inside GMRES after ' num2str(periods*10,'%6d') ' iterations, in block ' num2str(Blck_size,'%3d')]);
elseif flag1==2
disp(['Error in simul: Preconditioner is ill-conditioned, in block ' num2str(Blck_size,'%3d')]);
elseif flag1==3
disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block ' num2str(Blck_size,'%3d')]);
end
end
ilu_setup.droptol = ilu_setup.droptol/10;
reduced = 0;
else
dx = za - ya;
ya = ya + lambda*dx;
y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';
end
end
elseif stack_solve_algo==3
flag1=1;
while flag1>0
if preconditioner==2
[L1, U1]=ilu(g1a,ilu_setup);
[za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1);
elseif preconditioner==3
Size = Blck_size;
gss0 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
[L1, U1]=lu(gss0);
P1 = eye(size(gss0));
Q1 = eye(size(gss0));
L = kron(eye(periods),L1);
U = kron(eye(periods),U1);
P = kron(eye(periods),P1);
Q = kron(eye(periods),Q1);
[za,flag1] = bicgstab1(g1a,b,1e-7,Blck_size*periods,L,U, P, Q);
else
Size = Blck_size;
gss0 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
[L1, U1]=lu(gss0);
L1 = kron(eye(periods),L1);
U1 = kron(eye(periods),U1);
[za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1);
end
if flag1>0 || reduced
if verbose
if flag1==1
disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(periods*10,'%6d') ' iterations, in block ' num2str(Blck_size,'%3d')]);
elseif flag1==2
disp(['Error in simul: Preconditioner is ill-conditioned, in block ' num2str(Blck_size,'%3d')]);
elseif flag1==3
disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block ' num2str(Blck_size,'%3d')]);
end
end
ilu_setup.droptol = ilu_setup.droptol/10;
reduced = 0;
else
dx = za - ya;
ya = ya + lambda*dx;
y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';
end
end
elseif stack_solve_algo==4
ra = reshape(r(:, y_kmin+1:periods+y_kmin),periods*Blck_size, 1);
stpmx = 100 ;
stpmax = stpmx*max([sqrt(ya'*ya);size(y_index,2)]);
nn=1:size(ra,1);
g = (ra'*g1a)';
f = 0.5*ra'*ra;
p = -g1a\ra;
[yn,f,ra,check]=lnsrch1(ya,f,g,p,stpmax,'lnsrch1_wrapper_two_boundaries',nn,nn, options.solve_tolx, fname, y, y_index,x, params, steady_state, periods, y_kmin, Blck_size,options.periods);
dx = ya - yn;
y(1+y_kmin:periods+y_kmin,y_index)=reshape(yn',length(y_index),periods)';
end
end
iter=iter+1;
if verbose
disp(['iteration: ' num2str(iter,'%d') ' error: ' num2str(max_res,'%e')]);
end
end
if (iter>maxit_)
if verbose
printline(41)
%disp(['No convergence after ' num2str(iter,'%4d') ' iterations in Block ' num2str(Block_Num,'%d')])
end
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;
return
end
oo.deterministic_simulation.status = 1;
oo.deterministic_simulation.error = max_res;
oo.deterministic_simulation.iterations = iter;
oo.deterministic_simulation.block(Block_Num).status = 1;% Convergency obtained.
oo.deterministic_simulation.block(Block_Num).error = max_res;
oo.deterministic_simulation.block(Block_Num).iterations = iter;
|
