diff --git a/matlab/DsgeLikelihood.m b/matlab/DsgeLikelihood.m
index 8b1cc8a8f..831f0277b 100644
--- a/matlab/DsgeLikelihood.m
+++ b/matlab/DsgeLikelihood.m
@@ -1,5 +1,4 @@
 function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data)
-
 % function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data)
 % Evaluates the posterior kernel of a dsge model. 
 % 
@@ -22,9 +21,7 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
 %  
 % part of DYNARE, copyright Dynare Team (2004-2008)
 % Gnu Public License.
-
   global bayestopt_ estim_params_ options_ trend_coeff_ M_ oo_ xparam1_test
-
   fval		= [];
   ys		= [];
   trend_coeff	= [];
@@ -155,7 +152,7 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
   % 3. Initial condition of the Kalman filter
   %------------------------------------------------------------------------------
   if options_.lik_init == 1		% Kalman filter
-    Pstar = lyapunov_symm(T,R*Q*R');
+    Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium);
     Pinf	= [];
   elseif options_.lik_init == 2	% Old Diffuse Kalman filter
     Pstar = 10*eye(np);
@@ -165,7 +162,7 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
       Pstar = zeros(np,np);
       ivs = bayestopt_.restrict_var_list_stationary;
       R1 = R(ivs,:);
-      Pstar(ivs,ivs) = lyapunov_symm(T(ivs,ivs),R1*Q*R1');
+      Pstar(ivs,ivs) = lyapunov_symm(T(ivs,ivs),R1*Q*R1',options_.qz_criterium);
       %    Pinf  = bayestopt_.Pinf;
       % by M. Ratto
       RR=T(:,bayestopt_.restrict_var_list_nonstationary);
@@ -245,15 +242,20 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
         LIK = DiffuseLikelihoodH3(T,R,Q,H,Pinf,Pstar,data,trend,start);
       else
 	LIK = DiffuseLikelihoodH3corr(T,R,Q,H,Pinf,Pstar,data,trend,start);
-      end	
-    end	  
+      end
+    end 
   else
     if options_.kalman_algo == 1
        %nv = size(bayestopt_.Z,1);
        %LIK = kalman_filter(bayestopt_.Z,zeros(nv,nv),T,R,Q,data,zeros(size(T,1),1),Pstar,'u');
-      LIK = DiffuseLikelihood1(T,R,Q,Pinf,Pstar,data,trend,start);
-      % LIK = diffuse_likelihood1(T,R,Q,Pinf,Pstar,data-trend,start);
-      %if abs(LIK1-LIK)>0.0000000001
+       %tic
+       LIK = DiffuseLikelihood1(T,R,Q,Pinf,Pstar,data,trend,start);
+       %toc
+       %tic
+       %LIK1 = Diffuse_Likelihood1(T,R,Q,Pinf,Pstar,data,trend,start);
+       %toc
+       %LIK1-LIK
+       %if abs(LIK1-LIK)>0.0000000001
       %  disp(['LIK1 and LIK are not equal! ' num2str(abs(LIK1-LIK))])
       %end
       if isinf(LIK)
diff --git a/matlab/DsgeLikelihood_hh.m b/matlab/DsgeLikelihood_hh.m
index bbdebae7f..88d9addb3 100644
--- a/matlab/DsgeLikelihood_hh.m
+++ b/matlab/DsgeLikelihood_hh.m
@@ -160,7 +160,7 @@ function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend
   % 3. Initial condition of the Kalman filter
   %------------------------------------------------------------------------------
   if options_.lik_init == 1		% Kalman filter
-    Pstar = lyapunov_symm(T,R*Q*R');
+    Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium);
     Pinf	= [];
   elseif options_.lik_init == 2	% Old Diffuse Kalman filter
     Pstar = 10*eye(np);
@@ -183,7 +183,7 @@ function [fval,llik,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend
         R1(j,:) = R1(j,:)-RR(j,i).*R(ivd(i),:);
       end
     end
-    Pstar = lyapunov_symm(T0,R1*Q*R1');
+    Pstar = lyapunov_symm(T0,R1*Q*R1',options_.qz_criterium);
   end
   %------------------------------------------------------------------------------
   % 4. Likelihood evaluation
diff --git a/matlab/DsgeSmoother.m b/matlab/DsgeSmoother.m
index 37fc4245b..d5d5b59a7 100644
--- a/matlab/DsgeSmoother.m
+++ b/matlab/DsgeSmoother.m
@@ -1,6 +1,4 @@
 function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,Y)
-
-% function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,decomp] = DsgeSmoother(xparam1,gend,Y)
 % Estimation of the smoothed variables and innovations. 
 % 
 % INPUTS 
@@ -95,7 +93,7 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,
   H = M_.H;
   
   if options_.lik_init == 1		% Kalman filter
-    Pstar = lyapunov_symm(T,R*Q*transpose(R));
+    Pstar = lyapunov_symm(T,R*Q*transpose(R),options_.qz_criterium);
     Pinf	= [];
   elseif options_.lik_init == 2 % Old Diffuse Kalman filter
     Pstar = 10*eye(np);
@@ -105,7 +103,7 @@ function [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK,T,R,P,PK,d,
       Pstar = zeros(np,np);
       ivs = bayestopt_.restrict_var_list_stationary;
       R1 = R(ivs,:);
-      Pstar(ivs,ivs) = lyapunov_symm(T(ivs,ivs),R1*Q*R1');
+      Pstar(ivs,ivs) = lyapunov_symm(T(ivs,ivs),R1*Q*R1',options_.qz_criterium);
       %    Pinf  = bayestopt_.Pinf;
       % by M. Ratto
       RR=T(:,bayestopt_.restrict_var_list_nonstationary);
diff --git a/matlab/DsgeVarLikelihood.m b/matlab/DsgeVarLikelihood.m
index d0e5b1c41..757d2830b 100644
--- a/matlab/DsgeVarLikelihood.m
+++ b/matlab/DsgeVarLikelihood.m
@@ -122,7 +122,7 @@ end
 %------------------------------------------------------------------------------
 % 3. theoretical moments (second order)
 %------------------------------------------------------------------------------
-tmp0 = lyapunov_symm(T,R*Q*R');% I compute the variance-covariance matrix
+tmp0 = lyapunov_symm(T,R*Q*R',options_.qz_criterium);% I compute the variance-covariance matrix
 mf  = bayestopt_.mf1;          % of the restricted state vector.
 
 % Get the non centered second order moments
diff --git a/matlab/UnivariateSpectralDensity.m b/matlab/UnivariateSpectralDensity.m
index 6107b0cb9..2bfbb5e1c 100644
--- a/matlab/UnivariateSpectralDensity.m
+++ b/matlab/UnivariateSpectralDensity.m
@@ -7,6 +7,7 @@ function [omega,f] = UnivariateSpectralDensity(dr,var_list)
 % Adapted from th_autocovariances.m.  
 global options_ oo_ M_
 
+
 omega = []; f = [];
 
 if options_.order > 1
@@ -17,8 +18,8 @@ if options_.order > 1
 end
 
 pltinfo  = 1;%options_.SpectralDensity.Th.plot;
-cutoff   = 100;%options_.SpectralDensity.Th.cutoff;
-sdl      = 0.1;%options_.SepctralDensity.Th.sdl;
+cutoff   = 150;%options_.SpectralDensity.Th.cutoff;
+sdl      = 0.01;%options_.SepctralDensity.Th.sdl;
 omega    = (0:sdl:pi)';
 GridSize = length(omega);
 exo_names_orig_ord  = M_.exo_names_orig_ord;
@@ -76,34 +77,70 @@ for i=M_.maximum_lag:-1:2
   AS(:,j1) = AS(:,j1)+ghx(:,i1);
   i0 = i1;
 end
-b = ghu1*M_.Sigma_e*ghu1';
-[A,B] = kalman_transition_matrix(dr);
-% index of predetermined variables in A
-i_pred = [nstatic+(1:npred) M_.endo_nbr+1:length(A)];
-A = A(i_pred,i_pred);
-[vx, ns_var] =  lyapunov_symm(A,b);
-i_ivar = find(~ismember(ivar,dr.order_var(ns_var+nstatic)));
-ivar = ivar(i_ivar);
-iky = iv(ivar);  
+Gamma = zeros(nvar,cutoff+1);
+[A,B] = kalman_transition_matrix(dr,ikx',1:nx,dr.transition_auxiliary_variables,M_.exo_nbr);
+[vx, u] =  lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium);
+iky = iv(ivar);
+if ~isempty(u)
+    iky = iky(find(any(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));
+    ivar = dr.order_var(iky);
+end
 aa = ghx(iky,:);
 bb = ghu(iky,:);
-Gamma = zeros(nvar,cutoff+1);
-tmp = aa*vx*aa'+ bb*M_.Sigma_e*bb';
-k = find(abs(tmp) < 1e-12);
-tmp(k) = 0;
-Gamma(:,1) = diag(tmp);
-vxy = (A*vx*aa'+ghu1*M_.Sigma_e*bb');
-tmp = aa*vxy;
-k = find(abs(tmp) < 1e-12);
-tmp(k) = 0;
-Gamma(:,2) = diag(tmp);
-for i=2:cutoff
-  vxy = A*vxy;
-  tmp = aa*vxy;
-  k = find(abs(tmp) < 1e-12);
-  tmp(k) = 0;
-  Gamma(:,i+1) = diag(tmp);
+
+if options_.hp_filter == 0
+    tmp = aa*vx*aa'+ bb*M_.Sigma_e*bb';
+    k = find(abs(tmp) < 1e-12);
+    tmp(k) = 0;
+    Gamma(:,1) = diag(tmp);
+    vxy = (A*vx*aa'+ghu1*M_.Sigma_e*bb');
+    tmp = aa*vxy;
+    k = find(abs(tmp) < 1e-12);
+    tmp(k) = 0;
+    Gamma(:,2) = diag(tmp);
+    for i=2:cutoff
+        vxy = A*vxy;
+        tmp = aa*vxy;
+        k = find(abs(tmp) < 1e-12);
+        tmp(k) = 0;
+        Gamma(:,i+1) = diag(tmp);
+    end
+else
+    iky = iv(ivar);  
+    aa = ghx(iky,:);
+    bb = ghu(iky,:);
+    lambda = options_.hp_filter;
+    ngrid = options_.hp_ngrid;
+    freqs = 0 : ((2*pi)/ngrid) : (2*pi*(1 - .5/ngrid)); 
+    tpos  = exp( sqrt(-1)*freqs);
+    tneg  =  exp(-sqrt(-1)*freqs);
+    hp1 = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2);
+    mathp_col = [];
+    IA = eye(size(A,1));
+    IE = eye(M_.exo_nbr);
+    for ig = 1:ngrid
+        f_omega  =(1/(2*pi))*( [inv(IA-A*tneg(ig))*ghu1;IE]...
+                               *M_.Sigma_e*[ghu1'*inv(IA-A'*tpos(ig)) IE]); % state variables
+        g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % selected variables
+        f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
+        mathp_col = [mathp_col ; (f_hp(:))'];    % store as matrix row
+                                                 % for ifft
+    end;
+    imathp_col = real(ifft(mathp_col))*(2*pi);
+    tmp = reshape(imathp_col(1,:),nvar,nvar);
+    k = find(abs(tmp)<1e-12);
+    tmp(k) = 0;
+    Gamma(:,1) = diag(tmp);
+    sy = sqrt(Gamma(:,1));
+    sy = sy *sy';
+    for i=1:cutoff-1
+        tmp = reshape(imathp_col(i+1,:),nvar,nvar)./sy;
+        k = find(abs(tmp) < 1e-12);
+        tmp(k) = 0;
+        Gamma(:,i+1) = diag(tmp);
+    end
 end
+
 H = 1:cutoff;
 for i=1:nvar
   f(i,:) = Gamma(i,1)/(2*pi) + Gamma(i,H+1)*cos(H'*omega')/pi;
diff --git a/matlab/calib_obj.m b/matlab/calib_obj.m
index 018eb4726..8a5820121 100644
--- a/matlab/calib_obj.m
+++ b/matlab/calib_obj.m
@@ -1,7 +1,7 @@
 % targets and iy order: 1) variances 2) correlations 
 % 3) constraints on M_.Sigma_e itself 4) autocorrelations
 function f=calib_obj(M_.Sigma_e,A,ghu1,ghx,ghu,targets,var_weights,iy,nar)
-  global vx fold
+  global vx fold options_
   
   oo_.gamma_y = cell(nar+1,1);
 %  M_.Sigma_e = M_.Sigma_e'*M_.Sigma_e;
@@ -10,11 +10,11 @@ function f=calib_obj(M_.Sigma_e,A,ghu1,ghx,ghu,targets,var_weights,iy,nar)
   b=ghu1*M_.Sigma_e*ghu1';
   vx = [];
   if isempty(vx)
-    vx = lyapunov_symm(A,b);
+    vx = lyapunov_symm(A,b,options_.qz_criterium);
   else
     [vx,status] = bicgstab_(@f_var,b(:),vx(:),1e-8,50,A,nx);
     if status
-      vx = lyapunov_symm(A,b);
+      vx = lyapunov_symm(A,b,options_.qz_criterium);
     else
       vx=reshape(vx,nx,nx);
     end
diff --git a/matlab/calib_obj2.m b/matlab/calib_obj2.m
index 4b5993b72..4fbac55c2 100644
--- a/matlab/calib_obj2.m
+++ b/matlab/calib_obj2.m
@@ -1,14 +1,14 @@
 % targets and iy order: 1) variances 2) correlations 
 % 3) constraints on M_.Sigma_e itself 4) autocorrelations
 function objective=calib_obj2(M_.Sigma_e,A,ghu1,ghx,ghu,targets,var_weights,iy,nar)
