Skip to content
Snippets Groups Projects
Commit f6ce8926 authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files

Adapted posterior IRFs routines.

parent b11f6e25
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
matlab/PosteriorIRF.m
+ 35
− 37
View file @ f6ce8926
function PosteriorIRF(type)
% Builds posterior IRFs after the MH algorithm.
%
% INPUTS
% Builds posterior IRFs after the MH algorithm.
%
% INPUTS
% o type [char] string specifying the joint density of the
% deep parameters ('prior','posterior').
%
% OUTPUTS
% deep parameters ('prior','posterior').
%
% OUTPUTS
% None (oo_ and plots).
%
% SPECIAL REQUIREMENTS
......@@ -34,9 +34,10 @@ function PosteriorIRF(type)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_ estim_params_ oo_ M_ bayestopt_ dataset_
global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info
% Set the number of periods
if isempty(options_.irf) || ~options_.irf
if isempty(options_.irf) || ~options_.irf
options_.irf = 40;
end
% Set varlist if necessary
......@@ -58,7 +59,7 @@ end
% Get index of shocks for requested IRFs
irf_shocks_indx = getIrfShocksIndx();
% Set various parameters & Check or create directories
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
......@@ -68,8 +69,8 @@ np = estim_params_.np ;
npar = nvx+nvn+ncx+ncn+np;
offset = npar-np; clear('nvx','nvn','ncx','ncn','np');
nvobs = dataset_.info.nvobs;
gend = dataset_.info.ntobs;
nvobs = dataset_.vobs;
gend = dataset_.nobs;
MaxNumberOfPlotPerFigure = 9;
nn = sqrt(MaxNumberOfPlotPerFigure);
MAX_nirfs_dsge = ceil(options_.MaxNumberOfBytes/(options_.irf*nvar*M_.exo_nbr)/8);
......@@ -111,14 +112,14 @@ else% type = 'prior'
B = options_.prior_draws;
options_.B = B;
end
try
try
delete([MhDirectoryName filesep M_.fname '_irf_dsge*.mat'])
catch
catch
disp('No _IRFs (dsge) files to be deleted!')
end
try
try
delete([MhDirectoryName filesep M_.fname '_irf_bvardsge*.mat'])
catch
catch
disp('No _IRFs (bvar-dsge) files to be deleted!')
end
irun = 0;
......@@ -144,10 +145,6 @@ if MAX_nirfs_dsgevar
else
stock_irf_bvardsge = zeros(options_.irf,nvobs,M_.exo_nbr,B);
end
[mYY,mXY,mYX,mXX,Ydata,Xdata] = ...
var_sample_moments(options_.first_obs,options_.first_obs+options_.nobs-1,...
options_.dsge_varlag,-1,options_.datafile,options_.varobs,...
options_.xls_sheet,options_.xls_range);
NumberOfLags = options_.dsge_varlag;
NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
if options_.noconstant
......@@ -174,7 +171,7 @@ localVars.npar = npar;
localVars.type=type;
if strcmpi(type,'posterior')
while b<B
while b<B
b = b + 1;
x(b,:) = GetOneDraw(type);
end
......@@ -192,7 +189,7 @@ if options_.dsge_var
localVars.NumberOfLags = options_.dsge_varlag;
localVars.NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
localVars.Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
end
end
localVars.nvar=nvar;
localVars.IndxVariables=IndxVariables;
localVars.MAX_nirfs_dsgevar=MAX_nirfs_dsgevar;
......@@ -225,14 +222,15 @@ else
localVars.NumberOfIRFfiles_dsge=NumberOfIRFfiles_dsge;
localVars.NumberOfIRFfiles_dsgevar=NumberOfIRFfiles_dsgevar;
localVars.ifil2=ifil2;
globalVars = struct('M_',M_, ...
'options_', options_, ...
'bayestopt_', bayestopt_, ...
'estim_params_', estim_params_, ...
'oo_', oo_, ...
'dataset_',dataset_);
'dataset_',dataset_, ...
'dataset_info',dataset_info);
% which files have to be copied to run remotely
NamFileInput(1,:) = {'',[M_.fname '_static.m']};
NamFileInput(2,:) = {'',[M_.fname '_dynamic.m']};
......@@ -244,13 +242,13 @@ else
for j=1:length(fout),
nosaddle = nosaddle + fout(j).nosaddle;
end
end
% END first parallel section!
if nosaddle
disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))])
disp(['PosteriorIRF :: Percentage of discarded posterior draws = ' num2str(nosaddle/(B+nosaddle))])
end
ReshapeMatFiles('irf_dsge',type)
......@@ -323,7 +321,7 @@ if MAX_nirfs_dsgevar
MedianIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr);
VarIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr);
DistribIRFdsgevar = zeros(options_.irf,9,nvar,M_.exo_nbr);
HPDIRFdsgevar = zeros(options_.irf,2,nvar,M_.exo_nbr);
HPDIRFdsgevar = zeros(options_.irf,2,nvar,M_.exo_nbr);
fprintf('Estimation::mcmc: Posterior (bvar-dsge) IRFs...\n');
tit(M_.exo_names_orig_ord,:) = M_.exo_names;
kdx = 0;
......@@ -341,7 +339,7 @@ if MAX_nirfs_dsgevar
end
kdx = kdx + size(STOCK_IRF_BVARDSGE,1);
end
clear STOCK_IRF_BVARDSGE;
clear STOCK_IRF_BVARDSGE;
for i = irf_shocks_indx
for j = 1:nvar
name = [deblank(M_.endo_names(IndxVariables(j),:)) '_' deblank(tit(i,:))];
......@@ -384,7 +382,7 @@ localVars.MaxNumberOfPlotPerFigure=MaxNumberOfPlotPerFigure;
if options_.dsge_var
localVars.HPDIRFdsgevar=HPDIRFdsgevar;
localVars.MeanIRFdsgevar = MeanIRFdsgevar;
end
end
% The files .TeX are genereted in sequential way always!
......@@ -394,14 +392,14 @@ if options_.TeX
fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
fprintf(fidTeX,' \n');
titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex;
for i=irf_shocks_indx
NAMES = [];
TEXNAMES = [];
for j=1:nvar
if max(abs(MeanIRF(:,j,i))) > options_.impulse_responses.plot_threshold
if max(abs(MeanIRF(:,j,i))) > options_.impulse_responses.plot_threshold
name = deblank(varlist(j,:));
texname = deblank(varlist_TeX(j,:));
......@@ -413,7 +411,7 @@ if options_.TeX
TEXNAMES = char(TEXNAMES,['$' texname '$']);
end
end
end
fprintf(fidTeX,'\\begin{figure}[H]\n');
for jj = 1:size(TEXNAMES,1)
......@@ -430,10 +428,10 @@ if options_.TeX
fprintf(fidTeX,'\\end{figure}\n');
fprintf(fidTeX,' \n');
end
fprintf(fidTeX,'%% End of TeX file.\n');
fclose(fidTeX);
end
% The others file format are generated in parallel by PosteriorIRF_core2!
......@@ -453,7 +451,7 @@ if ~isoctave
else
globalVars = struct('M_',M_, ...
'options_', options_);
[fout] = masterParallel(options_.parallel, 1, M_.exo_nbr,NamFileInput,'PosteriorIRF_core2', localVars, globalVars, options_.parallel_info);
end
end
......
matlab/PosteriorIRF_core1.m
+ 16
− 34
View file @ f6ce8926
......@@ -41,7 +41,7 @@ function myoutput=PosteriorIRF_core1(myinputs,fpar,B,whoiam, ThisMatlab)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_ estim_params_ oo_ M_ bayestopt_ dataset_
global options_ estim_params_ oo_ M_ bayestopt_ dataset_ dataset_info
if nargin<4,
whoiam=0;
......@@ -189,24 +189,24 @@ while fpar<B
end
if MAX_nirfs_dsgevar
IRUN = IRUN+1;
[fval,junk1,junk2,cost_flag,info,PHI,SIGMAu,iXX] = dsge_var_likelihood(deep',dataset_,options_,M_,estim_params_,bayestopt_,oo_);
[fval,junk1,junk2,cost_flag,info,PHI,SIGMAu,iXX] = dsge_var_likelihood(deep',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,oo_);
dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
DSGE_PRIOR_WEIGHT = floor(dataset_.info.ntobs*(1+dsge_prior_weight));
SIGMA_inv_upper_chol = chol(inv(SIGMAu*dataset_.info.ntobs*(dsge_prior_weight+1)));
DSGE_PRIOR_WEIGHT = floor(dataset_.nobs*(1+dsge_prior_weight));
SIGMA_inv_upper_chol = chol(inv(SIGMAu*dataset_.nobs*(dsge_prior_weight+1)));
explosive_var = 1;
while explosive_var
% draw from the marginal posterior of SIGMA
SIGMAu_draw = rand_inverse_wishart(dataset_.info.nvobs, DSGE_PRIOR_WEIGHT-NumberOfParametersPerEquation, ...
SIGMAu_draw = rand_inverse_wishart(dataset_.vobs, DSGE_PRIOR_WEIGHT-NumberOfParametersPerEquation, ...
SIGMA_inv_upper_chol);
% draw from the conditional posterior of PHI
PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,dataset_.info.nvobs, PHI, ...
PHI_draw = rand_matrix_normal(NumberOfParametersPerEquation,dataset_.vobs, PHI, ...
chol(SIGMAu_draw)', chol(iXX)');
Companion_matrix(1:dataset_.info.nvobs,:) = transpose(PHI_draw(1:NumberOfLagsTimesNvobs,:));
Companion_matrix(1:dataset_.vobs,:) = transpose(PHI_draw(1:NumberOfLagsTimesNvobs,:));
% Check for stationarity
explosive_var = any(abs(eig(Companion_matrix))>1.000000001);
end
% Get the mean
mu = zeros(1,dataset_.info.nvobs);
mu = zeros(1,dataset_.vobs);
% Get rotation
if dsge_prior_weight > 0
Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e);
......@@ -216,24 +216,24 @@ while fpar<B
SIGMAu_chol = chol(SIGMAu_draw)';
SIGMAtrOMEGA = SIGMAu_chol*OMEGAstar';
PHIpower = eye(NumberOfLagsTimesNvobs);
irfs = zeros (options_.irf,dataset_.info.nvobs*M_.exo_nbr);
tmp3 = PHIpower(1:dataset_.info.nvobs,1:dataset_.info.nvobs)*SIGMAtrOMEGA;
irfs = zeros (options_.irf,dataset_.vobs*M_.exo_nbr);
tmp3 = PHIpower(1:dataset_.vobs,1:dataset_.vobs)*SIGMAtrOMEGA;
irfs(1,:) = tmp3(:)';
for t = 2:options_.irf
PHIpower = Companion_matrix*PHIpower;
tmp3 = PHIpower(1:dataset_.info.nvobs,1:dataset_.info.nvobs)*SIGMAtrOMEGA;
tmp3 = PHIpower(1:dataset_.vobs,1:dataset_.vobs)*SIGMAtrOMEGA;
irfs(t,:) = tmp3(:)'+kron(ones(1,M_.exo_nbr),mu);
end
tmp_dsgevar = kron(ones(options_.irf,1),mu);
for j = 1:(dataset_.info.nvobs*M_.exo_nbr)
for j = 1:(dataset_.vobs*M_.exo_nbr)
if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10
tmp_dsgevar(:,j) = (irfs(:,j));
end
end
if IRUN < MAX_nirfs_dsgevar
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.info.nvobs,M_.exo_nbr);
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.vobs,M_.exo_nbr);
else
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.info.nvobs,M_.exo_nbr);
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,dataset_.vobs,M_.exo_nbr);
instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ...
int2str(NumberOfIRFfiles_dsgevar) '.mat stock_irf_bvardsge;'];,
eval(['save ' instr]);
......@@ -242,7 +242,7 @@ while fpar<B
end
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
IRUN =0;
stock_irf_dsgevar = zeros(options_.irf,dataset_.info.nvobs,M_.exo_nbr,MAX_nirfs_dsgevar);
stock_irf_dsgevar = zeros(options_.irf,dataset_.vobs,M_.exo_nbr,MAX_nirfs_dsgevar);
end
end
if irun == MAX_nirfs_dsge || irun == B || fpar == B
......@@ -279,20 +279,6 @@ while fpar<B
ifil2 = ifil2 + 1;
irun2 = 0;
end
% if isoctave
% if (whoiam==0),
% printf(['Posterior IRF %3.f%% done\r'],(fpar/B*100));
% end
% elseif ~whoiam,
% waitbar(fpar/B,h);
% end
% if whoiam,
% if ~isoctave
% fprintf('Done! \n');
% end
% waitbarString = [ 'Subdraw ' int2str(fpar) '/' int2str(B) ' done.'];
% fMessageStatus((fpar-fpar0)/(B-fpar0),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
% end
dyn_waitbar((fpar-fpar0)/(B-fpar0),h);
end
......@@ -310,9 +296,5 @@ end
myoutput.OutputFileName = [OutputFileName_dsge;
OutputFileName_param;
OutputFileName_bvardsge];
myoutput.nosaddle = nosaddle;
myoutput.nosaddle = nosaddle;
matlab/PosteriorIRF_core2.m
+ 11
− 12
View file @ f6ce8926
......@@ -5,7 +5,7 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab)
% Perform in parallel execution a portion of the PosteriorIRF.m code.
% See also the comment in random_walk_metropolis_hastings_core.m funtion.
%
% INPUTS
% INPUTS
% See the comment in random_walk_metropolis_hastings_core.m funtion.
%
% OUTPUTS
......@@ -13,8 +13,8 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab)
% Contained:
% OutputFileName (i.e. the figures without the file .txt).
%
% ALGORITHM
% Portion of PosteriorIRF.m function code. Specifically the last 'for' cycle.
% ALGORITHM
% Portion of PosteriorIRF.m function code. Specifically the last 'for' cycle.
%
% SPECIAL REQUIREMENTS.
% None.
......@@ -36,7 +36,7 @@ function myoutput=PosteriorIRF_core2(myinputs,fpar,npar,whoiam, ThisMatlab)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_ M_
global options_ M_
if nargin<4,
whoiam=0;
......@@ -87,7 +87,7 @@ OutputFileName={};
subplotnum = 0;
for i=fpar:npar,
figunumber = 0;
for j=1:nvar
if max(abs(MeanIRF(:,j,i))) > options_.impulse_responses.plot_threshold
subplotnum = subplotnum+1;
......@@ -96,7 +96,7 @@ for i=fpar:npar,
elseif subplotnum == 1 && ~options_.relative_irf
hh = dyn_figure(options_,'Name',['Orthogonalized shock to ' tit(i,:)]);
end
set(0,'CurrentFigure',hh)
subplot(nn,nn,subplotnum);
if ~MAX_nirfs_dsgevar
......@@ -108,12 +108,12 @@ for i=fpar:npar,
set(h2,'FaceColor',[1 1 1]);
set(h2,'BaseValue',min(HPDIRF(:,1,j,i)));
plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',3)
% plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
% plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
box on
axis tight
xlim([1 options_.irf]);
hold off
else
else
h1 = area(1:options_.irf,HPDIRF(:,2,j,i));
set(h1,'FaceColor',[.9 .9 .9]);
set(h1,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,1,j,i))]));
......@@ -136,9 +136,9 @@ for i=fpar:npar,
else
if options_.debug
fprintf('POSTERIOR_IRF: The IRF of %s to %s is smaller than the irf_plot_threshold of %4.3f and will not be displayed.\n',deblank(varlist(j,:)),tit(i,:),options_.impulse_responses.plot_threshold)
end
end
end
if subplotnum == MaxNumberOfPlotPerFigure || (j == nvar && subplotnum> 0)
figunumber = figunumber+1;
dyn_saveas(hh,[DirectoryName '/' M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)],options_);
......@@ -154,9 +154,8 @@ for i=fpar:npar,
% fMessageStatus((i-fpar+1)/(npar-fpar+1),whoiam,waitbarString, waitbarTitle, Parallel(ThisMatlab));
dyn_waitbar((i-fpar+1)/(npar-fpar+1),[],waitbarString);
end
end% loop over exo_var
end% loop over exo_var
myoutput.OutputFileName = OutputFileName;
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment