Commit 750f92ac authored by Sébastien Villemot's avatar Sébastien Villemot
Browse files

Cleanups related to the removal of leads and lags in deterministic model

(introduced in 60f76786)
- remove simk.m and clean simul.m
- document the change in the reference manual
parent fbbc21d5
......@@ -1043,7 +1043,8 @@ end;
<itemizedlist>
<listitem><para>the timing convention of <xref linkend="predetermined_variables"/> will have been changed to the default Dynare timing convention; in other words, variables declared as predetermined will be lagged on period back,</para></listitem>
<listitem><para>the <link linkend="expectation_operator">expectation operators</link> will have been removed, replaced by auxiliary variables and new equations as explained in the documentation of the operator,</para></listitem>
<listitem><para>for stochastic models, variables with leads or lags greater or equal than two will have been removed, replaced by new auxiliary variables and equations.</para></listitem>
<listitem><para>endogenous variables with leads or lags greater or equal than two will have been removed, replaced by new auxiliary variables and equations,</para></listitem>
<listitem><para>for a stochastic model, exogenous variables with leads or lags will also have been replaced by new auxiliary variables and equations.</para></listitem>
</itemizedlist>
</para>
......@@ -2194,7 +2195,7 @@ steady;
<refsect1 id="aux_vars">
<title>Auxiliary variables for leads and lags</title>
<para>For a stochastic model, Dynare will perform a transformation of the model so that there is only one lead and one lag on endogenous, and no lead/lag on exogenous.</para>
<para>Dynare will perform a transformation of the model so that there is only one lead and one lag on endogenous, and, in the case of a stochastic model, no lead/lag on exogenous.</para>
<para>This transformation is achieved by the creation of auxiliary variables, and corresponding equations. For example, if <literal>x(+2)</literal> exists in the model, Dynare will create one auxiliary variable <literal>AUX_ENDO_LEAD = x(+1)</literal>, and replace <literal>x(+2)</literal> by <literal>AUX_ENDO_LEAD(+1)</literal>.</para>
<para>A similar transformation is done for lags greater than 2 on endogenous (auxiliary variables will have a name beginning with <literal>AUX_ENDO_LAG</literal>), and for exogenous with leads and lags (auxiliary variables will have a name beginning with <literal>AUX_EXO_LEAG</literal> or <literal>AUX_EXO_LAG</literal> respectively).</para>
<para>Once created, all auxiliary variables are included in the set of endogenous variables. The output of decision rules (see below) is such that auxiliary variable names are replaced by the original variables they refer to.</para>
......
function simk
% function simk
% performs deterministic simulations with lead or lag on more than one
% period
%
% INPUTS
% ...
% OUTPUTS
% ...
% 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.
%
% Copyright (C) 1996-2009 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/>.
global M_ options_ oo_
global it_ iyr0 broyden_
%func_name = [M_.fname '_static'];
nk = M_.maximum_endo_lag + M_.maximum_endo_lead + 1 ;
ny = size(M_.lead_lag_incidence,2) ;
icc1 = M_.lead_lag_incidence(nk,:) > 0;
for i = 1:M_.maximum_lead -1
icc1 = [M_.lead_lag_incidence(nk-i,:) | icc1(1,:); icc1] ;
end
icc1 = find(icc1') ;
iy = M_.lead_lag_incidence > 0 ;
isc = cumsum(sum(iy',1))' ;
iyr0 = find(M_.lead_lag_incidence') ;
ncc1 = size(icc1,1) ;
ncc = ncc1 + 1 ;
ncs = size(iyr0,1) ;
ky = zeros(ny,nk) ; % indices of variables at each lead or lag
lky = zeros(nk,1) ;
for i = 1:nk
j = find(M_.lead_lag_incidence(i,:))' ;
if isempty(j)
lky(i) = 0;
else
lky(i) = size(j,1) ;
ky(1:lky(i),i) = j ;
end
end
jwc = find(iy(2:M_.maximum_endo_lead+1,:)') ; % indices of columns for
% triangularization
% as many rows as lags in model
if isempty(jwc)
jwc = 0 ;
ljwc = 0 ;
temp = icc1 ;
else
ljwc = size(jwc,1) ; % length of each row in jwc
temp = union(jwc,icc1) ; % prepares next iteration
end
j1 = ky(1:lky(1),1) ;
lj1 = lky(1) ;
for i = 2:M_.maximum_endo_lag
[j1,lj1] = ffill(j1,lj1,selif(temp+(i-1)*ny,temp <= ny)) ;
if M_.maximum_endo_lead == 1
if lky(i+M_.maximum_endo_lead) > 0
[jwc,ljwc] = ffill(jwc,ljwc, ky(1:lky(i+M_.maximum_endo_lead),i+M_.maximum_endo_lead)+(M_.maximum_endo_lead-1)*ny) ;
if ljwc(i) == 0
temp = icc1;
else
temp = union(jwc(1:ljwc(i),i),icc1) ;
end
else
[jwc,ljwc] = ffill(jwc,ljwc,[]) ;
temp = icc1 ;
end
else
temp = temp(lj1(i)+1:size(temp,1),:) - ny ;
if lky(i+M_.maximum_endo_lead) > 0
[jwc,ljwc] = ffill(jwc,ljwc,[temp;ky(1:lky(i+M_.maximum_endo_lead),i+M_.maximum_endo_lead)+(M_.maximum_endo_lead-1)*ny]);
else
[jwc,ljwc] = ffill(jwc,ljwc,temp) ;
end
temp = union(jwc(1:ljwc(i),i),icc1) ;
end
end
[j1,lj1] = ffill(j1,lj1,selif(temp+M_.maximum_endo_lag*ny, temp <= ny)) ;
ltemp = zeros(M_.maximum_endo_lag,1) ;
jwc1 = zeros(ncc1,M_.maximum_endo_lag) ;
for i = 1:M_.maximum_endo_lag
temp = union(jwc(1:ljwc(i),i),icc1) ;
ltemp(i) = size(temp,1) ;
if ljwc(i) > 0
jwc(1:ljwc(i),i) = indnv(jwc(1:ljwc(i),i),temp) ;
end
jwc1(:,i) = indnv(icc1,temp) ;
end
h1 = clock ;
disp (['-----------------------------------------------------']) ;
disp ('MODEL SIMULATION') ;
fprintf ('\n') ;
for iter = 1:options_.maxit_
h2 = clock ;
oo_.endo_simul = oo_.endo_simul(:);
err_f = 0;
fid = fopen([M_.fname '.swp'],'w+') ;
it_ = 1+M_.maximum_lag ;
ic = [1:ny] ;
iyr = iyr0 ;
i = M_.maximum_endo_lag+1 ;
while (i>1) & (it_<=options_.periods+M_.maximum_endo_lag)
h3 = clock ;
if broyden_ & iter > 1
%d1_ = -feval(fh,oo_.endo_simul(iyr));
d1 = -feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),z,oo_.exo_simul, M_.params, it_);
else
%jacob(func_name,oo_.endo_simul(iyr)) ;
[d1,jacobian] = feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),oo_.exo_simul, M_.params, it_);
d1 = -d1 ;
end
err_f = max(err_f,max(abs(d1)));
if lky(i) ~= 0
j1i = ky(1:lky(i),i) ;
w0 = jacobian(:,isc(i-1)+1:isc(i)) ;
else
w0 = [];
end
ttemp = iy(i+1:i+M_.maximum_endo_lead,:)' ;
jwci = find(ttemp) ;
if ~ isempty(jwci)
w = jacobian(:,isc(i)+1:isc(i+M_.maximum_endo_lead)) ;
end
j = i ;
while j <= M_.maximum_endo_lag
if ~isempty(w0)
ofs = ((it_-M_.maximum_lag-M_.maximum_endo_lag+j-2)*ny)*ncc*8 ;
junk = fseek(fid,ofs,-1) ;
c = fread(fid,[ncc,ny],'float64')';
if isempty(jwci)
w = -w0*c(j1i,1:ncc1) ;
jwci = icc1 ;
else
iz = union(jwci,icc1) ;
ix = indnv(jwci,iz) ;
iy__ = indnv(icc1,iz) ;
temp = zeros(size(w,1),size(iz,1)) ;
temp(:,ix) = w;
temp(:,iy__) = temp(:,iy__)-w0*c(j1i,1:ncc1) ;
w = temp ;
jwci = iz ;
clear temp iz ix iy__ ;
end
d1 = d1-w0*c(j1i,ncc) ;
clear c ;
end
j = j + 1 ;
if isempty(jwci)
j1i = [];
if lky(j+M_.maximum_endo_lead) ~= 0
jwci = ky(1:lky(j+M_.maximum_endo_lead),j+M_.maximum_endo_lead) + (M_.maximum_endo_lead-1)*ny ;
w = jacobian(:,isc(j+M_.maximum_endo_lead-1)+1:isc(j+M_.maximum_endo_lead)) ;
else
jwci = [] ;
end
else
j1i = selif(jwci,jwci<(ny+1)) ;
w0 = w(:,1:size(j1i,1)) ;
if size(jwci,1) == size(j1i,1)
if lky(j+M_.maximum_endo_lead) ~= 0
jwci = ky(1:lky(j+M_.maximum_endo_lead),j+M_.maximum_endo_lead)+(M_.maximum_endo_lead-1)*ny ;
w = jacobian(:,isc(j+M_.maximum_endo_lead-1)+1:isc(j+M_.maximum_endo_lead)) ;
else
jwci = [] ;
end
else
jwci = jwci(size(j1i,1)+1:size(jwci,1),:)-ny ;
w = w(:,size(j1i,1)+1:size(w,2)) ;
if lky(j+M_.maximum_endo_lead) ~= 0
jwci = [ jwci; ky(1: lky(j+M_.maximum_endo_lead),j+M_.maximum_endo_lead)+(M_.maximum_endo_lead-1)*ny] ;
w = [w jacobian(:,isc(j+M_.maximum_endo_lead-1)+1:isc(j+M_.maximum_endo_lead))] ;
% else
% jwci = [] ;
end
end
end
end
jwci = [indnv(jwci,icc1);ncc] ;
w = [w d1] ;
c = zeros(ny,ncc) ;
c(:,jwci) = w0\w ;
clear w w0 ;
junk = fseek(fid,0,1) ;
fwrite(fid,c','float64') ;
clear c ;
it_ = it_ + 1;
ic = ic + ny ;
iyr = iyr + ny ;
i = i - 1 ;
end
icr0 = (it_-M_.maximum_lag-M_.maximum_endo_lag -1)*ny ;
while it_ <= options_.periods+M_.maximum_lag
if broyden_
%d1_ = -feval(fh,oo_.endo_simul(iyr));
d1 = -feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),z,oo_.exo_simul, M_.params, it_);
else
%jacob(func_name,oo_.endo_simul(iyr)) ;
[d1,jacobian] = feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),oo_.exo_simul, M_.params, it_);
d1 = -d1 ;
end
err_f = max(err_f,max(abs(d1)));
w0 = jacobian(:,1:isc(1)) ;
w = jacobian(:,isc(1)+1:isc(1+M_.maximum_endo_lead)) ;
j = 1 ;
while j <= M_.maximum_endo_lag
icr = j1(1:lj1(j),j)-(j-1)*ny ;
ofs = ((icr0+(j-1)*ny+1)-1)*ncc*8 ;
junk = fseek(fid,ofs,-1) ;
c = fread(fid,[ncc,ny],'float64')' ;
temp = zeros(ny,ltemp(j)) ;
if ljwc(j) > 0
temp(:,jwc(1:ljwc(j),j)) = w ;
end
temp(:,jwc1(:,j))=temp(:,jwc1(:,j))-w0*c(icr,1:ncc1) ;
w = temp ;
clear temp ;
d1 = d1-w0*c(icr,ncc) ;
clear c ;
j = j + 1 ;
w0 = w(:,1:lj1(j)) ;
if M_.maximum_endo_lead == 1
w = jacobian(:,isc(j+M_.maximum_endo_lead-1)+1:isc(j+M_.maximum_endo_lead)) ;
else
w = w(:,lj1(j)+1:size(w,2)) ;
if lky(j+M_.maximum_endo_lead) > 0
w = [w jacobian(:,isc(j+M_.maximum_endo_lead-1)+1:isc(j+M_.maximum_endo_lead))] ;
end
end
end
c = w0\[w d1] ;
d1 = [] ;
clear w w0 ;
junk = fseek(fid,0,1) ;
fwrite(fid,c','float64') ;
clear c ;
it_ = it_ + 1 ;
ic = ic + ny ;
iyr = iyr + ny ;
icr0 = icr0 + ny ;
end
if options_.terminal_condition == 1
ofs = (((it_-M_.maximum_lag-2)*ny+1)-1)*ncc*8 ;
junk = fseek(fid,ofs,-1) ;
c = fread(fid,[ncc,ny],'float64')';
for i = 1:M_.maximum_endo_lead
w = tril(triu(ones(ny,ny+ncc1))) ;
w(:,jwc1(:,M_.maximum_endo_lag)) = w(:,jwc1(:,M_.maximum_endo_lag))+c(:,1:ncc1) ;
c = [w(:,ny+1:size(w,2))' c(:,ncc)]/w(:,1:ny) ;
junk = fseek(fid,0,1) ;
fwrite(fid,c','float64') ;
it_ = it_+1 ;
ic = ic + ny ;
end
end
oo_.endo_simul = reshape(oo_.endo_simul,ny,options_.periods+M_.maximum_lag+M_.maximum_endo_lead) ;
if options_.terminal_condition == 1
hbacsup = clock ;
c = bksupk(ny,fid,ncc,icc1) ;
hbacsup = etime(clock,hbacsup) ;
c = reshape(c,ny,options_.periods+M_.maximum_endo_lead)' ;
y(:,1+M_.maximum_endo_lag:(options_.periods+M_.maximum_endo_lead+M_.maximum_endo_lag)) = y(:,1+M_.maximum_endo_lag:(options_.periods+M_.maximum_endo_lead+M_.maximum_endo_lag))+options_.slowc*c' ;
else
hbacsup = clock ;
c = bksupk(ny,fid,ncc,icc1) ;
hbacsup = etime(clock,hbacsup) ;
c = reshape(c,ny,options_.periods)' ;
oo_.endo_simul(:,1+M_.maximum_endo_lag:(options_.periods+M_.maximum_endo_lag)) = oo_.endo_simul(:,1+M_.maximum_endo_lag:(options_.periods+M_.maximum_endo_lag))+options_.slowc*c' ;
end
fclose(fid) ;
h2 = etime(clock,h2) ;
[junk,i1] = max(abs(c));
[junk,i2] = max(junk);
disp(['variable ' M_.endo_names(i2,:) ' period ' num2str(i1(i2))])
err = max(max(abs(c./options_.scalv'))) ;
disp ([num2str(iter) '- err = ' num2str(err)]) ;
disp (['err_f = ' num2str(err_f)])
disp ([' Time of this iteration : ' num2str(h2)]) ;
if options_.timing
disp ([' Back substitution : ' num2str(hbacsup)]) ;
end
if err < options_.dynatol
h1 = etime(clock,h1) ;
fprintf ('\n') ;
disp ([' Total time of simulation : ' num2str(h1)]) ;
fprintf ('\n') ;
disp ([' Convergence achieved.']) ;
disp (['-----------------------------------------------------']) ;
fprintf ('\n') ;
return ;
end
end
disp(['WARNING : the maximum number of iterations is reached.']) ;
fprintf ('\n') ;
disp (['-----------------------------------------------------']) ;
return ;
% 2/11/99 MJ took out reshapel
function simul(dr)
% function simul(dr)
% computes simulations
function simul
% function simul
% Computes deterministic simulations
%
% INPUTS
% dr: structure of decision rules for stochastic simulations
% None
%
% OUTPUTS
% ...
......@@ -12,7 +12,7 @@ function simul(dr)
% SPECIAL REQUIREMENTS
% none
% Copyright (C) 1996-2007 Dynare Team
% Copyright (C) 1996-2010 Dynare Team
%
% This file is part of Dynare.
%
......@@ -29,10 +29,10 @@ function simul(dr)
% 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_ ys0_
global M_ options_ oo_
if size(M_.lead_lag_incidence,2)-nnz(M_.lead_lag_incidence(M_.maximum_endo_lag+1,:)) > 0
mess = ['DYNARE: error in model specification : variable ' M_.endo_names(find(M_.lead_lag_incidence(M_.maximum_lag+1,:)==0),:)] ;
mess = ['SIMUL: error in model specification : variable ' M_.endo_names(find(M_.lead_lag_incidence(M_.maximum_lag+1,:)==0),:)] ;
mess = [mess ' doesn''t appear as current variable.'] ;
error (mess) ;
end
......@@ -82,18 +82,15 @@ else
oo_.endo_simul=bytecode('dynamic');
else
if M_.maximum_endo_lead == 0
error('simul doesn''t simulate purely backward models')
error('SIMUL: purely backward models are not supported')
elseif M_.maximum_endo_lag == 0
error('SIMUL: purely forward models are not supported')
elseif M_.maximum_endo_lag == 1 && M_.maximum_endo_lead == 1
sim1 ;
else
simk ;
error('SIMUL: internal error of Dynare, contact the developpers')
end;
end;
end;
dyn2vec;
% 6/18/01 MJ added dyn2vec if 40 variables or less
% 01/16/03 MJ use dyn2vec whatever the number of variables
% 02/18/03 MJ added oo_.steady_state for calling simult
% 05/24/03 MJ added options_ and options_.periods
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment