Skip to content
GitLab
Commit f904a1a1 authored by ferhat's avatar ferhat
Browse files

- new Incidence_Matrix class 

- lead and lag on exogenous variables
- corrections in dr1_sparse and dr11_sparse
- minor corrections in simulate

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@2255 ac1d8469-bf42-47a9-8791-bf33cf982152
parent a02e66e1
Expand all Hide whitespace changes
Inline Side-by-side
matlab/disp_dr_sparse.m 0 → 100644
View file @ f904a1a1
function disp_dr_sparse(dr,order,var_list)
% Copyright (C) 2001 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_
nx = 0;
nu = 0;
k = [];
klag = [];
k1 = [];
nspred = 0;
for i=1:length(M_.block_structure.block)
nspred = nspred + M_.block_structure.block(i).dr.nspred;
end;
ghu = zeros(M_.endo_nbr, M_.exo_nbr*(M_.maximum_exo_lag+M_.maximum_exo_lead+1));
ghx = zeros(M_.endo_nbr, nspred);
for i=1:length(M_.block_structure.block)
nx = nx + size(M_.block_structure.block(i).dr.ghx,2);
% M_.block_structure.block(i).dr.ghx
% M_.block_structure.block(i).equation
% M_.block_structure.block(i).variable
ghx(M_.block_structure.block(i).equation, M_.block_structure.block(i).variable(find(M_.block_structure.block(i).lead_lag_incidence(1: M_.block_structure.block(i).maximum_endo_lag,:))) ) = M_.block_structure.block(i).dr.ghx;
if(M_.block_structure.block(i).exo_nbr)
nu = nu + size(M_.block_structure.block(i).dr.ghu,2);
ghu(M_.block_structure.block(i).equation, M_.block_structure.block(i).exogenous) = M_.block_structure.block(i).dr.ghu;
end
k_tmp = find(M_.block_structure.block(i).dr.kstate(:,2) <= M_.block_structure.block(i).maximum_lag+1);
k = [k ; k_tmp];
klag = [klag ; M_.block_structure.block(i).dr.kstate(k_tmp,[1 2])];
k1 = [k1 ; M_.block_structure.block(i).variable(M_.block_structure.block(i).dr.order_var)'];
end
nvar = size(var_list,1);
if nvar == 0
nvar = length(k1);
ivar = [1:nvar];
else
ivar=zeros(nvar,1);
for i=1:nvar
i_tmp = strmatch(var_list(i,:),M_.endo_names(k1,:),'exact');
if isempty(i_tmp)
disp(var_list(i,:));
error (['One of the variable specified does not exist']) ;
else
ivar(i) = i_tmp;
end
end
end
disp('POLICY AND TRANSITION FUNCTIONS')
% variable names
str = ' ';
for i=1:nvar
str = [str sprintf('%16s',M_.endo_names(k1(ivar(i)),:))];
end
disp(str);
%
% constant
%
str = 'Constant ';
flag = 0;
for i=1:nvar
x = dr.ys(k1(ivar(i)));
if order > 1
x = x + dr.ghs2(ivar(i))/2;
end
if abs(x) > 1e-6
flag = 1;
str = [str sprintf('%16.6f',x)];
else
str = [str ' 0'];
end
end
if flag
disp(str)
end
if order > 1
str = '(correction) ';
flag = 0;
for i=1:nvar
x = dr.ghs2(ivar(i))/2;
if abs(x) > 1e-6
flag = 1;
str = [str sprintf('%16.6f',x)];
else
str = [str ' 0'];
end
end
if flag
disp(str)
end
end
%
% ghx
%
for k=1:nx
flag = 0;
str1 = sprintf('%s(%d)',M_.endo_names(k1(klag(k,1)),:),klag(k,2)-M_.maximum_lag-2);
str = sprintf('%-20s',str1);
for i=1:nvar
x = ghx(ivar(i),k);
if abs(x) > 1e-6
flag = 1;
str = [str sprintf('%16.6f',x)];
else
str = [str ' 0'];
end
end
if flag
disp(str)
end
end
%
% ghu
%
for k=1:nu
flag = 0;
str = sprintf('%-20s',M_.exo_names(k,:));
for i=1:nvar
x = ghu(ivar(i),k);
if abs(x) > 1e-6
flag = 1;
str = [str sprintf('%16.6f',x)];
else
str = [str ' 0'];
end
end
if flag
disp(str)
end
end
if order > 1
% ghxx
for k = 1:nx
for j = 1:k
flag = 0;
str1 = sprintf('%s(%d),%s(%d)',M_.endo_names(k1(klag(k,1)),:),klag(k,2)-M_.maximum_lag-2, ...
M_.endo_names(k1(klag(j,1)),:),klag(j,2)-M_.maximum_lag-2);
str = sprintf('%-20s',str1);
for i=1:nvar
if k == j
x = dr.ghxx(ivar(i),(k-1)*nx+j)/2;
else
x = dr.ghxx(ivar(i),(k-1)*nx+j);
end
if abs(x) > 1e-6
flag = 1;
str = [str sprintf('%16.6f',x)];
else
str = [str ' 0'];
end
end
if flag
disp(str)
end
end
end
%
% ghuu
%
for k = 1:nu
for j = 1:k
flag = 0;
str = sprintf('%-20s',[M_.exo_names(k,:) ',' M_.exo_names(j,:)] );
for i=1:nvar
if k == j
x = dr.ghuu(ivar(i),(k-1)*nu+j)/2;
else
x = dr.ghuu(ivar(i),(k-1)*nu+j);
end
if abs(x) > 1e-6
flag = 1;
str = [str sprintf('%16.6f',x)];
else
str = [str ' 0'];
end
end
if flag
disp(str)
end
end
end
%
% ghxu
%
for k = 1:nx
for j = 1:nu
flag = 0;
str1 = sprintf('%s(%d),%s',M_.endo_names(k1(klag(k,1)),:),klag(k,2)-M_.maximum_lag-2, ...
M_.exo_names(j,:));
str = sprintf('%-20s',str1);
for i=1:nvar
x = dr.ghxu(ivar(i),(k-1)*nu+j);
if abs(x) > 1e-6
flag = 1;
str = [str sprintf('%16.6f',x)];
else
str = [str ' 0'];
end
end
if flag
disp(str)
end
end
end
end
% $$$ dr.ghx
% $$$ dr.ghu
% $$$ dr.ghxx
% $$$ dr.ghuu
% $$$ dr.ghxu
% 01/08/2001 MJ added test for order in printing quadratic terms
% 02/21/2001 MJ pass all variable names through deblank()
% 02/21/2001 MJ changed from f to g format to write numbers
% 10/09/2002 MJ corrected error on constant whith subset of variables
matlab/dr11_sparse.m
View file @ f904a1a1
......@@ -17,7 +17,7 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
%task
info = 0;
klen = M_.maximum_endo_lag + M_.maximum_endo_lead + 1;
kstate = dr.kstate;
......@@ -46,7 +46,12 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi
m = m+length(k);
end
if M_.exo_nbr & task~=1
jacobia_
jacobia_(:,nz+1:end)
b
dr.ghu = -b\jacobia_(:,nz+1:end);
disp(['nz=' int2str(nz) ]);
dr.ghu
end
dr.eigval = eig(transition_matrix(dr,M_));
dr.rank = 0;
......@@ -156,12 +161,15 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi
k1 = find(kstate(n4:nd,2) == M_.maximum_endo_lag+1);
k2 = find(kstate(1:n3,2) == M_.maximum_endo_lag+2);
hx(k1,:)
gx(k2(nboth+1:end),:)
dr.ghx = [hx(k1,:); gx(k2(nboth+1:end),:)];
dr.ghx
%lead variables actually present in the model
j3 = nonzeros(kstate(:,3));
j4 = find(kstate(:,3));
% derivatives with respect to exogenous variables
disp(['M_.exo_nbr=' int2str(M_.exo_nbr)]);
if M_.exo_nbr
fu = aa(:,nz+(1:M_.exo_nbr));
a1 = b;
......@@ -216,7 +224,6 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi
dr.ghud{i} = -M2*dr.ghud{i-1}(end-nyf+1:end,:);
end
end
disp('end0');
if options_.order == 1
return
end
......@@ -479,5 +486,4 @@ function [dr,info,M_,options_,oo_] = dr11_sparse(dr,task,M_,options_,oo_, jacobi
end
end
disp('end');
end
\ No newline at end of file
matlab/dr1_sparse.m
View file @ f904a1a1
......@@ -210,23 +210,47 @@ function [dr,info,M_,options_,oo_] = dr1_sparse(dr,task,M_,options_,oo_)
if options_.order == 1
[junk,jacobia_] = feval([M_.fname '_dynamic'],ones(M_.maximum_lag+M_.maximum_lead+1,1)*dr.ys',[oo_.exo_simul ...
oo_.exo_det_simul], it_);
%full(jacobia_)
dr.eigval = [];
dr.nyf = 0;
dr.rank = 0;
first_col_exo = M_.endo_nbr * (M_.maximum_endo_lag + M_.maximum_endo_lead + 1);
for i=1:length(M_.block_structure.block)
%disp(['block = ' int2str(i)]);
M_.block_structure.block(i).dr.Null=0;
M_.block_structure.block(i).dr=set_state_space(M_.block_structure.block(i).dr,M_.block_structure.block(i));
col_selector=repmat(M_.block_structure.block(i).variable,1,M_.block_structure.block(i).maximum_endo_lag+M_.block_structure.block(i).maximum_endo_lead+1)+kron([M_.maximum_endo_lag-M_.block_structure.block(i).maximum_endo_lag:M_.maximum_endo_lag+M_.block_structure.block(i).maximum_endo_lead],M_.endo_nbr*ones(1,M_.block_structure.block(i).endo_nbr));
row_selector = M_.block_structure.block(i).equation;
%col_selector
jcb_=jacobia_(row_selector,col_selector);
jcb_ = jcb_(:,find(M_.block_structure.block(i).lead_lag_incidence')) ;
jcb_ = jcb_(:,find(M_.block_structure.block(i).lead_lag_incidence')) ;
if M_.block_structure.block(i).exo_nbr>0
col_selector = [ first_col_exo + ...
repmat(M_.block_structure.block(i).exogenous,1,M_.block_structure.block(i).maximum_exo_lag+M_.block_structure.block(i).maximum_exo_lead+1)+kron([M_.maximum_exo_lag-M_.block_structure.block(i).maximum_exo_lag:M_.maximum_exo_lag+M_.block_structure.block(i).maximum_exo_lead],M_.exo_nbr*ones(1,M_.block_structure.block(i).exo_nbr))];
end
%col_selector
jcb_ = [ jcb_ jacobia_(row_selector,col_selector)];
%full(jcb_)
hss_=0; %hessian(M_.block_structure.block(i).equation,M_.block_structure.block(i).variable);
dra = M_.block_structure.block(i).dr;
M_.block_structure.block(i).exo_nbr=M_.exo_nbr;
%M_.block_structure.block(i).exo_nbr=M_.exo_nbr;
[dra ,info,M_.block_structure.block(i),options_,oo_] = dr11_sparse(dra ,task,M_.block_structure.block(i),options_,oo_, jcb_, hss_);
M_.block_structure.block(i).dr = dra;
dr.eigval = [dr.eigval; dra.eigval];
dr.nyf = dr.nyf + nnz(dra.kstate(:,2)>M_.block_structure.block(i).maximum_endo_lag+1);
nyf = nnz(dra.kstate(:,2)>M_.block_structure.block(i).maximum_endo_lag+1);
n_explod = nnz(abs(dra.eigval) > options_.qz_criterium);
if nyf ~= n_explod
disp(['EIGENVALUES in block ' int2str(i) ':']);
[m_lambda,ii]=sort(abs(dra.eigval));
disp(sprintf('%16s %16s %16s\n','Modulus','Real','Imaginary'))
z=[m_lambda real(dra.eigval(ii)) imag(dra.eigval(ii))]';
disp(sprintf('%16.4g %16.4g %16.4g\n',z))
disp(['The rank condition is not satisfy in block ' int2str(i) ' :']);
disp([' ' int2str(nyf) ' forward-looking variable(s) for ' ...
int2str(n_explod) ' eigenvalue(s) larger than 1 in modulus']);
end
dr.nyf = dr.nyf + nyf;
dr.rank = dr.rank + dra.rank;
end;
end
......
matlab/dynare_m.exe
View file @ f904a1a1
No preview for this file type
matlab/stoch_simul_sparse.m 0 → 100644
View file @ f904a1a1
function info=stoch_simul(var_list)
% Copyright (C) 2001-2008 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_ it_
options_old = options_;
if options_.linear
options_.order = 1;
end
if (options_.order == 1)
options_.replic = 1;
end
TeX = options_.TeX;
iter_ = max(options_.periods,1);
if M_.exo_nbr > 0
oo_.exo_simul= ones(iter_ + M_.maximum_lag + M_.maximum_lead,1) * oo_.exo_steady_state';
end
check_model;
[oo_.dr, info] = resol(oo_.steady_state,0);
if info(1)
options_ = options_old;
print_info(info);
return
end
oo_dr_kstate = [];
oo_dr_nstatic = 0;
for i=1:length(M_.block_structure.block)
oo_dr_kstate = [ oo_dr_kstate ; M_.block_structure.block(i).dr.kstate];
oo_dr_nstatic = oo_dr_nstatic + M_.block_structure.block(i).dr.nstatic;
end
if ~options_.noprint
disp(' ')
disp('MODEL SUMMARY')
disp(' ')
disp([' Number of variables: ' int2str(M_.endo_nbr)])
disp([' Number of stochastic shocks: ' int2str(M_.exo_nbr)])
disp([' Number of state variables: ' ...
int2str(length(find(oo_dr_kstate(:,2) <= M_.maximum_lag+1)))])
disp([' Number of jumpers: ' ...
int2str(length(find(oo_dr_kstate(:,2) == M_.maximum_lag+2)))])
disp([' Number of static variables: ' int2str(oo_dr_nstatic)])
my_title='MATRIX OF COVARIANCE OF EXOGENOUS SHOCKS';
labels = deblank(M_.exo_names);
headers = strvcat('Variables',labels);
lh = size(labels,2)+2;
table(my_title,headers,labels,M_.Sigma_e,lh,10,6);
disp(' ')
disp_dr_sparse(oo_.dr,options_.order,var_list);
end
if options_.simul == 0 & options_.nomoments == 0
disp_th_moments(oo_.dr,var_list);
elseif options_.simul == 1
if options_.periods == 0
error('STOCH_SIMUL error: number of periods for the simulation isn''t specified')
end
if options_.periods < options_.drop
disp(['STOCH_SIMUL error: The horizon of simulation is shorter' ...
' than the number of observations to be DROPed'])
options_ =options_old;
return
end
oo_.endo_simul = simult(repmat(oo_.dr.ys,1,M_.maximum_lag),oo_.dr);
dyn2vec;
if options_.nomoments == 0
disp_moments(oo_.endo_simul,var_list);
end
end
if options_.irf
n = size(var_list,1);
if n == 0
n = M_.endo_nbr;
ivar = [1:n]';
var_list = M_.endo_names;
if TeX
var_listTeX = M_.endo_names_tex;
end
else
ivar=zeros(n,1);
if TeX
var_listTeX = [];
end
for i=1:n
i_tmp = strmatch(var_list(i,:),M_.endo_names,'exact');
if isempty(i_tmp)
error (['One of the specified variables does not exist']) ;
else
ivar(i) = i_tmp;
if TeX
var_listTeX = strvcat(var_listTeX,deblank(M_.endo_names_tex(i_tmp,:)));
end
end
end
end
if TeX
fidTeX = fopen([M_.fname '_IRF.TeX'],'w');
fprintf(fidTeX,'%% TeX eps-loader file generated by stoch_simul.m (Dynare).\n');
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
end
olditer = iter_;% Est-ce vraiment utile ? Il y a la mme ligne dans irf...
SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord)=M_.Sigma_e+1e-14*eye(M_.exo_nbr);
cs = transpose(chol(SS));
tit(M_.exo_names_orig_ord,:) = M_.exo_names;
if TeX
titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex;
end
for i=1:M_.exo_nbr
if SS(i,i) > 1e-13
y=irf(oo_.dr,cs(M_.exo_names_orig_ord,i), options_.irf, options_.drop, ...
options_.replic, options_.order);
if options_.relative_irf
y = 100*y/cs(i,i);
end
irfs = [];
mylist = [];
if TeX
mylistTeX = [];
end
for j = 1:n
assignin('base',[deblank(M_.endo_names(ivar(j),:)) '_' deblank(M_.exo_names(i,:))],...
y(ivar(j),:)');
eval(['oo_.irfs.' deblank(M_.endo_names(ivar(j),:)) '_' ...
deblank(M_.exo_names(i,:)) ' = y(ivar(j),:);']);
if max(y(ivar(j),:)) - min(y(ivar(j),:)) > 1e-10
irfs = cat(1,irfs,y(ivar(j),:));
mylist = strvcat(mylist,deblank(var_list(j,:)));
if TeX
mylistTeX = strvcat(mylistTeX,deblank(var_listTeX(j,:)));
end
end
end
if options_.nograph == 0
number_of_plots_to_draw = size(irfs,1);
[nbplt,nr,nc,lr,lc,nstar] = pltorg(number_of_plots_to_draw);
if nbplt == 0
elseif nbplt == 1
if options_.relative_irf
hh = figure('Name',['Relative response to' ...
' orthogonalized shock to ' tit(i,:)]);
else
hh = figure('Name',['Orthogonalized shock to' ...
' ' tit(i,:)]);
end
for j = 1:number_of_plots_to_draw
subplot(nr,nc,j);
plot(1:options_.irf,transpose(irfs(j,:)),'-k','linewidth',1);
hold on
plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
hold off
xlim([1 options_.irf]);
title(deblank(mylist(j,:)),'Interpreter','none');
end
eval(['print -depsc2 ' M_.fname '_IRF_' deblank(tit(i,:)) '.eps']);
if ~exist('OCTAVE_VERSION')
eval(['print -dpdf ' M_.fname '_IRF_' deblank(tit(i,:))]);
saveas(hh,[M_.fname '_IRF_' deblank(tit(i,:)) '.fig']);
end
if TeX
fprintf(fidTeX,'\\begin{figure}[H]\n');
for j = 1:number_of_plots_to_draw
fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{$%s$}\n'],deblank(mylist(j,:)),deblank(mylistTeX(j,:)));
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_IRF_%s}\n',M_.fname,deblank(tit(i,:)));
fprintf(fidTeX,'\\caption{Impulse response functions (orthogonalized shock to $%s$).}',titTeX(i,:));
fprintf(fidTeX,'\\label{Fig:IRF:%s}\n',deblank(tit(i,:)));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
% close(hh)
else
for fig = 1:nbplt-1
if options_.relative_irf == 1
hh = figure('Name',['Relative response to orthogonalized shock' ...
' to ' tit(i,:) ' figure ' int2str(fig)]);
else
hh = figure('Name',['Orthogonalized shock to ' tit(i,:) ...
' figure ' int2str(fig)]);
end
for plt = 1:nstar
subplot(nr,nc,plt);
plot(1:options_.irf,transpose(irfs((fig-1)*nstar+plt,:)),'-k','linewidth',1);
hold on
plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
hold off
xlim([1 options_.irf]);
title(deblank(mylist((fig-1)*nstar+plt,:)),'Interpreter','none');
end
eval(['print -depsc2 ' M_.fname '_IRF_' deblank(tit(i,:)) int2str(fig) '.eps']);
if ~exist('OCTAVE_VERSION')
eval(['print -dpdf ' M_.fname '_IRF_' deblank(tit(i,:)) int2str(fig)]);
saveas(hh,[M_.fname '_IRF_' deblank(tit(i,:)) int2str(fig) '.fig']);
end
if TeX
fprintf(fidTeX,'\\begin{figure}[H]\n');
for j = 1:nstar
fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{$%s$}\n'],deblank(mylist((fig-1)*nstar+j,:)),deblank(mylistTeX((fig-1)*nstar+j,:)));
end
fprintf(fidTeX,'\\centering \n');
fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_IRF_%s%s}\n',M_.fname,deblank(tit(i,:)),int2str(fig));
if options_.relative_irf
fprintf(fidTeX,['\\caption{Relative impulse response' ...
' functions (orthogonalized shock to $%s$).}'],deblank(titTeX(i,:)));
else
fprintf(fidTeX,['\\caption{Impulse response functions' ...
' (orthogonalized shock to $%s$).}'],deblank(titTeX(i,:)));
end
fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s:%s}\n',deblank(tit(i,:)),int2str(fig));
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
% close(hh);
end
hh = figure('Name',['Orthogonalized shock to ' tit(i,:) ' figure ' int2str(nbplt) '.']);
m = 0;
for plt = 1:number_of_plots_to_draw-(nbplt-1)*nstar;
m = m+1;
subplot(lr,lc,m);
plot(1:options_.irf,transpose(irfs((nbplt-1)*nstar+plt,:)),'-k','linewidth',1);
hold on
plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
hold off