-  global vx fold
+  global vx fold options_
   
   objective = cell (nar+3);
   oo_.gamma_y = cell(nar+1,1);
   M_.Sigma_e=diag(M_.Sigma_e);
   nx = size(ghx,2);
   b=ghu1*M_.Sigma_e*ghu1';
-  vx = lyapunov_symm(A,b);
+  vx = lyapunov_symm(A,b,options_.qz_criterium);
   oo_.gamma_y{1} = ghx*vx*ghx'+ ghu*M_.Sigma_e*ghu';
   if ~isempty(targets{1})
     objective{1} = sqrt(oo_.gamma_y{1}(iy{1}));
diff --git a/matlab/check_posterior_analysis_data.m b/matlab/check_posterior_analysis_data.m
new file mode 100644
index 000000000..fca13bdda
--- /dev/null
+++ b/matlab/check_posterior_analysis_data.m
@@ -0,0 +1,80 @@
+function [info,description] = check_posterior_analysis_data(type,M_)
+% part of DYNARE, copyright Dynare Team (2008)
+% Gnu Public License.
+    
+    info = 0;
+    if nargout>1
+        description = '';
+    end
+    
+    %% Get informations about mcmc files. 
+    if ~exist([ M_.dname '/metropolis'],'dir')
+        disp('check_posterior_analysis_data:: Can''t find any mcmc file!')
+        return
+    end
+    mhname = get_name_of_the_last_mh_file(M_);
+    mhdate = get_date_of_a_file(mhname);
+    
+    %% Get informations about _posterior_draws files.
+    if ~exist([ M_.dname '/metropolis/' M_.fname '_posterior_draws.mat'])
+        info = 1; % select_posterior_draws has to be called first.
+        if nargout>1
+            description = 'select_posterior_draws has to be called.';
+        end
+        return
+    else
+        pddate = get_date_of_a_file([ M_.dname '/metropolis/' M_.fname '_posterior_draws.mat']);
+        if pddate<mhdate
+            info = 2; % _posterior_draws files have to be updated.
+            if nargout>1
+                description = 'posterior draws files have to be updated.';
+            end
+            return
+        else
+            info = 3; % Ok!
+            if nargout>1
+                description = 'posterior draws files are up to date.';
+            end
+        end
+    end
+    
+    %% Get informations about posterior data files.
+    switch type
+      case 'variance'
+        generic_post_data_file_name = 'Posterior2ndOrderMoments';
+      case 'decomposition'
+        generic_post_data_file_name = 'PosteriorVarianceDecomposition';
+      case 'dynamic_decomposition'
+        generic_post_data_file_name = 'PosteriorDynamicVarianceDecomposition';
+      otherwise
+        disp('This feature is not yest implemented!')
+    end
+    pdfinfo = dir([ M_.dname '/metropolis/' M_.fname '_' generic_post_data_file_name '*'])
+    if isempty(pdfinfo)
+        info = 4; % posterior draws have to be processed.
+        if nargout>1
+            description = 'posterior draws files have to be processed.';
+        end
+        return
+    else
+        number_of_the_last_post_data_file = length(pdfinfo);
+        name_of_the_last_post_data_file = ...
+            [ M_.dname ...
+              '/metropolis/' ...
+              M_.dname ...
+              generic_post_data_file_name ...
+              int2str(number_of_the_last_post_data_file) ...
+              '.mat' ];
+        pdfdate = get_date_of_a_file(name_of_the_last_post_data_file);
+        if pdfdate<pddate
+            info = 5; % posterior data files have to be updated.
+            if nargout>1
+                description = 'posterior data files have to be updated.';
+            end            
+        else
+            info = 6; % Ok (nothing to do ;-)
+            if nargout>1
+                description = 'There is nothing to do';
+            end        
+        end
+    end
\ No newline at end of file
diff --git a/matlab/compute_model_moments.m b/matlab/compute_model_moments.m
index 57399adab..c8b551bd9 100644
--- a/matlab/compute_model_moments.m
+++ b/matlab/compute_model_moments.m
@@ -1,11 +1,9 @@
-function moments=compute_model_moments(dr,options)
-
-% function compute_model_moments(options)
-% Computes posterior filter smoother and forecasts
+function moments=compute_model_moments(dr,M_,options,_)
 %
 % INPUTS
 %    dr:        structure describing model solution
-%    options:   structure of Dynare options
+%    M_:   structure of Dynare options
+%     options_
 %    
 % OUTPUTS
 %    moments: a cell array containing requested moments
@@ -16,7 +14,4 @@ function moments=compute_model_moments(dr,options)
 % part of DYNARE, copyright Dynare Team (2008)
 % Gnu Public License.
 
-% subset of variables
-    varlist = options.varlist;
-
-    moments = th_autocovariances(dr,varlist);
+    moments = th_autocovariances(dr,options_.varlist,M_,options_);
diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m
index 240fc0ca8..56683339c 100644
--- a/matlab/disp_th_moments.m
+++ b/matlab/disp_th_moments.m
@@ -19,7 +19,7 @@ function disp_th_moments(dr,var_list)
     end
   end
   
-  [oo_.gamma_y,ivar] = th_autocovariances(dr,ivar);
+  [oo_.gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_);
   m = dr.ys(ivar);
 
   
diff --git a/matlab/dsge_posterior_theoretical_covariance.m b/matlab/dsge_posterior_theoretical_covariance.m
index dc45a3a41..feaebe018 100644
--- a/matlab/dsge_posterior_theoretical_covariance.m
+++ b/matlab/dsge_posterior_theoretical_covariance.m
@@ -1,10 +1,12 @@
-function dsge_posterior_theoretical_covariance()
+function [nvar,vartan,CovarFileNumber] = dsge_posterior_theoretical_covariance(SampleSize,M_,options_,oo_)
 % This function estimates the posterior density of the endogenous
 % variables second order moments. 
 % 
 % INPUTS 
