From 47ae59505945d9ac74eb39ac0a34d22b0bfbd510 Mon Sep 17 00:00:00 2001
From: adjemian <adjemian@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Fri, 15 Dec 2006 11:37:24 +0000
Subject: [PATCH] Bug corrections : TeX output, oo_ structure (posterior
 density), posterior IRF. Added : Posterior IRFs for BVAR-DSGE, conditional
 forecasts

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1119 ac1d8469-bf42-47a9-8791-bf33cf982152
---
 matlab/DsgeLikelihood.m                   |   2 +-
 matlab/GetPosteriorParametersStatistics.m | 332 ++++++----
 matlab/PlotPosteriorDistributions.m       |   4 +-
 matlab/PosteriorIRF.m                     | 747 ++++++++++++++--------
 matlab/ReshapeMatFiles.m                  |  26 +-
 matlab/dynare_estimation.m                |  20 +-
 matlab/imcforecast.m                      | 119 ++++
 matlab/mcforecast3.m                      |  13 +
 matlab/plot_icforecast.m                  |  32 +
 9 files changed, 865 insertions(+), 430 deletions(-)
 create mode 100755 matlab/imcforecast.m
 create mode 100755 matlab/mcforecast3.m
 create mode 100755 matlab/plot_icforecast.m

diff --git a/matlab/DsgeLikelihood.m b/matlab/DsgeLikelihood.m
index 4ce33bcb91..cb242a8022 100644
--- a/matlab/DsgeLikelihood.m
+++ b/matlab/DsgeLikelihood.m
@@ -192,4 +192,4 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
   % Adds prior if necessary
   % ------------------------------------------------------------------------------
   lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4);
-  fval    = (likelihood-lnprior);
+  fval    = (likelihood-lnprior);
\ No newline at end of file
diff --git a/matlab/GetPosteriorParametersStatistics.m b/matlab/GetPosteriorParametersStatistics.m
index a947d3a12f..347b69b8ae 100644
--- a/matlab/GetPosteriorParametersStatistics.m
+++ b/matlab/GetPosteriorParametersStatistics.m
@@ -1,5 +1,22 @@
-function GetPosteriorParametersStatistics()
-% stephane.adjemian@ens.fr [09-09-2005]
+function get_posterior_parameters_statistics()
+% This function prints and saves posterior estimates after the mcmc
+% (+updates of oo_ & TeX output). 
+% 
+% INPUTS 
+%   None.
+%  
+% OUTPUTS 
+%   None.  
+%
+% ALGORITHM 
+%   None.       
+%
+% SPECIAL REQUIREMENTS
+%   None.
+%  
+%  
+% part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
+% Gnu Public License.
 global estim_params_ M_ options_ bayestopt_ oo_
 
 TeX   	= options_.TeX;
@@ -12,6 +29,8 @@ np    	= estim_params_.np ;
 nx    	= nvx+nvn+ncx+ncn+np;
 
 DirectoryName = CheckPath('metropolis');
+OutputDirectoryName = CheckPath('Output');
+
 load([ DirectoryName '/'  M_.fname '_mh_history'])
 FirstMhFile = record.KeepedDraws.FirstMhFile;
 FirstLine = record.KeepedDraws.FirstLine; ifil = FirstLine;
@@ -21,138 +40,199 @@ FirstMhFile = record.KeepedDraws.FirstMhFile;
 NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
 clear record;
 
-disp(' ')
-disp(' ')
-disp('ESTIMATION RESULTS')
-disp(' ')
-disp(sprintf('Log data density is %f.',oo_.MarginalDensity.ModifiedHarmonicMean))
 pnames=['     ';'beta ';'gamm ';'norm ';'invg ';'unif ';'invg2'];
 tit2 = sprintf('%10s %7s %10s %14s %4s %6s\n',' ','prior mean','post. mean','conf. interval','prior','pstdev');
+pformat = '%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f';
+
+disp(' ');disp(' ');disp('ESTIMATION RESULTS');disp(' ');
+disp(sprintf('Log data density is %f.',oo_.MarginalDensity.ModifiedHarmonicMean))
 if np
-  disp(' ')
-  disp('parameters')
-  disp(tit2)
-  ip = nvx+nvn+ncx+ncn+1;
-  for i=1:np
-    Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
-    [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
-    name = bayestopt_.name{ip};
-    disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
-		 name, ...
-		 bayestopt_.pmean(ip),post_mean,hpd_interval, ...
-		 pnames(bayestopt_.pshape(ip)+1,:), ...
-		 bayestopt_.pstdev(ip)));
-    eval(['oo_.posterior_mean.parameters.' name ' = post_mean;']);
-    eval(['oo_.posterior_hpdinf.parameters.' name ' = hpd_interval(1);']); 
-    eval(['oo_.posterior_hpdsup.parameters.' name ' = hpd_interval(2);']);
-    eval(['oo_.posterior_median.' name ' = post_median;']);
-    eval(['oo_.posterior_variance.' name ' = post_var;']);
-    eval(['oo_.posterior_deciles.' name ' = post_deciles;']);
-    eval(['oo_.posterior_density.' name ' = density;']);    
-    ip = ip+1;
-  end
+    if TeX
+        fid = TeXBegin(OutputDirectoryName,M_.fname,1,'parameters');
+    end
+    disp(' ')
+    disp('parameters')
+    disp(tit2)
+    ip = nvx+nvn+ncx+ncn+1;
+    for i=1:np
+        Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
+        [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
+        name = bayestopt_.name{ip};
+        disp(sprintf(pformat,name,bayestopt_.pmean(ip),post_mean,hpd_interval, ...
+		 pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
+        oo_ = Filloo(oo_,name,'parameters',post_mean,hpd_interval,post_median,post_var,post_deciles,density);    
+        if TeX
+            TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
+                    bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);
+        end
+        ip = ip+1;
+    end
+    if TeX
+        TeXEnd(fid,1,'parameters');
+    end
 end
 if nvx
-  ip = 1;
-  disp(' ')
-  disp('standard deviation of shocks')
-  disp(tit2)
-  for i=1:nvx
-    Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
-    [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
-    k = estim_params_.var_exo(i,1);
-    name = deblank(M_.exo_names(k,:));
-    disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
-		 name,bayestopt_.pmean(ip),post_mean, ...
-		 hpd_interval,pnames(bayestopt_.pshape(ip)+1,:), ...
-		 bayestopt_.pstdev(ip))); 
-    M_.Sigma_e(k,k) = post_mean*post_mean;
-    eval(['oo_.posterior_mean.shocks_std.' name ' = post_mean;']);
-    eval(['oo_.posterior_hpdinf.shocks_std.' name ' = hpd_interval(1);']); 
-    eval(['oo_.posterior_hpdsup.shocks_std.' name ' = hpd_interval(2);']);
-    eval(['oo_.posterior_median.shocks_std.' name ' = post_median;']);
-    eval(['oo_.posterior_variance.shocks_std.' name ' = post_var;']);
-    eval(['oo_.posterior_deciles.shocks_std.' name ' = post_deciles;']);
-    eval(['oo_.posterior_density.shocks_std.' name ' = density;']);
-    ip = ip+1;
-  end
+    if TeX
+        fid = TeXBegin(OutputDirectoryName,M_.fname,2,'standard deviation of structural shocks');
+    end
+    ip = 1;
+    disp(' ')
+    disp('standard deviation of shocks')
+    disp(tit2)
+    for i=1:nvx
+        Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
+        [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
+        k = estim_params_.var_exo(i,1);
+        name = deblank(M_.exo_names(k,:));
+        disp(sprintf(pformat,name,bayestopt_.pmean(ip),post_mean,hpd_interval,...
+                     pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip))); 
+        M_.Sigma_e(k,k) = post_mean*post_mean;
+        oo_ = Filloo(oo_,name,'shocks_std',post_mean,hpd_interval,post_median,post_var,post_deciles,density);
+        if TeX
+            TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
+                    bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);
+        end
+        ip = ip+1;
+    end
+    if TeX
+        TeXEnd(fid,2,'standard deviation of structural shocks');        
+    end
 end
 if nvn
