sim1_lbj.m 4.5 KB
Newer Older
Houtan Bastani's avatar
Houtan Bastani committed
1
function sim1_lbj()
2
3
4
% function sim1_lbj
% performs deterministic simulations with lead or lag on one period
% using the historical LBJ algorithm
5
6
7
8
9
10
%
% INPUTS
%   ...
% OUTPUTS
%   ...
% ALGORITHM
11
12
13
14
%   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.
15
16
17
18
%
% SPECIAL REQUIREMENTS
%   None.

Sébastien Villemot's avatar
Sébastien Villemot committed
19
% Copyright (C) 1996-2012 Dynare Team
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
%
% 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/>.

global M_ options_ oo_

lead_lag_incidence = M_.lead_lag_incidence;

40
ny = size(oo_.endo_simul,1) ;
41
42
nyp = nnz(lead_lag_incidence(1,:)) ;
nyf = nnz(lead_lag_incidence(3,:)) ;
43
nrs = ny+nyp+nyf+1 ;
44
nrc = nyf+1 ;
45
46
iyf = find(lead_lag_incidence(3,:)>0) ;
iyp = find(lead_lag_incidence(1,:)>0) ;
47
48
49
50
51
52
53
54
55
56
57
isp = [1:nyp] ;
is = [nyp+1:ny+nyp] ;
isf = iyf+nyp ;
isf1 = [nyp+ny+1:nyf+nyp+ny+1] ;
stop = 0 ;
iz = [1:ny+nyp+nyf];

disp (['-----------------------------------------------------']) ;
disp (['MODEL SIMULATION :']) ;
fprintf('\n') ;

58
it_init = M_.maximum_lag+1 ;
59
60

h1 = clock ;
61
for iter = 1:options_.simul.maxit
62
63
    h2 = clock ;
    
64
65
66
67
68
    if options_.terminal_condition == 0
        c = zeros(ny*options_.periods,nrc) ;
    else
        c = zeros(ny*(options_.periods+1),nrc) ;
    end
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
    it_ = it_init ;
    z = [oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1)] ;
    [d1,jacobian] = feval([M_.fname '_dynamic'],z,oo_.exo_simul, M_.params, oo_.steady_state,it_);
    jacobian = [jacobian(:,iz) -d1] ;
    ic = [1:ny] ;
    icp = iyp ;
    c (ic,:) = jacobian(:,is)\jacobian(:,isf1) ;
    for it_ = it_init+(1:options_.periods-1)
        z = [oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1)] ;
        [d1,jacobian] = feval([M_.fname '_dynamic'],z,oo_.exo_simul, ...
                              M_.params, oo_.steady_state, it_);
        jacobian = [jacobian(:,iz) -d1] ;
        jacobian(:,[isf nrs]) = jacobian(:,[isf nrs])-jacobian(:,isp)*c(icp,:) ;
        ic = ic + ny ;
        icp = icp + ny ;
        c (ic,:) = jacobian(:,is)\jacobian(:,isf1) ;
    end
    
    if options_.terminal_condition == 1
        s = eye(ny) ;
        s(:,isf) = s(:,isf)+c(ic,1:nyf) ;
        ic = ic + ny ;
        c(ic,nrc) = s\c(ic,nrc) ;
        c = bksup1(c,ny,nrc,iyf,options_.periods) ;
        c = reshape(c,ny,options_.periods+1) ;
        oo_.endo_simul(:,it_init+(0:options_.periods)) = oo_.endo_simul(:,it_init+(0:options_.periods))+options_.slowc*c ;
    else
        c = bksup1(c,ny,nrc,iyf,options_.periods) ;
        c = reshape(c,ny,options_.periods) ;
        oo_.endo_simul(:,it_init+(0:options_.periods-1)) = oo_.endo_simul(:,it_init+(0:options_.periods-1))+options_.slowc*c ;
100
    end
101
102
103
104
    
    err = max(max(abs(c./options_.scalv')));
    disp([num2str(iter) ' -     err = ' num2str(err)]) ;
    disp(['     Time of iteration       :' num2str(etime(clock,h2))]) ;
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
    
    if err < options_.dynatol.f
        stop = 1 ;
        fprintf('\n') ;
        disp([' Total time of simulation        :' num2str(etime(clock,h1))]) ;
        fprintf('\n') ;
        disp([' Convergency obtained.']) ;
        fprintf('\n') ;
        oo_.deterministic_simulation.status = 1;% Convergency obtained.
        oo_.deterministic_simulation.error = err;
        oo_.deterministic_simulation.iterations = iter;
        break
    end
end

if ~stop
    fprintf('\n') ;
    disp(['     Total time of simulation        :' num2str(etime(clock,h1))]) ;
    fprintf('\n') ;
124
    disp(['WARNING : maximum number of iterations is reached (modify options_.simul.maxit).']) ;
125
126
127
128
    fprintf('\n') ;
    oo_.deterministic_simulation.status = 0;% more iterations are needed.
    oo_.deterministic_simulation.error = err;
    oo_.deterministic_simulation.errors = c/abs(err);    
129
    oo_.deterministic_simulation.iterations = options_.simul.maxit;
130
131
132
end
disp (['-----------------------------------------------------']) ;