-%   None.
-%  
+%   SampleSize   [integer]
+%   
+%
+%
 % OUTPUTS 
 %   None.
 %
@@ -21,52 +23,30 @@ function dsge_posterior_theoretical_covariance()
 % part of DYNARE, copyright Dynare Team (2007-2008)
 % Gnu Public License.
 
-global M_ options_ oo_
-
-type = 'posterior';% To be defined as a input argument later...
-NumberOfSimulations = 1000;% To be defined in a global structure...
-
+type = 'posterior';
+    
 % Set varlist (vartan)
 [ivar,vartan] = set_stationary_variables_list;
-
-% Set various parameters & Check or create files and directories
-if strcmpi(type,'posterior')
-    MhDirectoryName = CheckPath('metropolis');
-else
-    MhDirectoryName = CheckPath('prior');
-end
-fname = [ MhDirectoryName '/' M_.fname];
-%save([fname '_Posterior2ndOrder'],'varlist');
-DrawsFiles = dir([fname '_' type '_draws*' ]);
-if ~rows(DrawsFiles)
-    if strcmpi(type,'posterior')
-        drsize = size_of_the_reduced_form_model(oo_.dr);
-        if drsize*NumberOfSimulations>101 %Too big model!
-            drsize=0;
-        end
-        SampleAddress = selec_posterior_draws(NumberOfSimulations,drsize);
-    else
-        % (samples from the prior) To be done later...
-    end
-    DrawsFiles = dir([fname '_' type '_draws*']);
-end
-
-% Get the number of stationary endogenous variables.
 nvar = length(ivar);
 
-nar = options_.ar;% Saves size of the auto-correlation function.
-options_.ar = 0;% Set the size of the auto-correlation function to zero.
+% Set the size of the auto-correlation function to zero.
+nar = options_.ar;
+options_.ar = 0;    
 
-NumberOfDrawsFiles = rows(DrawsFiles);
+% Get informations about the _posterior_draws files.
+DrawsFiles = dir([M_.dname '/metropolis/' M_.fname '_' type '_draws*' ]);
+NumberOfDrawsFiles = length(DrawsFiles);
+
+% Number of lines in posterior data files.
 MaXNumberOfCovarLines = ceil(options_.MaxNumberOfBytes/(nvar*(nvar+1)/2)/8);
 
-if NumberOfSimulations<=MaXNumberOfCovarLines
+if SampleSize<=MaXNumberOfCovarLines
     Covariance_matrix = zeros(NumberOfSimulations,nvar*(nvar+1)/2);
     NumberOfCovarFiles = 1;
 else
     Covariance_matrix = zeros(MaXNumberOfCovarLines,nvar*(nvar+1)/2);
-    NumberOfLinesInTheLastCovarFile = mod(NumberOfSimulations,MaXNumberOfCovarLines);
-    NumberOfCovarFiles = ceil(NumberOfSimulations/MaXNumberOfCovarLines);
+    NumberOfLinesInTheLastCovarFile = mod(SampleSize,MaXNumberOfCovarLines);
+    NumberOfCovarFiles = ceil(SampleSize/MaXNumberOfCovarLines);
 end
 
 NumberOfCovarLines = rows(Covariance_matrix);
@@ -86,10 +66,10 @@ for file = 1:NumberOfDrawsFiles
             set_parameters(pdraws{linee,1});
             [dr,info] = resol(oo_.steady_state,0);
         end
-        tmp = th_autocovariances(dr,ivar);
+        tmp = th_autocovariances(dr,ivar,M_,options_);
         for i=1:nvar
             for j=i:nvar
-                Covariance_matrix(linea,idx(i,j,nvar)) = tmp{1}(i,j);
+                Covariance_matrix(linea,symmetric_matrix_index(i,j,nvar)) = tmp{1}(i,j);
             end
         end
         if linea == NumberOfCovarLines
@@ -109,45 +89,12 @@ for file = 1:NumberOfDrawsFiles
         end
     end
 end
-options_.ar = nar; clear('pdraws','tmp');
 