-  disp(' ')
-  disp('standard deviation of measurement errors')
-  disp(tit2)
-  ip = nvx+1;
-  for i=1:nvn
-    Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
-    [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);    
-    name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
-    disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', ...
-		 name,...
-		 bayestopt_.pmean(ip), ...
-		 post_mean,hpd_interval, ...
-		 pnames(bayestopt_.pshape(ip)+1,:), ...
-		 bayestopt_.pstdev(ip)));
-    eval(['oo_.posterior_mean.measurement_errors_std.' name  ' = post_mean;']);
-    eval(['oo_.posterior_hpdinf.measurement_errors_std.' name ' = hpd_interval(1);']); 
-    eval(['oo_.posterior_hpdsup.measurement_errors_std.' name ' = hpd_interval(2);']);		      
-    eval(['oo_.posterior_median.measurement_errors_std.' name ' = post_median;']);
-    eval(['oo_.posterior_variance.measurement_errors_std.' name ' = post_var;']);
-    eval(['oo_.posterior_deciles.measurement_errors_std.' name ' = post_deciles;']);
-    eval(['oo_.posterior_density.measurement_errors_std.' name ' = density;']);    
-    ip = ip+1;
-  end
+    if TeX
+        fid = TeXBegin(OutputDirectoryName,M_.fname,3,'standard deviation of measurement errors')
+    end
+    disp(' ')
+    disp('standard deviation of measurement errors')
+    disp(tit2)
+    ip = nvx+1;
+    for i=1:nvn
+        Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
+        [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);    
+        name = deblank(options_.varobs(estim_params_.var_endo(i,1),:));
+        disp(sprintf(pformat,name,bayestopt_.pmean(ip),post_mean,hpd_interval, ...
+                     pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
+        oo_ = Filloo(oo_,name,'measurement_errors_std',post_mean,hpd_interval,post_median,...
+                     post_var,post_deciles,density);
+        if TeX
+            TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
+                    bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);
+        end
+        ip = ip+1;
+    end
+    if TeX
+        TeXEnd(fid,3,'standard deviation of measurement errors');        
+    end
 end
 if ncx
-  disp(' ')
-  disp('correlation of shocks')
-  disp(tit2)
-  ip = nvx+nvn+1;
-  for i=1:ncx
-    Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
-    [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
-    k1 = estim_params_.corrx(i,1);
-    k2 = estim_params_.corrx(i,2);
-    name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
-    NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
-    disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', name, ...
-		 bayestopt_.pmean(ip),post_mean,hpd_interval, ...
-		 pnames(bayestopt_.pshape(ip)+1,:), ...
-		 bayestopt_.pstdev(ip)));
-    eval(['oo_.posterior_mean.shocks_corr.' NAME ' = post_mean;']);
-    eval(['oo_.posterior_hpdinf.shocks_corr.' NAME ' = hpd_interval(1);']); 
-    eval(['oo_.posterior_hpdsup.shocks_corr.' NAME ' = hpd_interval(2);']);
-    eval(['oo_.posterior_median.shocks_corr.' NAME ' = post_median;']);
-    eval(['oo_.posterior_variance.shocks_corr.' NAME ' = post_var;']);
-    eval(['oo_.posterior_deciles.shocks_corr.' NAME ' = post_deciles;']);
-    eval(['oo_.posterior_density.shocks_corr.' NAME ' = density;']);	  
-    M_.Sigma_e(k1,k2) = post_mean*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
-    M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
-    ip = ip+1;
-  end
+    if TeX
+        fid = TeXBegin(OutputDirectoryName,M_.fname,4,'correlation of structural shocks');
+    end
+    disp(' ')
+    disp('correlation of shocks')
+    disp(tit2)
+    ip = nvx+nvn+1;
+    for i=1:ncx
+        Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
+        [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
+        k1 = estim_params_.corrx(i,1);
+        k2 = estim_params_.corrx(i,2);
+        name = [deblank(M_.exo_names(k1,:)) ',' deblank(M_.exo_names(k2,:))];
+        NAME = [deblank(M_.exo_names(k1,:)) '_' deblank(M_.exo_names(k2,:))];
+        disp(sprintf(pformat, name,bayestopt_.pmean(ip),post_mean,hpd_interval, ...
+                     pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
+        oo_ = Filloo(oo_,NAME,'shocks_corr',post_mean,hpd_interval,post_median,post_var,post_deciles,density);	  
+        M_.Sigma_e(k1,k2) = post_mean*sqrt(M_.Sigma_e(k1,k1)*M_.Sigma_e(k2,k2));
+        M_.Sigma_e(k2,k1) = M_.Sigma_e(k1,k2);
+        if TeX
+            TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
+                    bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);
+        end
+        ip = ip+1;
+    end
+    if TeX
+        TeXEnd(fid,4,'correlation of structural shocks');
+    end
 end
 if ncn
-  disp(' ')
-  disp('correlation of measurement errors')
-  disp(tit2)
-  ip = nvx+nvn+ncx+1;
-  for i=1:ncn
-    Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
-    [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
-    k1 = estim_params_.corrn(i,1);
-    k2 = estim_params_.corrn(i,2);
-    name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];
-    NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
-    disp(sprintf('%12s %7.3f %8.4f %7.4f %7.4f %4s %6.4f', name, ...
-		 bayestopt_.pmean(ip),post_mean,hpd_interval, ...
-		 pnames(bayestopt_.pshape(ip)+1,:), ...
-		 bayestopt_.pstdev(ip))); 
-    eval(['oo_.posterior_mean.measurement_errors_corr.' NAME ' = post_mean;']);
-    eval(['oo_.posterior_hpdinf.measurement_errors_corr.' NAME ' = hpd_interval(1);']); 
-    eval(['oo_.posterior_hpdsup.measurement_errors_corr.' NAME ' = hpd_interval(2);']);      
-    eval(['oo_.posterior_median.measurement_errors_corr.' NAME ' = post_median;']);
-    eval(['oo_.posterior_variance.measurement_errors_corr.' NAME ' = post_var;']);
-    eval(['oo_.posterior_deciles.measurement_errors_corr.' NAME ' = post_decile;']);
-    eval(['oo_.posterior_density.measurement_errors_corr.' NAME ' = density;']);
-    ip = ip+1;
-  end
-end
\ No newline at end of file
+    if TeX
+        fid = TeXBegin(OutputDirectoryName,M_.fname,5,'correlation of measurement errors');
+    end
+    disp(' ')
+    disp('correlation of measurement errors')
+    disp(tit2)
+    ip = nvx+nvn+ncx+1;
+    for i=1:ncn
+        Draws = GetAllPosteriorDraws(ip,FirstMhFile,FirstLine,TotalNumberOfMhFiles,NumberOfDraws);
+        [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = posterior_moments(Draws,1);
+        k1 = estim_params_.corrn(i,1);
+        k2 = estim_params_.corrn(i,2);
+        name = [deblank(M_.endo_names(k1,:)) ',' deblank(M_.endo_names(k2,:))];
+        NAME = [deblank(M_.endo_names(k1,:)) '_' deblank(M_.endo_names(k2,:))];
+        disp(sprintf(pformat, name,bayestopt_.pmean(ip),post_mean,hpd_interval, ...
+                     pnames(bayestopt_.pshape(ip)+1,:),bayestopt_.pstdev(ip)));
+        oo_ = Filloo(oo_,NAME,'measurement_errors_corr',post_mean,hpd_interval,...
+                     post_median,post_var,post_deciles,density);
+        if TeX
+            TeXCore(fid,name,deblank(pnames(bayestopt_.pshape(ip)+1,:)),bayestopt_.pmean(ip),...
+                    bayestopt_.pstdev(ip),post_mean,sqrt(post_var),hpd_interval);            
+        end
+        ip = ip+1;
+    end
+    if TeX
+        TeXEnd(fid,5,'correlation of measurement errors');        
+    end
+end
+
+
+%
+%% subfunctions:
+%
+function fid = TeXBegin(directory,fname,fnum,title)
+    TeXfile = [directory '/' fname '_Posterior_Mean_' int2str(fnum) '.TeX'];
+    fidTeX = fopen(TeXfile,'w');
+    fprintf(fidTeX,'%% TeX-table generated by Dynare.\n');
+    fprintf(fidTeX,['%% RESULTS FROM METROPOLIS HASTINGS (' title ')\n']);
+    fprintf(fidTeX,['%% ' datestr(now,0)]);
+    fprintf(fidTeX,' \n');
+    fprintf(fidTeX,' \n');
+    fprintf(fidTeX,'\\begin{table}\n');
+    fprintf(fidTeX,'\\centering\n');
+    fprintf(fidTeX,'\\begin{tabular}{l|lcccccc} \n');
+    fprintf(fidTeX,'\\hline\\hline \\\\ \n');
+    fprintf(fidTeX,['  & Prior distribution & Prior mean  & Prior ' ...
+                    's.d. & Posterior mean & Posterior s.d.  & HPD inf & HPD sup\\\\ \n']);
+    fprintf(fidTeX,'\\hline \\\\ \n');
+    fid = fidTeX;
+
+    
+function TeXCore(fid,name,shape,priormean,priorstd,postmean,poststd,hpd)
+    fprintf(fid,['$%s$ & %s & %7.3f & %6.4f & %7.3f& %6.4f & %7.4f & %7.4f \\\\ \n'],... 
+            name,...
+            shape,...
+            priormean,...
+            priorstd,...
+            postmean,...
+            poststd,...
+            hpd(1),...
+            hpd(2));
+    
+
+function TeXEnd(fid,fnum,title)
+    fprintf(fid,'\\hline\\hline \n');
+    fprintf(fid,'\\end{tabular}\n ');    
+    fprintf(fid,['\\caption{Results from Metropolis-Hastings (' title ')}\n ']);
+    fprintf(fid,['\\label{Table:MHPosterior:' int2str(fnum)  '}\n']);
+    fprintf(fid,'\\end{table}\n');
+    fprintf(fid,'%% End of TeX file.\n');
+    fclose(fid);
+    
+    
+function oo = Filloo(oo,name,type,postmean,hpdinterval,postmedian,postvar,postdecile,density)
+    eval(['oo.posterior_mean.' type '.' name ' = postmean;']);
+    eval(['oo.posterior_hpdinf.' type '.' name ' = hpdinterval(1);']); 
+    eval(['oo.posterior_hpdsup.' type '.' name ' = hpdinterval(2);']);      
+    eval(['oo.posterior_median.' type '.' name ' = postmedian;']);
+    eval(['oo.posterior_variance.' type '.' name ' = postvar;']);
+    eval(['oo.posterior_deciles.' type '.' name ' = postdecile;']);
+    eval(['oo.posterior_density.' type '.' name ' = density;']);
\ No newline at end of file
diff --git a/matlab/PlotPosteriorDistributions.m b/matlab/PlotPosteriorDistributions.m
index aefa2a5d12..c4ed888a28 100644
--- a/matlab/PlotPosteriorDistributions.m
+++ b/matlab/PlotPosteriorDistributions.m
@@ -88,8 +88,8 @@ for i=1:npar
   else
     j = i - (nvx+nvn+ncx+ncn);
     name = deblank(M_.param_names(estim_params_.param_vals(j,1),:));
-    eval(['x1 = oo_.posterior_density.' name '(:,1);'])
-    eval(['f1 = oo_.posterior_density.' name '(:,2);'])
+    eval(['x1 = oo_.posterior_density.parameters.' name '(:,1);'])
+    eval(['f1 = oo_.posterior_density.parameters.' name '(:,2);'])
     if options_.posterior_mode_estimation
       eval(['pmode = oo_.posterior_mode.parameters.' name ';'])
     end
diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m
index 63f3bcc795..a0f7b762c1 100644
--- a/matlab/PosteriorIRF.m
+++ b/matlab/PosteriorIRF.m
@@ -1,281 +1,466 @@
-function posterior_irf(type)
-% Metropolis-Hastings algorithm. 
-% 
-% INPUTS 
-%   o type       [char]     string specifying the joint density of the
-%   deep parameters ('prior','posterior'). 
-%  
-% OUTPUTS 
-%   None (oo_ and plots).
-%
-%
-% ALGORITHM 
-%   None.       
-%
-% SPECIAL REQUIREMENTS
-%   None.
-%  
-%  
-% part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
-% Gnu Public License.
-global options_ estim_params_ oo_ M_ dsge_prior_weight
-nvx  = estim_params_.nvx;
-nvn  = estim_params_.nvn;
-ncx  = estim_params_.ncx;
-ncn  = estim_params_.ncn;
-np   = estim_params_.np ;
-npar = nvx+nvn+ncx+ncn+np;
-offset = npar-np;
-%
-MaxNumberOfPlotPerFigure = 9;% The square root must be an integer!
-nn = sqrt(MaxNumberOfPlotPerFigure);
-DirectoryName = CheckPath('Output');
-if strcmpi(type,'posterior')
-  MhDirectoryName = CheckPath('metropolis');
-elseif strcmpi(type,'gsa')
-  MhDirectoryName = CheckPath('GSA');
-else
-  MhDirectoryName = CheckPath('prior');
-end  
-MAX_nirfs = ceil(options_.MaxNumberOfBytes/(options_.irf*length(oo_.steady_state)*M_.exo_nbr)/8)+50;
-MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
-
-if strcmpi(type,'posterior')
-  load([ MhDirectoryName '/'  M_.fname '_mh_history'])
-  TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
-  NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
-elseif strcmpi(type,'gsa')
-  load([ MhDirectoryName '/'  M_.fname '_prior'],'lpmat','istable')
-  x=lpmat(istable,:);
-  clear lpmat istable
-  NumberOfDraws=size(x,1);
-  B=NumberOfDraws; options_.B = B;
-else% type = 'prior'
-  NumberOfDraws = 500;
-end
-if ~strcmpi(type,'gsa')
-  B = min([round(.5*NumberOfDraws),500]); options_.B = B;
-end
-try delete([MhDirectoryName '\' M_.fname '_IRFs*']);
-catch disp('No _IRFs files to be deleted!')
-end
-irun = 0;
-irun2 = 0;
-NumberOfIRFfiles = 1;
-ifil2 = 1;
-if strcmpi(type,'posterior')
-  h = waitbar(0,'Bayesian (posterior) IRFs...');
-elseif strcmpi(type,'gsa')
-  h = waitbar(0,'GSA (prior) IRFs...');
-else
-  h = waitbar(0,'Bayesian (prior) IRFs...');
-end
-if B <= MAX_nruns
-  stock_param = zeros(B, npar);
-else
-  stock_param = zeros(MAX_nruns, npar);
-end
-if B >= MAX_nirfs
-  stock_irf = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,MAX_nirfs);
-else
-  stock_irf = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,B);
-end
-for b=1:B
-  irun = irun+1;
-  irun2 = irun2+1;
-  if ~strcmpi(type,'gsa')
-    deep = GetOneDraw(type);
-  else
-    deep = x(b,:);
-  end
-  stock_param(irun2,:) = deep;  
-  set_parameters(deep);
-  dr = resol(oo_.steady_state,0);
-  SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr);
-  cs = transpose(chol(SS));
-  for i = 1:M_.exo_nbr
-    if SS(i,i) > 1e-13
-      y=irf(dr,SS(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
-      for j = 1:M_.endo_nbr
-        if max(y(j,:)) - min(y(j,:)) > 1e-10 
-          stock_irf(:,j,i,irun) = transpose(y(j,:));
-        end
-      end
-    end
-  end
-  if irun == MAX_nirfs | irun == B | b == B
-    if b == B
-      stock_irf = stock_irf(:,:,:,1:irun);
-    end
-    save([MhDirectoryName '/' M_.fname '_irf' int2str(NumberOfIRFfiles)],'stock_irf');
-    NumberOfIRFfiles = NumberOfIRFfiles+1;
-    irun = 0;
-  end
-  if irun2 == MAX_nruns | b == B
-    if b == B
-      stock_param = stock_param(1:irun2,:);
-    end
-    stock = stock_param;
-    save([MhDirectoryName '/' M_.fname '_param_irf' int2str(ifil2)],'stock');
-    ifil2 = ifil2 + 1;
-    irun2 = 0;
-  end
-  waitbar(b/B,h);
-end
-NumberOfIRFfiles = NumberOfIRFfiles-1;
-ifil2 = ifil2-1;
-close(h);
-
-ReshapeMatFiles('irf',type)
-
-if strcmpi(type,'gsa')
-  return
-end  
-
-
-varlist = options_.varlist;
-if isempty(varlist)
-  varlist = M_.endo_names;
-  SelecVariables = transpose(1:M_.endo_nbr);
-  nvar = M_.endo_nbr;
-else
-  nvar = size(varlist,1);
-  SelecVariables = [];
-  for i=1:nvar
-    if ~isempty(strmatch(varlist(i,:),M_.endo_names,'exact'))
-      SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
-    end
-  end
-end
-
-MeanIRF = zeros(options_.irf,nvar,M_.exo_nbr);
-MedianIRF = zeros(options_.irf,nvar,M_.exo_nbr);
-StdIRF = zeros(options_.irf,nvar,M_.exo_nbr);
-DistribIRF = zeros(options_.irf,9,nvar,M_.exo_nbr);
-HPDIRF = zeros(options_.irf,2,nvar,M_.exo_nbr);
-
-if options_.TeX
-  varlist_TeX = [];
-  for i=1:nvar
-    varlist_TeX = strvcat(varlist_TeX,M_.endo_names_tex(SelecVariables(i),:));
-  end
-end
-
-fprintf('MH: Posterior IRFs...\n');
-tit(M_.exo_names_orig_ord,:) = M_.exo_names;
-kdx = 0;
-for file = 1:NumberOfIRFfiles
-  load([MhDirectoryName '/' M_.fname '_IRFs' int2str(file)]);
-  for i = 1:M_.exo_nbr
-    for j = 1:nvar
-      for k = 1:size(STOCK_IRF,1)
-        kk = k+kdx;
-        [MeanIRF(kk,j,i),MedianIRF(kk,j,i),VarIRF(kk,j,i),HPDIRF(kk,:,j,i),DistribIRF(kk,:,j,i)] = ...
-          posterior_moments(squeeze(STOCK_IRF(k,SelecVariables(j),i,:)),0);
-      end
-    end
-  end
-  kdx = kdx + size(STOCK_IRF,1);
-end
-clear STOCK_IRF;
-
-for i = 1:M_.exo_nbr
-  for j = 1:nvar
-    name = [deblank(M_.endo_names(SelecVariables(j),:)) '_' deblank(tit(i,:))];
-    eval(['oo_.PosteriorIRF.Mean.' name ' = MeanIRF(:,j,i);']);
-    eval(['oo_.PosteriorIRF.Median.' name ' = MedianIRF(:,j,i);']);
-    eval(['oo_.PosteriorIRF.Var.' name ' = VarIRF(:,j,i);']);
-    eval(['oo_.PosteriorIRF.Distribution.' name ' = DistribIRF(:,:,j,i);']);
-    eval(['oo_.PosteriorIRF.HPDinf.' name ' = HPDIRF(:,1,j,i);']);
-    eval(['oo_.PosteriorIRF.HPDsup.' name ' = HPDIRF(:,2,j,i);']);
-  end
-end
-%%
-%% 	Finally i build the plots.
-%%
-if options_.TeX
-  fidTeX = fopen([DirectoryName '/' M_.fname '_BayesianIRF.TeX'],'w');
-  fprintf(fidTeX,'%% TeX eps-loader file generated by PosteriorIRF.m (Dynare).\n');
-  fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
-  fprintf(fidTeX,' \n');
-  titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex;
-end
-%%
-subplotnum = 0;
-for i=1:M_.exo_nbr
-  NAMES = [];
-  if options_.TeX
-    TEXNAMES = [];
-  end
-  figunumber = 0;
-  for j=1:nvar
-    if max(abs(MeanIRF(:,j,i))) > 10^(-6)
-      subplotnum = subplotnum+1;
-      if options_.nograph
-        if subplotnum == 1 & options_.relative_irf
-          hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)],'Visible','off');
-        elseif subplotnum == 1 & ~options_.relative_irf
-          hh = figure('Name',['Orthogonalized shock to ' tit(i,:)],'Visible','off');
-        end
-      else
-        if subplotnum == 1 & options_.relative_irf
-          hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)]);
-        elseif subplotnum == 1 & ~options_.relative_irf
-          hh = figure('Name',['Orthogonalized shock to ' tit(i,:)]);
-        end
-      end
-      set(0,'CurrentFigure',hh)
-      subplot(nn,nn,subplotnum);
-      plot([1 options_.irf],[0 0],'-r','linewidth',0.5);
-      hold on
-      for k = 1:9
-        plot(1:options_.irf,DistribIRF(:,k,j,i),'-g','linewidth',0.5)
-      end
-      plot(1:options_.irf,MeanIRF(:,j,i),'-k','linewidth',1)
-      xlim([1 options_.irf]);
-      hold off
-      name = deblank(varlist(j,:));
-      NAMES = strvcat(NAMES,name);
-      if options_.TeX
-        texname = deblank(varlist_TeX(j,:));
-        TEXNAMES = strvcat(TEXNAMES,['$' texname '$']);
-      end
-      title(name,'Interpreter','none')
-    end
-    if subplotnum == MaxNumberOfPlotPerFigure | (j == nvar  & subplotnum>0)
-      figunumber = figunumber+1;
-      set(hh,'visible','on')
-      eval(['print -depsc2 ' DirectoryName '/'  M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)]);
-      eval(['print -dpdf ' DirectoryName '/' M_.fname  '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)]);
-      saveas(hh,[DirectoryName '/' M_.fname  '_Bayesian_IRF_' deblank(tit(i,:))  '_' int2str(figunumber) '.fig']);
-      set(hh,'visible','off')
-      if options_.nograph, close(hh), end
-      if options_.TeX
-        fprintf(fidTeX,'\\begin{figure}[H]\n');
-        for jj = 1:size(TEXNAMES,1)
-          fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
-        end    
-        fprintf(fidTeX,'\\centering \n');
-        fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_Bayesian_IRF_%s}\n',M_.fname,deblank(tit(i,:)));
-        if options_.relative_irf
-          fprintf(fidTeX,['\\caption{Bayesian relative IRF.}']);
-        else
-          fprintf(fidTeX,'\\caption{Bayesian IRF.}');
-        end
-        fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s}\n',deblank(tit(i,:)));
-        fprintf(fidTeX,'\\end{figure}\n');
-        fprintf(fidTeX,' \n');
-      end
-      subplotnum = 0;
-    end
-  end% loop over selected endo_var
-end% loop over exo_var  
-%%
-if options_.TeX
-  fprintf(fidTeX,'%% End of TeX file.\n');
-  fclose(fidTeX);
-end
-fprintf('MH: Posterior IRFs, done!\n');
+function posterior_irf(type)
+% Builds posterior IRFs after the MH algorithm. 
+% 
+% INPUTS 
+%   o type       [char]     string specifying the joint density of the
+%   deep parameters ('prior','posterior'). 
+%  
+% OUTPUTS 
+%   None (oo_ and plots).
+%
+%
+% ALGORITHM
+%   None.       
+%
+% SPECIAL REQUIREMENTS
+%   None.
+%  
+%  
+% part of DYNARE, copyright S. Adjemian, M. Juillard (2006)
+% Gnu Public License.
+global options_ estim_params_ oo_ M_ bayestopt_
+if isempty(options_.irf) | ~options_.irf 
+    options_.irf = 40;
+end
+nvx  = estim_params_.nvx;
+nvn  = estim_params_.nvn;
+ncx  = estim_params_.ncx;
+ncn  = estim_params_.ncn;
+np   = estim_params_.np ;
+npar = nvx+nvn+ncx+ncn+np;
+offset = npar-np;
+nvobs = size(options_.varobs,1);
+gend = options_.nobs;
+%
+MaxNumberOfPlotPerFigure = 9;% The square root must be an integer!
+nn = sqrt(MaxNumberOfPlotPerFigure);
+DirectoryName = CheckPath('Output');
+if strcmpi(type,'posterior')
+  MhDirectoryName = CheckPath('metropolis');
+elseif strcmpi(type,'gsa')
+  MhDirectoryName = CheckPath('GSA');
+else
+  MhDirectoryName = CheckPath('prior');
+end
+MAX_nirfs_dsge = ceil(options_.MaxNumberOfBytes/(options_.irf*length(oo_.steady_state)*M_.exo_nbr)/8)+50;
+MAX_nruns = ceil(options_.MaxNumberOfBytes/(npar+2)/8);
+if ~isempty(strmatch('dsge_prior_weight',M_.param_names))
+    MAX_nirfs_dsgevar = ceil(options_.MaxNumberOfBytes/(options_.irf*nvobs*M_.exo_nbr)/8)+50;
+else
+    MAX_nirfs_dsgevar = 0;
+end
+if strcmpi(type,'posterior')
+  load([ MhDirectoryName '/'  M_.fname '_mh_history'])
+  TotalNumberOfMhDraws = sum(record.MhDraws(:,1));
+  NumberOfDraws = TotalNumberOfMhDraws-floor(options_.mh_drop*TotalNumberOfMhDraws);
+elseif strcmpi(type,'gsa')
+  load([ MhDirectoryName '/'  M_.fname '_prior'],'lpmat','istable')
+  x=lpmat(istable,:);
+  clear lpmat istable
+  NumberOfDraws=size(x,1);
+  B=NumberOfDraws; options_.B = B;
+else% type = 'prior'
+  NumberOfDraws = 500;
+end
+if ~strcmpi(type,'gsa')
+  B = min([round(.5*NumberOfDraws),500]); options_.B = B;
+end
+try delete([MhDirectoryName '/' M_.fname '_irf_dsge*.mat'])
+catch disp('No _IRFs (dsge) files to be deleted!')
+end
+try delete([MhDirectoryName '/' M_.fname '_irf_bvardsge*.mat'])
+catch disp('No _IRFs (bvar-dsge) files to be deleted!')
+end
+irun = 0;
+IRUN = 0;
+irun2 = 0;
+NumberOfIRFfiles_dsge = 1;
+NumberOfIRFfiles_dsgevar = 1;
+ifil2 = 1;
+if strcmpi(type,'posterior')
+  h = waitbar(0,'Bayesian (posterior) IRFs...');
+elseif strcmpi(type,'gsa')
+  h = waitbar(0,'GSA (prior) IRFs...');
+else
+  h = waitbar(0,'Bayesian (prior) IRFs...');
+end
+if B <= MAX_nruns
+  stock_param = zeros(B, npar);
+else
+  stock_param = zeros(MAX_nruns, npar);
+end
+if B >= MAX_nirfs_dsge
+  stock_irf_dsge = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,MAX_nirfs_dsge);
+else
+  stock_irf_dsge = zeros(options_.irf,M_.endo_nbr,M_.exo_nbr,B);
+end
+if MAX_nirfs_dsgevar
+    if B >= MAX_nirfs_dsgevar
+        stock_irf_bvardsge = zeros(options_.irf,nvobs,M_.exo_nbr,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_.varlag,-1);
+    NumberOfLags = options_.varlag;
+    NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
+    COMP_draw = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
+end
+for b=1:B
+  irun = irun+1;
+  irun2 = irun2+1;
+  if ~strcmpi(type,'gsa')
+    deep = GetOneDraw(type);
+  else
+    deep = x(b,:);
+  end
+  stock_param(irun2,:) = deep;  
+  set_parameters(deep);
+  dr = resol(oo_.steady_state,0);
+  SS(M_.exo_names_orig_ord,M_.exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr);
+  SS = transpose(chol(SS));
+  for i = 1:M_.exo_nbr
+    if SS(i,i) > 1e-13
+      y=irf(dr,SS(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
+      for j = 1:M_.endo_nbr
+        if max(y(j,:)) - min(y(j,:)) > 1e-10 
+          stock_irf_dsge(:,j,i,irun) = transpose(y(j,:));
+        end
+      end
+    end
+  end
+  if MAX_nirfs_dsgevar
+      IRUN = IRUN+1;
+      tmp_dsgevar = zeros(options_.irf,nvobs*M_.exo_nbr);
+      [fval,cost_flag,ys,trend_coeff,info,PHI,SIGMAu,iXX] =  DsgeVarLikelihood(deep',gend);
+      dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
+      DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight));
+      tmp1 = SIGMAu*gend*(dsge_prior_weight+1);
+      val  = 1;
+      tmp1 = chol(inv(tmp1))'; 
+      while val;
+          % draw from the marginal posterior of sig
+          tmp2 = tmp1*randn(nvobs,DSGE_PRIOR_WEIGHT-NumberOfLagsTimesNvobs);
+          SIGMAu_draw = inv(tmp2*tmp2');
+          % draw from the conditional posterior of PHI
+          VARvecPHI = kron(SIGMAu_draw,iXX);
+          PHI_draw  = PHI(:) + chol(VARvecPHI)'*randn(nvobs*NumberOfLagsTimesNvobs,1);
+          COMP_draw(1:nvobs,:) = reshape(PHI_draw,NumberOfLagsTimesNvobs,nvobs)';
+          % Check for stationarity
+          tests = find(abs(eig(COMP_draw))>0.9999999999);
+          if isempty(tests)
+	      val=0;
+          end
+      end
+      % Get rotation
+      if dsge_prior_weight > 0
+          Atheta(oo_.dr.order_var,M_.exo_names_orig_ord) = oo_.dr.ghu*sqrt(M_.Sigma_e);
+          A0 = Atheta(bayestopt_.mfys,:);
+          [OMEGAstar,SIGMAtr] = qr2(A0');
+      end
+      SIGMAu_chol = chol(SIGMAu_draw)';
+      SIGMAtrOMEGA = SIGMAu_chol*OMEGAstar';
+      PHIpower = eye(NumberOfLagsTimesNvobs);
+      irfs = zeros (options_.irf,nvobs*M_.exo_nbr);
+      tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
+      irfs(1,:) = tmp3(:)';
+      for t = 2:options_.irf
+          PHIpower = COMP_draw*PHIpower;
+          tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
+          irfs(t,:)  = tmp3(:)';
+      end            
+      for j = 1:(nvobs*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,nvobs,M_.exo_nbr);
+      else
+          stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,nvobs,M_.exo_nbr);
+          instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) ' stock_irf_bvardsge;'];,
+          eval(['save ' instr]);
+          NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1; 
+          IRUN =0;
+          stock_irf_dsgevar = zeros(options_.irf,nvobs,M_.exo_nbr,MAX_nirfs_dsgevar);
+      end
+  end
+  if irun == MAX_nirfs_dsge | irun == B | b == B
+    if b == B
+      stock_irf_dsge = stock_irf_dsge(:,:,:,1:irun);
+      if MAX_nirfs_dsgevar & (b == B | IRUN == B)
+          stock_irf_bvardsge = stock_irf_bvardsge(:,:,:,1:IRUN);
+          instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) ' stock_irf_bvardsge;'];,
+          eval(['save ' instr]);
+          NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
+          irun = 0;
+      end
+    end
+    save([MhDirectoryName '/' M_.fname '_irf_dsge' int2str(NumberOfIRFfiles_dsge)],'stock_irf_dsge');
+    NumberOfIRFfiles_dsge = NumberOfIRFfiles_dsge+1;
+    irun = 0;
+  end
+  if irun2 == MAX_nruns | b == B
+    if b == B
+      stock_param = stock_param(1:irun2,:);
+    end
+    stock = stock_param;
+    save([MhDirectoryName '/' M_.fname '_param_irf' int2str(ifil2)],'stock');
+    ifil2 = ifil2 + 1;
+    irun2 = 0;
+  end
+  waitbar(b/B,h);
+end
+NumberOfIRFfiles_dsge = NumberOfIRFfiles_dsge-1;
+NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar-1;
+
+ifil2 = ifil2-1;
+close(h);
+
+ReshapeMatFiles('irf_dsge')
+if MAX_nirfs_dsgevar
+    ReshapeMatFiles('irf_bvardsge')
+end
+
+if strcmpi(type,'gsa')
+  return
+end  
+
+
+varlist = options_.varlist;
+
+if isempty(varlist)
+  if MAX_nirfs_dsgevar
+      varlist = options_.varobs;
+      nvar = size(varlist,1);
+      SelecVariables = [];
+      for i=1:nvar
+          if ~isempty(strmatch(varlist(i,:),M_.endo_names,'exact'))
+              SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
+          end
+      end
+  else
+      varlist = M_.endo_names;
+      SelecVariables = transpose(1:M_.endo_nbr);
+      nvar = M_.endo_nbr;
+  end
+else
+  nvar = size(varlist,1);
+  SelecVariables = [];
+  for i=1:nvar
+    if ~isempty(strmatch(varlist(i,:),M_.endo_names,'exact'))
+      SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
+    end
+  end
+  if MAX_nirfs_dsgevar% Here I test if declared variables are observed variables.
+      SelecVariables = [];
+      varlistbis = [];
+      for i=1:nvar
+          if ~isempty(strmatch(varlist(i,:),options_.varobs,'exact'))
+              SelecVariables = [SelecVariables;strmatch(varlist(i,:),M_.endo_names,'exact')];
+              varlistbis = strvcat(varlistbis,varlist(i,:));
+          else
+              disp(' ')
+              disp(['Warning :: ' varlist(i,:) 'is not an observed variable!'])
+              disp(['This variable will not be considered for the IRFs.'])
+              disp(' ')
+          end
+      end
+      varlist = varlistbis;
+  end 
+end
+
+
+MeanIRF = zeros(options_.irf,nvar,M_.exo_nbr);
+MedianIRF = zeros(options_.irf,nvar,M_.exo_nbr);
+VarIRF = zeros(options_.irf,nvar,M_.exo_nbr);
+DistribIRF = zeros(options_.irf,9,nvar,M_.exo_nbr);
+HPDIRF = zeros(options_.irf,2,nvar,M_.exo_nbr);
+
+if options_.TeX
+  varlist_TeX = [];
+  for i=1:nvar
+    varlist_TeX = strvcat(varlist_TeX,M_.endo_names_tex(SelecVariables(i),:));
+  end
+end
+
+fprintf('MH: Posterior (dsge) IRFs...\n');
+tit(M_.exo_names_orig_ord,:) = M_.exo_names;
+kdx = 0;
+for file = 1:NumberOfIRFfiles_dsge
+  load([MhDirectoryName '/' M_.fname '_IRF_DSGEs' int2str(file)]);
+  for i = 1:M_.exo_nbr
+    for j = 1:nvar
+      for k = 1:size(STOCK_IRF_DSGE,1)
+        kk = k+kdx;
+        [MeanIRF(kk,j,i),MedianIRF(kk,j,i),VarIRF(kk,j,i),HPDIRF(kk,:,j,i),DistribIRF(kk,:,j,i)] = ...
+          posterior_moments(squeeze(STOCK_IRF_DSGE(k,SelecVariables(j),i,:)),0);
+      end
+    end
+  end
+  kdx = kdx + size(STOCK_IRF_DSGE,1);
+end
+clear STOCK_IRF_DSGE;
+
+for i = 1:M_.exo_nbr
+  for j = 1:nvar
+    name = [deblank(M_.endo_names(SelecVariables(j),:)) '_' deblank(tit(i,:))];
+    eval(['oo_.PosteriorIRF.dsge.Mean.' name ' = MeanIRF(:,j,i);']);
+    eval(['oo_.PosteriorIRF.dsge.Median.' name ' = MedianIRF(:,j,i);']);
+    eval(['oo_.PosteriorIRF.dsge.Var.' name ' = VarIRF(:,j,i);']);
+    eval(['oo_.PosteriorIRF.dsge.Distribution.' name ' = DistribIRF(:,:,j,i);']);
+    eval(['oo_.PosteriorIRF.dsge.HPDinf.' name ' = HPDIRF(:,1,j,i);']);
+    eval(['oo_.PosteriorIRF.dsge.HPDsup.' name ' = HPDIRF(:,2,j,i);']);
+  end
+end
+
+
+if MAX_nirfs_dsgevar
+    MeanIRFdsgevar = zeros(options_.irf,nvar,M_.exo_nbr);
+    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);    
+    fprintf('MH: Posterior (bvar-dsge) IRFs...\n');
+    tit(M_.exo_names_orig_ord,:) = M_.exo_names;
+    kdx = 0;
+    for file = 1:NumberOfIRFfiles_dsgevar
+        load([MhDirectoryName '/' M_.fname '_IRF_BVARDSGEs' int2str(file)]);
+        for i = 1:M_.exo_nbr
+            for j = 1:nvar
+                for k = 1:size(STOCK_IRF_BVARDSGE,1)
+                    kk = k+kdx;
+                    [MeanIRFdsgevar(kk,j,i),MedianIRFdsgevar(kk,j,i),VarIRFdsgevar(kk,j,i),...
+                     HPDIRFdsgevar(kk,:,j,i),DistribIRFdsgevar(kk,:,j,i)] = ...
+                        posterior_moments(squeeze(STOCK_IRF_BVARDSGE(k,j,i,:)),0);%SelecVariables(j)
+                end
+            end
+        end
+        kdx = kdx + size(STOCK_IRF_BVARDSGE,1);
+    end
+    clear STOCK_IRF_BVARDSGE; 
+    for i = 1:M_.exo_nbr
+        for j = 1:nvar
+            name = [deblank(M_.endo_names(SelecVariables(j),:)) '_' deblank(tit(i,:))];
+            eval(['oo_.PosteriorIRF.bvardsge.Mean.' name ' = MeanIRFdsgevar(:,j,i);']);
+            eval(['oo_.PosteriorIRF.bvardsge.Median.' name ' = MedianIRFdsgevar(:,j,i);']);
+            eval(['oo_.PosteriorIRF.bvardsge.Var.' name ' = VarIRFdsgevar(:,j,i);']);
+            eval(['oo_.PosteriorIRF.bvardsge.Distribution.' name ' = DistribIRFdsgevar(:,:,j,i);']);
+            eval(['oo_.PosteriorIRF.bvardsge.HPDinf.' name ' = HPDIRFdsgevar(:,1,j,i);']);
+            eval(['oo_.PosteriorIRF.bvardsge.HPDsup.' name ' = HPDIRFdsgevar(:,2,j,i);']);
+        end
+    end
+end
+%%
+%% 	Finally I build the plots.
+%%
+if options_.TeX
+  fidTeX = fopen([DirectoryName '/' M_.fname '_BayesianIRF.TeX'],'w');
+  fprintf(fidTeX,'%% TeX eps-loader file generated by PosteriorIRF.m (Dynare).\n');
+  fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
+  fprintf(fidTeX,' \n');
+  titTeX(M_.exo_names_orig_ord,:) = M_.exo_names_tex;
+end
+%%
+subplotnum = 0;
+for i=1:M_.exo_nbr
+  NAMES = [];
+  if options_.TeX
+    TEXNAMES = [];
+  end
+  figunumber = 0;
+  for j=1:nvar
+    if max(abs(MeanIRF(:,j,i))) > 10^(-6)
+      subplotnum = subplotnum+1;
+      if options_.nograph
+        if subplotnum == 1 & options_.relative_irf
+          hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)],'Visible','off');
+        elseif subplotnum == 1 & ~options_.relative_irf
+          hh = figure('Name',['Orthogonalized shock to ' tit(i,:)],'Visible','off');
+        end
+      else
+        if subplotnum == 1 & options_.relative_irf
+          hh = figure('Name',['Relative response to orthogonalized shock to ' tit(i,:)]);
+        elseif subplotnum == 1 & ~options_.relative_irf
+          hh = figure('Name',['Orthogonalized shock to ' tit(i,:)]);
+        end
+      end
+      set(0,'CurrentFigure',hh)
+      subplot(nn,nn,subplotnum);
+      if ~MAX_nirfs_dsgevar
+          h1 = area(1:options_.irf,HPDIRF(:,2,j,i));
+          set(h1,'FaceColor',[.9 .9 .9]);
+          set(h1,'BaseValue',min(HPDIRF(:,1,j,i)));
+          hold on
+          h2 = area(1:options_.irf,HPDIRF(:,1,j,i),'FaceColor',[1 1 1],'BaseValue',min(HPDIRF(:,1,j,i)));
+          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);          
+          box on
+          axis tight
+          xlim([1 options_.irf]);
+          hold off
+      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))]));
+          hold on;
+          h2 = area(1:options_.irf,HPDIRF(:,1,j,i));
+          set(h2,'FaceColor',[1 1 1]);
+          set(h2,'BaseValue',min([min(HPDIRF(:,1,j,i)),min(HPDIRFdsgevar(:,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,MeanIRFdsgevar(:,j,i),'--k','linewidth',2)
+          plot(1:options_.irf,HPDIRFdsgevar(:,1,j,i),'--k','linewidth',1)
+          plot(1:options_.irf,HPDIRFdsgevar(:,2,j,i),'--k','linewidth',1)
+          box on
+          axis tight
+          xlim([1 options_.irf]);
+          hold off
+      end
+      name = deblank(varlist(j,:));
+      NAMES = strvcat(NAMES,name);
+      if options_.TeX
+        texname = deblank(varlist_TeX(j,:));
+        TEXNAMES = strvcat(TEXNAMES,['$' texname '$']);
+      end
+      title(name,'Interpreter','none')
+    end
+    if subplotnum == MaxNumberOfPlotPerFigure | (j == nvar  & subplotnum> 0)
+      figunumber = figunumber+1;
+      set(hh,'visible','on')
+      eval(['print -depsc2 ' DirectoryName '/'  M_.fname '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)]);
+      eval(['print -dpdf ' DirectoryName '/' M_.fname  '_Bayesian_IRF_' deblank(tit(i,:)) '_' int2str(figunumber)]);
+      saveas(hh,[DirectoryName '/' M_.fname  '_Bayesian_IRF_' deblank(tit(i,:))  '_' int2str(figunumber) '.fig']);
+      set(hh,'visible','off')
+      if options_.nograph, close(hh), end
+      if options_.TeX
+        fprintf(fidTeX,'\\begin{figure}[H]\n');
+        for jj = 1:size(TEXNAMES,1)
+          fprintf(fidTeX,['\\psfrag{%s}[1][][0.5][0]{%s}\n'],deblank(NAMES(jj,:)),deblank(TEXNAMES(jj,:)));
+        end    
+        fprintf(fidTeX,'\\centering \n');
+        fprintf(fidTeX,'\\includegraphics[scale=0.5]{%s_Bayesian_IRF_%s}\n',M_.fname,deblank(tit(i,:)));
+        if options_.relative_irf
+          fprintf(fidTeX,['\\caption{Bayesian relative IRF.}']);
+        else
+          fprintf(fidTeX,'\\caption{Bayesian IRF.}');
+        end
+        fprintf(fidTeX,'\\label{Fig:BayesianIRF:%s}\n',deblank(tit(i,:)));
+        fprintf(fidTeX,'\\end{figure}\n');
+        fprintf(fidTeX,' \n');
+      end
+      subplotnum = 0;
+    end
+  end% loop over selected endo_var
+end% loop over exo_var  
+%%
+if options_.TeX
+  fprintf(fidTeX,'%% End of TeX file.\n');
+  fclose(fidTeX);
+end
+fprintf('MH: Posterior IRFs, done!\n');
\ No newline at end of file
diff --git a/matlab/ReshapeMatFiles.m b/matlab/ReshapeMatFiles.m
index 7f8b9d4496..133ecd487b 100644
--- a/matlab/ReshapeMatFiles.m
+++ b/matlab/ReshapeMatFiles.m
@@ -26,10 +26,14 @@ else
 end  
 end
   switch type
-   case 'irf'
-    CAPtype  = 'IRF';
+   case 'irf_dsge'
+    CAPtype  = 'IRF_DSGE';
     TYPEsize = [ options_.irf , M_.endo_nbr , M_.exo_nbr ];
-    TYPEarray = 4; 
+    TYPEarray = 4;    
+   case 'irf_bvardsge'
+    CAPtype  = 'IRF_BVARDSGE';
+    TYPEsize = [ options_.irf , size(options_.varobs,1) , M_.exo_nbr ];
+    TYPEarray = 4;      
    case 'smooth'
     CAPtype  = 'SMOOTH';
     TYPEsize = [ M_.endo_nbr , options_.nobs ];
@@ -59,7 +63,7 @@ end
     return
   end
   
-  TYPEfiles = dir([MhDirectoryName M_.fname '_' type '*']);
+  TYPEfiles = dir([MhDirectoryName M_.fname '_' type '*.mat']);
   NumberOfTYPEfiles = length(TYPEfiles);
   B = options_.B;
   
@@ -98,9 +102,10 @@ end
       %eval(['STOCK_' CAPtype ' = sort(stock_' type ',4);'])
       save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(1)],['STOCK_' CAPtype ]);
     end
-    for file = 1:NumberOfTYPEfiles
-      delete([MhDirectoryName M_.fname '_' type int2str(file) '.mat'])
-    end
+    % Original file format may be useful in some cases...
+    % for file = 1:NumberOfTYPEfiles
+    %  delete([MhDirectoryName M_.fname '_' type int2str(file) '.mat'])
+    % end
    case 3
     if NumberOfTYPEfiles>1
       NumberOfPeriodsPerTYPEfiles = ceil( TYPEsize(2)/NumberOfTYPEfiles );
@@ -132,7 +137,8 @@ end
       %eval(['STOCK_' CAPtype ' = sort(stock_' type ',3);'])
       save([MhDirectoryName M_.fname '_' CAPtype 's' int2str(1)],['STOCK_' CAPtype ]);      
     end
-    for file = 1:NumberOfTYPEfiles
-      delete([MhDirectoryName M_.fname '_' type  int2str(file) '.mat'])
-    end
+    % Original file format may be useful in some cases...
+    % for file = 1:NumberOfTYPEfiles
+    %   delete([MhDirectoryName M_.fname '_' type  int2str(file) '.mat'])
+    % end
   end
\ No newline at end of file
diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m
index f26bdea652..314940f2c5 100644
--- a/matlab/dynare_estimation.m
+++ b/matlab/dynare_estimation.m
@@ -650,7 +650,7 @@ OutputDirectoryName = CheckPath('Output');
 
 if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior mode) Latex output
   if np
-    filename = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_1.TeX'];
+    filename = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_1.TeX'];
     fidTeX = fopen(filename,'w');
     fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
     fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (parameters)\n');
@@ -670,7 +670,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
 	      M_.param_names_tex(estim_params_.param_vals(i,1),:),...
 	      deblank(pnames(bayestopt_.pshape(ip)+1,:)),...
 	      bayestopt_.pmean(ip),...
-	      estim_params_.param_vals(i,6),...
+	      bayestopt_.pstdev(ip),...
 	      xparam1(ip),...
 	      stdh(ip));
       ip = ip + 1;    
@@ -685,7 +685,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
     fclose(fidTeX);
   end
   if nvx
-    TeXfile = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_2.TeX'];
+    TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_2.TeX'];
     fidTeX = fopen(TeXfile,'w');
     fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
     fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (standard deviation of structural shocks)\n');
@@ -706,7 +706,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
 	      deblank(M_.exo_names_tex(k,:)),...
 	      deblank(pnames(bayestopt_.pshape(ip)+1,:)),...
 	      bayestopt_.pmean(ip),...
-	      estim_params_.var_exo(i,7),...
+	      bayestopt_.pstdev(ip),...
 	      xparam1(ip), ...
 	      stdh(ip)); 
       ip = ip+1;
@@ -721,7 +721,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
     fclose(fidTeX);
   end
   if nvn
-    TeXfile = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_3.TeX'];
+    TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_3.TeX'];
     fidTeX  = fopen(TeXfile,'w');
     fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
     fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (standard deviation of measurement errors)\n');
@@ -741,7 +741,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
 	      deblank(M_.endo_names_tex(idx,:)), ...
 	      deblank(pnames(bayestopt_.pshape(ip)+1,:)), ...        
 	      bayestopt_.pmean(ip), ...
-	      estim_params_.var_endo(i,7),...        
+	      bayestopt_.pstdev(ip),...        
 	      xparam1(ip),...
 	      stdh(ip)); 
       ip = ip+1;
@@ -755,7 +755,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
     fclose(fidTeX);
   end
   if ncx
-    TeXfile = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_4.TeX'];
+    TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_4.TeX'];
     fidTeX = fopen(TeXfile,'w');
     fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
     fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (correlation of structural shocks)\n');
@@ -776,7 +776,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
 	      [deblank(M_.exo_names_tex(k1,:)) ',' deblank(M_.exo_names_tex(k2,:))], ...
 	      deblank(pnames(bayestopt_.pshape(ip)+1,:)), ...
 	      bayestopt_.pmean(ip), ...
-	      estim_params_.corrx(i,8), ...
+	      bayestopt_.pstdev(ip), ...
 	      xparam1(ip), ...
 	      stdh(ip));
       ip = ip+1;
@@ -790,7 +790,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
     fclose(fidTeX);
   end
   if ncn
-    TeXfile = [OutputDirectoryName '\' M_.fname '_Posterior_Mode_5.TeX'];
+    TeXfile = [OutputDirectoryName '/' M_.fname '_Posterior_Mode_5.TeX'];
     fidTeX = fopen(TeXfile,'w');
     fprintf(fidTeX,'%% TeX-table generated by dynare_estimation (Dynare).\n');
     fprintf(fidTeX,'%% RESULTS FROM POSTERIOR MAXIMIZATION (correlation of measurement errors)\n');
@@ -811,7 +811,7 @@ if any(bayestopt_.pshape > 0) & options_.TeX %% Bayesian estimation (posterior m
 	      [deblank(M_.endo_names_tex(k1,:)) ',' deblank(M_.endo_names_tex(k2,:))], ...
 	      pnames(bayestopt_.pshape(ip)+1,:), ...
 	      bayestopt_.pmean(ip), ...
-	      estim_params_.corrn(i,8), ...
+	      bayestopt_.pstdev(ip), ...
 	      xparam1(ip), ...
 	      stdh(ip));
       ip = ip+1;
diff --git a/matlab/imcforecast.m b/matlab/imcforecast.m
new file mode 100755
index 0000000000..655f6e4b5d
--- /dev/null
+++ b/matlab/imcforecast.m
@@ -0,0 +1,119 @@
+function icforecast(ptype,cV,cS,cL,H,mcValue,B,ci)
+% stephane.adjemian@ens.fr
+global options_ oo_ M_
+   
+xparam = get_posterior_parameters(ptype);
+gend   = options_.nobs;
+
+% Read and demean data 
+rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
+rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:);
+if options_.loglinear == 1 & ~options_.logdata
+  rawdata = log(rawdata);
+end
+if options_.prefilter == 1
+  bayestopt_.mean_varobs = mean(rawdata,1);
+  data = transpose(rawdata-ones(gend,1)*bayestopt_.mean_varobs);
+else
+  data = transpose(rawdata);
+end
+
+set_parameters(xparam);
+
+[atT,innov,measurement_error,filtered_state_vector,ys,trend_coeff] = DsgeSmoother(xparam,gend,data);
+
+InitState(:,1) = atT(:,end);
+[T,R,ys,info] = dynare_resolve;
+
+sQ = sqrt(M_.Sigma_e);
+
+NumberOfStates = length(InitState);
+FORCS1 = zeros(NumberOfStates,H+1,B);
+
+for b=1:B
+    FORCS1(:,1,b) = InitState;
+end
+
+EndoSize = size(M_.endo_names,1);
+ExoSize = size(M_.exo_names,1);
+
+n1 = size(cV,1);
+n2 = size(cS,1);
+
+if n1 ~= n2
+    disp('imcforecast :: Error!')
+    disp(['imcforecast :: The number of variables doesn''t match the number of shocks'])
+    return
+end
+
+idx = [];
+jdx = [];
+
+for i = 1:n1
+    idx = [idx ; oo_.dr.inv_order_var(strmatch(deblank(cV(i,:)),M_.endo_names,'exact'))];
+    jdx = [jdx ; strmatch(deblank(cS(i,:)),M_.exo_names,'exact')];
+end
+mv = zeros(n1,NumberOfStates);
+mu = zeros(ExoSize,n2);
+for i=1:n1
+    mv(i,idx(i)) = 1;
+    mu(jdx(i),i) = 1;
+end
+
+
+if (size(mcValue,2) == 1);
+        mcValue = mcValue*ones(1,cL);
+else
+    cL = size(mcValue,2);
+end
+
+randn('state',0);
+
+for b=1:B
+    shocks = sQ*randn(ExoSize,H);
+    shocks(jdx,:) = zeros(length(jdx),H);
+    FORCS1(:,:,b) = mcforecast3(cL,H,mcValue,shocks,FORCS1(:,:,b),T,R,mv, mu);
+end
+
+mFORCS1 = mean(FORCS1,3);
+
+tt = (1-ci)/2;
+t1 = round(B*tt);
+t2 = round(B*(1-tt));
+
+forecasts.controled_variables = cV;
+forecasts.instruments = cS;
+
+for i = 1:EndoSize
+    eval(['forecasts.cond.mean.' deblank(M_.endo_names(oo_.dr.order_var(i),:)) ' = mFORCS1(i,:)'';']);
+    tmp = sort(squeeze(FORCS1(i,:,:))');
+    eval(['forecasts.cond.ci.' deblank(M_.endo_names(oo_.dr.order_var(i),:)) ...
+          ' = [tmp(t1,:)'' ,tmp(t2,:)'' ]'';']);
+end
+
+clear FORCS1;
+
+FORCS2 = zeros(NumberOfStates,H+1,B);
+for b=1:B
+    FORCS2(:,1,b) = InitState;
+end
+
+randn('state',0);
+
+for b=1:B
+    shocks = sQ*randn(ExoSize,H);
+    shocks(jdx,:) = zeros(length(jdx),H);
+    FORCS2(:,:,b) = mcforecast3(0,H,mcValue,shocks,FORCS2(:,:,b),T,R,mv, mu);
+end
+
+mFORCS2 = mean(FORCS2,3);
+
+
+for i = 1:EndoSize
+    eval(['forecasts.uncond.mean.' deblank(M_.endo_names(oo_.dr.order_var(i),:)) ' = mFORCS2(i,:)'';']);
+    tmp = sort(squeeze(FORCS2(i,:,:))');
+    eval(['forecasts.uncond.ci.' deblank(M_.endo_names(oo_.dr.order_var(i),:)) ...
+          ' = [tmp(t1,:)'' ,tmp(t2,:)'' ]'';']);
+end
+
+save('conditional_forecasts.mat','forecasts');
\ No newline at end of file
diff --git a/matlab/mcforecast3.m b/matlab/mcforecast3.m
new file mode 100755
index 0000000000..149b46c1d0
--- /dev/null
+++ b/matlab/mcforecast3.m
@@ -0,0 +1,13 @@
+function forcs = mcforecast3(cL,H,mcValue,shocks,forcs,T,R,mv,mu)
+% stephane.adjemian@ens.fr [06-11-2006]
+
+if cL
+    e = zeros(size(mcValue,1),cL);
+    for t = 1:cL
+        e(:,t) = inv(mv*R*mu)*(mcValue(:,t)-mv*T*forcs(:,t)-mv*R*shocks(:,t));
+        forcs(:,t+1) = T*forcs(:,t)+R*(mu*e(:,t)+shocks(:,t));
+    end
+end
+for t = cL+1:H
+    forcs(:,t+1) = T*forcs(:,t)+R*shocks(:,t);
+end
\ No newline at end of file
diff --git a/matlab/plot_icforecast.m b/matlab/plot_icforecast.m
new file mode 100755
index 0000000000..40daa3417d
--- /dev/null
+++ b/matlab/plot_icforecast.m
@@ -0,0 +1,32 @@
+function plot_icforecast(Variable)
+% stephane.adjemian@ens.fr
+    
+load conditional_forecasts;
+
+
+eval(['ci1 = forecasts.cond.ci.' Variable ';'])
+eval(['m1 = forecasts.cond.mean.' Variable ';'])
+eval(['ci2 = forecasts.uncond.ci.' Variable ';'])
+eval(['m2 = forecasts.uncond.mean.' Variable ';'])
+
+
+
+H = length(m1);
+
+% area(1:H,ci1(2,:),'FaceColor',[.9 .9 .9],'BaseValue',min([min(ci1(1,:)),min(ci2(1,:))]))
+
+h1 = area(1:H,ci1(2,1:H))
+set(h1,'BaseValue',min([min(ci1(1,:)),min(ci2(1,:))]))
+set(h1,'FaceColor',[.9 .9 .9])
+
+hold on
+% area(1:H,ci1(1,:),'FaceColor',[1 1 1],'BaseValue',min([min(ci1(1,:)),min(ci2(1,:))]))
+h2 = area(1:H,ci1(1,1:H));
+set(h2,'BaseValue',min([min(ci1(1,:)),min(ci2(1,:))]))
+set(h2,'FaceColor',[1 1 1])
+plot(1:H,m1,'-k','linewidth',3)
+plot(1:H,m2,'--k','linewidth',3)
+plot(1:H,ci2(1,:),'--k','linewidth',1)
+plot(1:H,ci2(2,:),'--k','linewidth',1)
+axis tight
+hold off
\ No newline at end of file
-- 
GitLab