-% Compute statistics and save in oo_
-for i=1:nvar
-    for j=i:nvar
-        i1 = 1;
-        tmp = zeros(NumberOfSimulations,1);
-        for file = 1:CovarFileNumber
-            load([fname '_Posterior2ndOrderMoments' int2str(file)]);
-            i2 = i1 + rows(Covariance_matrix) - 1;
-            tmp(i1:i2) = Covariance_matrix(:,idx(i,j,nvar));
-            i1 = i2+1;
-        end
-        [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
-            posterior_moments(tmp,1,options_.mh_conf_sig);
-        name = fieldname(i,j,vartan);
-        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.mean.' name ' = post_mean;']);
-        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.median.' name ' = post_median;']);
-        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.variance.' name ' = post_var;']);
-        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.hpdinf.' name ' = hpd_interval(1);']);
-        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.hpdsup.' name ' = hpd_interval(2);']);
-        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.deciles.' name ' = post_deciles;']);
-        eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.density.' name ' = density;']);
-    end 
-end
+options_.ar = nar;
 
 
-
-    
-function k = idx(i,j,n)
-    k = (i-1)*n+j-i*(i-1)/2;
-
 function r = rows(M)
     r = size(M,1);
 
 function c = cols(M)
-    c = size(M,2);
-    
-function name = fieldname(i,j,vlist)
-    n1 = deblank(vlist(i,:));
-    n2 = deblank(vlist(j,:));
-    name = [n1 '.' n2];
\ No newline at end of file
+    c = size(M,2);
\ No newline at end of file
diff --git a/matlab/dsge_posterior_theoretical_variance_decomposition.m b/matlab/dsge_posterior_theoretical_variance_decomposition.m
index e8914c061..2749d0bc5 100644
--- a/matlab/dsge_posterior_theoretical_variance_decomposition.m
+++ b/matlab/dsge_posterior_theoretical_variance_decomposition.m
@@ -1,4 +1,5 @@
-function dsge_posterior_theoretical_variance_decomposition()
+function [nvar,vartan,NumberOfDecompFiles] = ...
+        dsge_posterior_theoretical_variance_decomposition(SampleSize,M_,options_,oo_)
 % This function estimates the posterior distribution of the variance
 % decomposition of the observed endogenous variables.
 % 
@@ -20,42 +21,22 @@ function dsge_posterior_theoretical_variance_decomposition()
 % part of DYNARE, copyright Dynare Team (2007-2008).
 % Gnu Public License.
 
-global M_ options_ oo_
-
 type = 'posterior';% To be defined as a input argument later...
-NumberOfSimulations = 800;% To be defined in a global structure...
 
 % Set varlist (vartan)
 [ivar,vartan] = set_stationary_variables_list;
-ivar
-vartan
+nvar = length(ivar);
 
-% Set various parameters, Check or create files and directories &
-% initialize arrays.
-if strcmpi(type,'posterior')
-    MhDirectoryName = CheckPath('metropolis');
-else
-    MhDirectoryName = CheckPath('prior');
-end
-fname = [ MhDirectoryName '/' M_.fname];
-DrawsFiles = dir([fname '_' type '_draws*' ]);
-if ~rows(DrawsFiles)
-    if strcmpi(type,'posterior')
-        drsize = size_of_the_reduced_form_model(oo_.dr);
-        if drsize*NumberOfSimulations>101%Big model!
-            drsize=0;
-        end
-        SampleAddress = selec_posterior_draws(NumberOfSimulations,drsize);
-    else% (samples from the prior) To be done later...
-    end
-    DrawsFiles = dir([fname '_' type '_draws*']);
-end
+% Set the size of the auto-correlation function to zero.
+nar = options_.ar;
+options_.ar = 0;    
 
-nar = options_.ar;% Saves size of the auto-correlation function.
-options_.ar = 0;% Set the size of the auto-correlation function.
+% Get informations about the _posterior_draws files.
+DrawsFiles = dir([M_.dname '/metropolis/' M_.fname '_' type '_draws*' ]);
+NumberOfDrawsFiles = length(DrawsFiles);
 
 nexo = M_.exo_nbr;
-nvar = length(ivar);
+
 
 NumberOfDrawsFiles = rows(DrawsFiles);
 NumberOfSavedElementsPerSimulation = nvar*(nexo+1);
@@ -89,10 +70,11 @@ for file = 1:NumberOfDrawsFiles
             set_parameters(pdraws{linee,1});
             [dr,info] = resol(oo_.steady_state,0);
         end
-        tmp = th_autocovariances(dr,ivar);
-        for i=1:nvar
-            Decomposition_array(linea,i) = tmp{1}(i,i);
-        end
+        tmp = th_autocovariances(dr,ivar,M_,options_);
+        %for i=1:nvar
+        %    Decomposition_array(linea,i) = tmp{1}(i,i);
+        %end
+        Decomposition_array(linea,:) = transpose(tmp{1});
         for i=1:nvar
             for j=1:nexo
                 Decomposition_array(linea,nvar+(i-1)*nexo+j) = tmp{2}(i,j);
@@ -115,62 +97,11 @@ for file = 1:NumberOfDrawsFiles
         end
     end
 end
-options_.ar = nar; clear('pdraws','tmp');
 
-% Compute statistics and save in oo_
-
-for i=1:nvar
-    for j=1:nexo
-        i1 = 1;
-        tmp = zeros(NumberOfSimulations,1);
-        for file = 1:DecompFileNumber
-            load([fname '_PosteriorVarianceDecomposition' int2str(file)]);
-            i2 = i1 + rows(Decomposition_array) - 1;
-            tmp(i1:i2) = Decomposition_array(:,nvar+(i-1)*nexo+j);
-            i1 = i2+1;
-        end
-        name = [ deblank(vartan(i,:)) '.' deblank(M_.exo_names(j,:)) ];
-        t1 = min(tmp); t2 = max(tmp);
-        t3 = t2-t1;% How to normalize ? t1 and t2 may be zero...
-        if t3<1.0e-12
-            if t1<1.0e-12
-                t1 = 0;
-            end
-            if abs(t1-1)<1.0e-12
-                t1 = 1;
-            end 
-            post_mean = t1;
-            post_median = t1;
-            post_var = 0;
-            hpd_interval = NaN(2,1);
-            post_deciles = NaN(9,1);
-            density = NaN;
-        else
-            [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
-                posterior_moments(tmp,1,options_.mh_conf_sig);
-        end
-        
-        eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.mean.' name ' = post_mean;']);
-        eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.median.' name ' = post_median;']);
-        eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.variance.' name ' = post_var;']);
-        eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.hpdinf.' name ' = hpd_interval(1);']);
-        eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.hpdsup.' name ' = hpd_interval(2);']);
-        eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.deciles.' name ' = post_deciles;']);
-        eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.density.' name ' = density;']);
-    end
-end
-
-
-function k = idx(i,j,n)
-    k = (i-1)*n+j-i*(i-1)/2;
+options_.ar = nar;% Useless because options_ is not a global anymore...
 
 function r = rows(M)
     r = size(M,1);
 
 function c = cols(M)
-    c = size(M,2);
-    
-function name = fieldname(i,j,vlist)
-    n1 = deblank(vlist(i,:));
-    n2 = deblank(vlist(j,:));
-    name = [n1 '.' n2];
\ No newline at end of file
+    c = size(M,2);
\ No newline at end of file
diff --git a/matlab/dynare_resolve.m b/matlab/dynare_resolve.m
index 1c6c9fecd..abd9586d1 100644
--- a/matlab/dynare_resolve.m
+++ b/matlab/dynare_resolve.m
@@ -32,7 +32,7 @@ function [A,B,ys,info] = dynare_resolve(iv,ic,aux)
   global oo_ M_
   
   [oo_.dr,info] = resol(oo_.steady_state,0);
-  
+
   if info(1) > 0
     A = [];
     B = [];
@@ -52,5 +52,5 @@ function [A,B,ys,info] = dynare_resolve(iv,ic,aux)
     aux(k,2) = aux(k,2) + oo_.dr.nfwrd;
   end
   
-  [A,B] = kalman_transition_matrix(oo_.dr,iv,ic,aux);
+  [A,B] = kalman_transition_matrix(oo_.dr,iv,ic,aux,M_.exo_nbr);
   ys = oo_.dr.ys;
\ No newline at end of file
diff --git a/matlab/forcst.m b/matlab/forcst.m
index d299da5ed..986f60989 100644
--- a/matlab/forcst.m
+++ b/matlab/forcst.m
@@ -30,7 +30,7 @@ function [yf,int_width]=forcst(dr,y0,horizon,var_list)
     nc = size(dr.ghx,2);
     endo_nbr = M_.endo_nbr;
     inv_order_var = dr.inv_order_var;
-    [A,B] = kalman_transition_matrix(dr,nstatic+(1:npred),1:nc,dr.transition_auxiliary_variables);
+    [A,B] = kalman_transition_matrix(dr,nstatic+(1:npred),1:nc,dr.transition_auxiliary_variables,M_.exo_nbr);
     
     nvar = size(var_list,1);
     if nvar == 0
diff --git a/matlab/get_date_of_a_file.m b/matlab/get_date_of_a_file.m
new file mode 100644
index 000000000..30cf5de33
--- /dev/null
+++ b/matlab/get_date_of_a_file.m
@@ -0,0 +1,8 @@
+function [d1,d2] = get_date_of_a_file(filename)
+% part of DYNARE, copyright Dynare Team (2008)
+% Gnu Public License.
+    info = dir(filename);
+    d1 = info.datenum;
+    if nargout>1
+        d2 = info.date;
+    end
\ No newline at end of file
diff --git a/matlab/get_innovation_contemporaneous_impact.m b/matlab/get_innovation_contemporaneous_impact.m
index dfb822892..dd1d5ecc8 100644
--- a/matlab/get_innovation_contemporaneous_impact.m
+++ b/matlab/get_innovation_contemporaneous_impact.m
@@ -1,4 +1,4 @@
-function F = get_innovation_contemporaneous_impact('type')
+function B = get_innovation_contemporaneous_impact(type,info)
 
 % function F = get_innovation_contemporaneous_impact('type')
 % The approximated reduced form model is 
@@ -15,7 +15,8 @@ function F = get_innovation_contemporaneous_impact('type')
 % given by F = Z*B. Matrix F is returned by this function.      
 %
 % INPUTS 
-%   o type = "mode","mean"
+%   o type   [string]      "mode","mean"
+%   o info   [integer]     if equal to 1, matrix B is saved in a mat file.
 %  
 % OUTPUTS 
 %   o F (F is also saved in a file)
@@ -26,17 +27,21 @@ function F = get_innovation_contemporaneous_impact('type')
 % part of DYNARE, copyright Dynare Team (2006-2008)
 % Gnu Public License.
 
-global oo_ M_ bayestopt_
+global oo_ M_ bayestopt_ options_
 
 if nargin == 0
     type = 'mode';
 end
 
+if nargin == 1
+    info = 0;
+end
+
 get_posterior_parameters(type);
 
-[dr,info]=dr1(oo_.dr,0);
+[dr,info,M_,options_,oo_]=dr1(oo_.dr,0,M_,options_,oo_);
 
 B(dr.order_var,M_.exo_names_orig_ord) = dr.ghu*sqrt(M_.Sigma_e);
-F = B(bayestopt_.mfys,:);
+B = B(bayestopt_.mfys,:);
 
-save([M_.fname '_InnovImpact',F]);
\ No newline at end of file
+save([M_.fname '_InnovImpact'],'B');
\ No newline at end of file
diff --git a/matlab/get_moments_size.m b/matlab/get_moments_size.m
index fd07ef3ca..eae42c5ca 100644
--- a/matlab/get_moments_size.m
+++ b/matlab/get_moments_size.m
@@ -1,5 +1,4 @@
 function s=get_moments_size(options)
-
 % function PosteriorFilterSmootherAndForecast(Y,gend, type)
 % Computes posterior filter smoother and forecasts
 %
diff --git a/matlab/get_name_of_the_last_mh_file.m b/matlab/get_name_of_the_last_mh_file.m
new file mode 100644
index 000000000..7079a8c10
--- /dev/null
+++ b/matlab/get_name_of_the_last_mh_file.m
@@ -0,0 +1,17 @@
+function mhname = get_name_of_the_last_mh_file(M_)
+% part of DYNARE, copyright Dynare Team (2008)
+% Gnu Public License.
+    model_name = M_.fname ;
+    mcmc_directory = M_.dname ;
+    load([ mcmc_directory '/metropolis/' model_name '_mh_history']) ;
+    mh_number = record.LastFileNumber ;
+    bk_number = record.Nblck ;
+    clear('record') ;
+    mhname = [ mcmc_directory ...
+               '/metropolis/' ...
+               model_name ...
+               '_mh' ...
+               int2str(mh_number) ...
+               '_blck' ...
+               int2str(bk_number) ...
+               '.mat'] ;
\ No newline at end of file
diff --git a/matlab/get_variance_of_endogenous_variables.m b/matlab/get_variance_of_endogenous_variables.m
index 2585bd264..cafc99ee2 100644
--- a/matlab/get_variance_of_endogenous_variables.m
+++ b/matlab/get_variance_of_endogenous_variables.m
@@ -29,9 +29,9 @@ function [vx1,i_ns] = get_variance_of_endogenous_variables(dr,i_var)
   nc = size(ghx,2);
   n = length(i_var);
   
-  [A,B] = kalman_transition_matrix(dr,nstatic+(1:npred),1:nc,dr.transition_auxiliary_variables);
+  [A,B] = kalman_transition_matrix(dr,nstatic+(1:npred),1:nc,dr.transition_auxiliary_variables,M_.exo_nbr);
   
-  [vx,u] = lyapunov_symm(A,B*Sigma_e*B');
+  [vx,u] = lyapunov_symm(A,B*Sigma_e*B',options_.qz_criterium);
   
   if size(u,2) > 0
     i_stat = find(any(abs(ghx*u) < options_.Schur_vec_tol,2));
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 9d8acb7f0..76708d449 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -145,6 +145,8 @@ function global_initialization()
   options_.cutoff = 1e-12;
   options_.student_degrees_of_freedom = 3;
   options_.subdraws = [];
+  options_.PosteriorSampleSize = 1000;
+  options_.MaximumNumberOfMegaBytes = 111;
   
   % Misc
   options_.conf_sig = 0.6;
diff --git a/matlab/kalman_transition_matrix.m b/matlab/kalman_transition_matrix.m
index 2a26a1671..b43512f80 100644
--- a/matlab/kalman_transition_matrix.m
+++ b/matlab/kalman_transition_matrix.m
@@ -1,14 +1,13 @@
-function [A,B] = kalman_transition_matrix(dr,iv,ic,aux)
-
-% function [A,B] = kalman_transition_matrix(dr,iv,ic,aux)
-% makes transition matrices out of ghx and ghu for Kalman filter
+function [A,B] = kalman_transition_matrix(dr,iv,ic,aux,exo_nbr)
+% Builds the transition equation of the state space representation out of ghx and ghu for Kalman filter
 % 
 % INPUTS
 %    dr:      structure of decisions rules for stochastic simulations
 %    iv:      selected variables
 %    ic:      state variables position in the transition matrix columns
 %    aux:     indices for auxiliary equations
-%
+%    exo_nbr: number of exogenous variables
+%    
 % OUTPUTS
 %    A:       matrix of predetermined variables effects in linear solution (ghx)
 %    B:       matrix of shocks effects in linear solution (ghu)
@@ -18,15 +17,13 @@ function [A,B] = kalman_transition_matrix(dr,iv,ic,aux)
 %  
 % part of DYNARE, copyright Dynare Team (2003-2008)
 % Gnu Public License.
-
-  global M_
   
   n_iv = length(iv);
   n_ir1 = size(aux,1);
   nr = n_iv + n_ir1;
   
   A = zeros(nr,nr);
-  B = zeros(nr,M_.exo_nbr);
+  B = zeros(nr,exo_nbr);
   
   i_n_iv = 1:n_iv;
   A(i_n_iv,ic) = dr.ghx(iv,:);
diff --git a/matlab/lyapunov_symm.m b/matlab/lyapunov_symm.m
index 1b787c457..8227a22cd 100644
--- a/matlab/lyapunov_symm.m
+++ b/matlab/lyapunov_symm.m
@@ -1,7 +1,5 @@
-function [x,u]=lyapunov_symm(a,b)
-
-% function [x,u]=lyapunov_symm(a,b)
-% solves the Lyapunov equation x-a*x*a' = b, for b (and then x) symmetrical
+function [x,u]=lyapunov_symm(a,b,qz_criterium)
+% Solves the Lyapunov equation x-a*x*a' = b, for b (and then x) symmetrical
 % if a has some unit roots, the function computes only the solution of the stable subsystem
 %  
 % INPUTS:
@@ -21,9 +19,7 @@ function [x,u]=lyapunov_symm(a,b)
 %   needs Matlab version with ordeig function
 %
 % part of DYNARE, copyright Dynare Team (2006-2008)
-% Gnu Public License  
-
-  global options_ 
+% Gnu Public License   
 
   ns_var = [];
   u = [];
@@ -35,9 +31,9 @@ function [x,u]=lyapunov_symm(a,b)
   end
   [u,t] = schur(a);
   if exist('ordeig','builtin')
-    e1 = abs(ordeig(t)) > 2-options_.qz_criterium;
+    e1 = abs(ordeig(t)) > 2-qz_criterium;
   else
-    e1 = abs(my_ordeig(t)) > 2-options_.qz_criterium;
+    e1 = abs(my_ordeig(t)) > 2-qz_criterium;
   end
   k = sum(e1);
   if exist('ordschur','builtin')
diff --git a/matlab/posterior_analysis.m b/matlab/posterior_analysis.m
new file mode 100644
index 000000000..a3ff246c4
--- /dev/null
+++ b/matlab/posterior_analysis.m
@@ -0,0 +1,137 @@
+function posterior_analysis(type,arg1,arg2,options_,M_,oo_)  
+% part of DYNARE, copyright Dynare Team (2008)
+% Gnu Public License.
+    
+    info = check_posterior_analysis_data(type,M_);
+
+    switch info
+      case 0
+        disp('check_posterior_analysis_data:: Can''t find any mcmc file!')
+        error('Check the options of the estimation command...')
+      case {1,2}
+        SampleSize = options_.PosteriorSampleSize;
+        MaxMegaBytes = options_.MaximumNumberOfMegaBytes;
+        drsize = size_of_the_reduced_form_model(oo_.dr);
+        if drsize*SampleSize>MaxMegaBytes
+            drsize=0;
+        end
+        SampleAddress = selec_posterior_draws(SampleSize,drsize);
+      case {4,5}
+        switch type
+          case 'variance'
+            [nvar,vartan,NumberOfFiles] = ...
+                dsge_posterior_theoretical_covariance(SampleSize,M_,options_,oo_);
+          case 'decomposition'
+            [nvar,vartan,NumberOfFiles] = ...
+                dsge_posterior_theoretical_variance_decomposition(SampleSize,M_,options_,oo_);
+          otherwise
+            disp('Not yet implemented')
+        end
+      case 6
+        switch type
+          case 'variance'
+            covariance_posterior_analysis(NumberOfFiles,SampleSize,M_.dname,M_.fname,...
+                                          vartan,nvar,arg1,arg2);
+          case 'decomposition'
+            variance_decomposition_posterior_analysis(NumberOfFiles,SampleSize,M_.dname,M_.fname,...
+                                                      M_.exo_names,arg2,vartan,nvar,arg1);
+          otherwise
+            disp('Not yet implemented')
+        end
+      otherwise
+        error(['posterior_analysis:: Check_posterior_analysis_data gave a meaningless output!'])
+    end
+    
+        
+    
+    
+    
+function covariance_posterior_analysis(NumberOfFiles,NumberOfSimulations,dname,fname,vartan,nvar,var1,var2)
+    indx1 = check_name(vartan,var1)
+    if isempty(indx1)
+        disp(['posterior_analysis:: ' var1 ' is not a stationary endogenous variable!'])
+        return
+    end
+    if ~isempty(var2)
+        indx2 = check_name(vartan,var2)
+        if isempty(indx2)
+            disp(['posterior_analysis:: ' var2 ' is not a stationary endogenous variable!'])
+            return
+        end
+    else
+        indx2 = indx1;
+        var2 = var1;
+    end
+    i1 = 1; tmp = zeros(NumberOfSimulations,1);
+    for file = 1:CovarFileNumber
+        load([ dname '/metropolis/'  fname '_Posterior2ndOrderMoments' int2str(file)]);
+        i2 = i1 + rows(Covariance_matrix) - 1;
+        tmp(i1:i2) = Covariance_matrix(:,symmetric_matrix_index(indx1,indx2,nvar));
+        i1 = i2+1;
+    end
+    [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
+        posterior_moments(tmp,1,options_.mh_conf_sig);
+    if strcmpi(var1,var2)
+        name = var1;
+    else
+        name = [var1 '.' var2];
+    end
+    eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.mean.' name ' = post_mean;']);
+    eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.median.' name ' = post_median;']);
+    eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.variance.' name ' = post_var;']);
+    eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.hpdinf.' name ' = hpd_interval(1);']);
+    eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.hpdsup.' name ' = hpd_interval(2);']);
+    eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.deciles.' name ' = post_deciles;']);
+    eval(['oo_.PosteriorTheoreticalMoments.dsge.covariance.density.' name ' = density;']);
+    
+    
+function variance_decomposition_posterior_analysis(NumberOfFiles,NumberOfSimulations,dname,fname, ...
+                                          exonames,exo,vartan,nvar,var)
+    indx1 = check_name(vartan,var)
+    if isempty(indx)
+        disp(['posterior_analysis:: ' var ' is not a stationary endogenous variable!'])
+        return
+    end
+    jndx = check_name(exonames,exo);
+    if isempty(jndx)
+        disp(['posterior_analysis:: ' exo ' is not a declared exogenous variable!'])
+        return
+    end
+    i1 = 1; tmp = zeros(NumberOfSimulations,1);
+    for file = 1:NumberOfFiles
+        load([fname '_PosteriorVarianceDecomposition' int2str(file)]);
+        i2 = i1 + rows(Decomposition_array) - 1;
+        tmp(i1:i2) = Decomposition_array(:,nvar+(i-1)*nexo+j);
+        i1 = i2+1;
+    end
+    name = [ var '.' exo ];
+        t1 = min(tmp); t2 = max(tmp);
+        t3 = t2-t1;% How to normalize ? t1 and t2 may be zero...
+        if t3<1.0e-12
+            if t1<1.0e-12
+                t1 = 0;
+            end
+            if abs(t1-1)<1.0e-12
+                t1 = 1;
+            end 
+            post_mean = t1;
+            post_median = t1;
+            post_var = 0;
+            hpd_interval = NaN(2,1);
+            post_deciles = NaN(9,1);
+            density = NaN;
+        else
+            [post_mean, post_median, post_var, hpd_interval, post_deciles, density] = ...
+                posterior_moments(tmp,1,options_.mh_conf_sig);
+        end
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.mean.' name ' = post_mean;']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.median.' name ' = post_median;']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.variance.' name ' = post_var;']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.hpdinf.' name ' = hpd_interval(1);']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.hpdsup.' name ' = hpd_interval(2);']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.deciles.' name ' = post_deciles;']);
+        eval(['oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.density.' name ' = density;']);
+
+    
+function n = check_name(vartan,varname)
+    n = strmatch(varname,vartan,'exact')
\ No newline at end of file
diff --git a/matlab/prior_posterior_statistics.m b/matlab/prior_posterior_statistics.m
index 4ef5b7f1f..a48264dd0 100644
--- a/matlab/prior_posterior_statistics.m
+++ b/matlab/prior_posterior_statistics.m
@@ -127,12 +127,11 @@ for b=1:B
   set_all_parameters(deep);
   dr = resol(oo_.steady_state,0);
   if moments_varendo
-      stock_moments{irun(8)} = compute_model_moments(dr,options_);
+      stock_moments{irun(8)} = compute_model_moments(dr,M_,options_);
   end
   if run_smoother
       [alphahat,etahat,epsilonhat,ahat,SteadyState,trend_coeff,aK] = ...
           DsgeSmoother(deep,gend,Y);
-      
       if options_.loglinear
           stock_smooth(dr.order_var,:,irun(1)) = alphahat(1:endo_nbr,:)+ ...
               repmat(log(dr.ys(dr.order_var)),1,gend);
diff --git a/matlab/symmetric_matrix_index.m b/matlab/symmetric_matrix_index.m
new file mode 100644
index 000000000..b9da53d9a
--- /dev/null
+++ b/matlab/symmetric_matrix_index.m
@@ -0,0 +1,4 @@
+function k = symmetric_matrix_index(i,j,n)
+% part of DYNARE, copyright Dynare Team (2007-2008).
+% Gnu Public License.    
+    k = (i-1)*n+j-i*(i-1)/2;
\ No newline at end of file
diff --git a/matlab/th_autocovariances.m b/matlab/th_autocovariances.m
index 99cf42aeb..6fab02c85 100644
--- a/matlab/th_autocovariances.m
+++ b/matlab/th_autocovariances.m
@@ -1,14 +1,14 @@
-function [Gamma_y,ivar]=th_autocovariances(dr,ivar)
-
-% function [Gamma_y,ivar]=th_autocovariances(dr,ivar)
-% computes the theoretical auto-covariances, Gamma_y, for an AR(p) process 
+function [Gamma_y,ivar]=th_autocovariances(dr,ivar,M_,options_)
+% Computes the theoretical auto-covariances, Gamma_y, for an AR(p) process 
 % with coefficients dr.ghx and dr.ghu and shock variances Sigma_e_
 % for a subset of variables ivar (indices in lgy_)
 %  
 % INPUTS
 %   dr:         structure of decisions rules for stochastic simulations
 %   ivar:       subset of variables
-%  
+%   M_
+%   options_
+%    
 % OUTPUTS
 %   Gamma_y:    theoritical auto-covariances
 %   ivar:       subset of variables
@@ -19,9 +19,6 @@ function [Gamma_y,ivar]=th_autocovariances(dr,ivar)
 % part of DYNARE, copyright Dynare Team (2001-2008)
 % Gnu Public License.
 
-
-  global M_ options_
-
   exo_names_orig_ord  = M_.exo_names_orig_ord;
   if sscanf(version('-release'),'%d') < 13
     warning off
@@ -68,9 +65,9 @@ function [Gamma_y,ivar]=th_autocovariances(dr,ivar)
 
   ipred = nstatic+(1:npred)';
   % state space representation for state variables only
-  [A,B] = kalman_transition_matrix(dr,ipred,1:nx,dr.transition_auxiliary_variables);
+  [A,B] = kalman_transition_matrix(dr,ipred,1:nx,dr.transition_auxiliary_variables,M_.exo_nbr);
   if options_.order == 2 | options_.hp_filter == 0
-    [vx, u] =  lyapunov_symm(A,B*M_.Sigma_e*B');
+    [vx, u] =  lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium);
     iky = iv(ivar);
     if ~isempty(u)
       iky = iky(find(all(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));
@@ -113,10 +110,10 @@ function [Gamma_y,ivar]=th_autocovariances(dr,ivar)
       b1 = b1*cs;
       b2(:,exo_names_orig_ord) = ghu(iky,:);
       b2 = b2*cs;
-      vx  = lyapunov_symm(A,b1*b1');
+      vx  = lyapunov_symm(A,b1*b1',options_.qz_criterium);
       vv = diag(aa*vx*aa'+b2*b2');
       for i=1:M_.exo_nbr
-	vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)');
+	vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.qz_criterium);
 	Gamma_y{nar+2}(:,i) = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'))./vv;
       end
     end
@@ -161,9 +158,9 @@ function [Gamma_y,ivar]=th_autocovariances(dr,ivar)
     end
     
     %variance decomposition
-    if M_.exo_nbr > 1 
+    if M_.exo_nbr > 1
       Gamma_y{nar+2} = zeros(nvar,M_.exo_nbr);
-      SS(exo_names_orig_ord,exo_names_orig_ord)=M_.Sigma_e+1e-14*eye(M_.exo_nbr);
+      SS(exo_names_orig_ord,exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr);
       cs = chol(SS)';
       SS = cs*cs';
       b1(:,exo_names_orig_ord) = ghu1;