diff --git a/matlab/check_name.m b/matlab/check_name.m
new file mode 100644
index 0000000000000000000000000000000000000000..2b77b2938228974167322a65dcb9c06edee7b743
--- /dev/null
+++ b/matlab/check_name.m
@@ -0,0 +1,2 @@
+function n = check_name(vartan,varname)
+    n = strmatch(varname,vartan,'exact');
\ No newline at end of file
diff --git a/matlab/check_posterior_analysis_data.m b/matlab/check_posterior_analysis_data.m
index fca13bdda00681d5bc5495eca20843a491ddbf0b..0f7946e16e9a2dd7ad02335eeb9a6dd0e3af6da0 100644
--- a/matlab/check_posterior_analysis_data.m
+++ b/matlab/check_posterior_analysis_data.m
@@ -49,7 +49,7 @@ function [info,description] = check_posterior_analysis_data(type,M_)
       otherwise
         disp('This feature is not yest implemented!')
     end
-    pdfinfo = dir([ M_.dname '/metropolis/' M_.fname '_' generic_post_data_file_name '*'])
+    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
@@ -59,9 +59,9 @@ function [info,description] = check_posterior_analysis_data(type,M_)
     else
         number_of_the_last_post_data_file = length(pdfinfo);
         name_of_the_last_post_data_file = ...
-            [ M_.dname ...
+            [ './' M_.dname ...
               '/metropolis/' ...
-              M_.dname ...
+              M_.fname '_' ... 
               generic_post_data_file_name ...
               int2str(number_of_the_last_post_data_file) ...
               '.mat' ];
diff --git a/matlab/cols.m b/matlab/cols.m
new file mode 100644
index 0000000000000000000000000000000000000000..fa262da17da7ac1db49287f7ab75bcab4ba59eba
--- /dev/null
+++ b/matlab/cols.m
@@ -0,0 +1,2 @@
+function c = cols(M)
+    c = size(M,2);
\ No newline at end of file
diff --git a/matlab/covariance_posterior_analysis.m b/matlab/covariance_posterior_analysis.m
new file mode 100644
index 0000000000000000000000000000000000000000..aec229222e9b3697f490c5da0631945ca5801e4e
--- /dev/null
+++ b/matlab/covariance_posterior_analysis.m
@@ -0,0 +1,56 @@
+function oo_ = covariance_posterior_analysis(NumberOfSimulations,dname,fname,vartan,nvar,var1,var2,mh_conf_sig,oo_)
+    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
+    if isfield(oo_,'PosteriorTheoreticalMoments')
+        if isfield(oo_.PosteriorTheoreticalMoments,'dsge')
+            if isfield(oo_.PosteriorTheoreticalMoments.dsge,'covariance')
+                if isfield(oo_.PosteriorTheoreticalMoments.dsge.covariance.mean,var1)
+                    eval(['s1 = oo_.PosteriorTheoreticalMoments.dsge.covariance.mean' '.' var1 ';'])  
+                    if isfield(s1,var2)
+                        % Nothing to do.
+                        return
+                    end
+                else
+                    if isfield(oo_.PosteriorTheoreticalMoments.dsge.covariance.mean,var2)
+                        eval(['s2 = oo_.PosteriorTheoreticalMoments.dsge.covariance.mean' '.' var2 ';'])
+                        if isfield(s1,var1)
+                            % Nothing to do (the covariance matrix is symmetric!).
+                            return
+                        end
+                    end
+                end         
+            end
+        end
+    end
+    tmp = dir([ dname '/metropolis/'  fname '_Posterior2ndOrderMoments*.mat']);
+    NumberOfFiles= length(tmp);
+    i1 = 1; tmp = zeros(NumberOfSimulations,1);
+    for file = 1:NumberOfFiles
+        load([ dname '/metropolis/'  fname '_Posterior2ndOrderMoments' int2str(file) '.mat']);
+        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,mh_conf_sig);
+    name = [var1 '.' var2];
+    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;']);
\ No newline at end of file
diff --git a/matlab/dsge_posterior_theoretical_covariance.m b/matlab/dsge_posterior_theoretical_covariance.m
index feaebe0188e3c978e96243b04266c5649abadfa9..eae3552df54f69d87c27d6ec8f4b316aeee77c6a 100644
--- a/matlab/dsge_posterior_theoretical_covariance.m
+++ b/matlab/dsge_posterior_theoretical_covariance.m
@@ -41,7 +41,7 @@ NumberOfDrawsFiles = length(DrawsFiles);
 MaXNumberOfCovarLines = ceil(options_.MaxNumberOfBytes/(nvar*(nvar+1)/2)/8);
 
 if SampleSize<=MaXNumberOfCovarLines
-    Covariance_matrix = zeros(NumberOfSimulations,nvar*(nvar+1)/2);
+    Covariance_matrix = zeros(SampleSize,nvar*(nvar+1)/2);
     NumberOfCovarFiles = 1;
 else
     Covariance_matrix = zeros(MaXNumberOfCovarLines,nvar*(nvar+1)/2);
@@ -55,7 +55,7 @@ CovarFileNumber = 1;
 % Compute 2nd order moments and save them in *_Posterior2ndOrderMoments* files
 linea = 0;
 for file = 1:NumberOfDrawsFiles
-    load([MhDirectoryName '/' DrawsFiles(file).name]);
+    load([M_.dname '/metropolis/' DrawsFiles(file).name]);
     NumberOfDraws = rows(pdraws);
     isdrsaved = cols(pdraws)-1;
     for linee = 1:NumberOfDraws
@@ -73,7 +73,7 @@ for file = 1:NumberOfDrawsFiles
             end
         end
         if linea == NumberOfCovarLines
-            save([fname '_Posterior2ndOrderMoments' int2str(CovarFileNumber)],'Covariance_matrix');
+            save([ M_.dname '/metropolis/' M_.fname '_Posterior2ndOrderMoments' int2str(CovarFileNumber)],'Covariance_matrix');
             CovarFileNumber = CovarFileNumber + 1;
             linea = 0;
             test = CovarFileNumber-NumberOfCovarFiles;
@@ -90,11 +90,4 @@ for file = 1:NumberOfDrawsFiles
     end
 end
 
-options_.ar = nar;
-
-
-function r = rows(M)
-    r = size(M,1);
-
-function c = cols(M)
-    c = size(M,2);
\ No newline at end of file
+options_.ar = nar;
\ 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 2749d0bc56cf5aed61395f13418daad248218620..7bab2b570db6aff1d471c3ff475d702902379e6d 100644
--- a/matlab/dsge_posterior_theoretical_variance_decomposition.m
+++ b/matlab/dsge_posterior_theoretical_variance_decomposition.m
@@ -42,13 +42,13 @@ NumberOfDrawsFiles = rows(DrawsFiles);
 NumberOfSavedElementsPerSimulation = nvar*(nexo+1);
 MaXNumberOfDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation/8);
 
-if NumberOfSimulations<=MaXNumberOfDecompLines
-    Decomposition_array = zeros(NumberOfSimulations,nvar*(nexo+1));
+if SampleSize<=MaXNumberOfDecompLines
+    Decomposition_array = zeros(SampleSize,nvar*nexo);
     NumberOfDecompFiles = 1;
 else
-    Decomposition_array = zeros(MaXNumberOfDecompLines,nvar*(nexo+1));
-    NumberOfLinesInTheLastDecompFile = mod(NumberOfSimulations,MaXNumberOfDecompLines);
-    NumberOfDecompFiles = ceil(NumberOfSimulations/MaXNumberOfDecompLines);
+    Decomposition_array = zeros(MaXNumberOfDecompLines,nvar*nexo);
+    NumberOfLinesInTheLastDecompFile = mod(SampleSize,MaXNumberOfDecompLines);
+    NumberOfDecompFiles = ceil(SampleSize/MaXNumberOfDecompLines);
 end
 
 NumberOfDecompLines = rows(Decomposition_array);
@@ -59,7 +59,7 @@ DecompFileNumber = 1;
 % implied by each structural shock.
 linea = 0;
 for file = 1:NumberOfDrawsFiles
-    load([MhDirectoryName '/' DrawsFiles(file).name]);
+    load([M_.dname '/metropolis/' DrawsFiles(file).name]);
     isdrsaved = cols(pdraws)-1;
     NumberOfDraws = rows(pdraws);
     for linee = 1:NumberOfDraws
@@ -71,26 +71,22 @@ for file = 1:NumberOfDrawsFiles
             [dr,info] = resol(oo_.steady_state,0);
         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);
+                Decomposition_array(linea,(i-1)*nexo+j) = tmp{2}(i,j);
             end
         end
         if linea == NumberOfDecompLines
-            save([fname '_PosteriorVarianceDecomposition' int2str(DecompFileNumber)],'Decomposition_array');
+            save([M_.dname '/metropolis/' M_.fname '_PosteriorVarianceDecomposition' int2str(DecompFileNumber)],'Decomposition_array');
             DecompFileNumber = DecompFileNumber + 1;
             linea = 0;
             test = DecompFileNumber-NumberOfDecompFiles;
             if ~test% Prepare the last round...
-                Decomposition_array = zeros(NumberOfLinesInTheLastDecompFile,nvar*(nexo+1));
+                Decomposition_array = zeros(NumberOfLinesInTheLastDecompFile,nvar*nexo);
                 NumberOfDecompLines = NumberOfLinesInTheLastDecompFile;
                 DecompFileNumber = DecompFileNumber - 1;
             elseif test<0;
-                Decomposition_array = zeros(MaXNumberOfDecompLines,nvar*(nexo+1));
+                Decomposition_array = zeros(MaXNumberOfDecompLines,nvar*nexo);
             else
                 clear('Decomposition_array');
             end
@@ -98,10 +94,4 @@ for file = 1:NumberOfDrawsFiles
     end
 end
 
-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);
\ No newline at end of file
+options_.ar = nar;
\ No newline at end of file
diff --git a/matlab/posterior_analysis.m b/matlab/posterior_analysis.m
index a3ff246c445d79d4e1f44afa74966f90427f7df4..bcfa7f266dbdae9dbd6c0def0e4355ce4379a0be 100644
--- a/matlab/posterior_analysis.m
+++ b/matlab/posterior_analysis.m
@@ -1,15 +1,14 @@
-function posterior_analysis(type,arg1,arg2,options_,M_,oo_)  
+function oo_ = 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_);
-
+    SampleSize = options_.PosteriorSampleSize;
     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
@@ -21,117 +20,29 @@ function posterior_analysis(type,arg1,arg2,options_,M_,oo_)
           case 'variance'
             [nvar,vartan,NumberOfFiles] = ...
                 dsge_posterior_theoretical_covariance(SampleSize,M_,options_,oo_);
+            oo_ = covariance_posterior_analysis(SampleSize,M_.dname,M_.fname,...
+                                                vartan,nvar,arg1,arg2,options_.mh_conf_sig,oo_);          
           case 'decomposition'
             [nvar,vartan,NumberOfFiles] = ...
                 dsge_posterior_theoretical_variance_decomposition(SampleSize,M_,options_,oo_);
+            oo_ = variance_decomposition_posterior_analysis(SampleSize,M_.dname,M_.fname,...
+                                                            M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_);
           otherwise
             disp('Not yet implemented')
         end
       case 6
+        [ivar,vartan] = set_stationary_variables_list;
+        nvar = length(ivar);
         switch type
-          case 'variance'
-            covariance_posterior_analysis(NumberOfFiles,SampleSize,M_.dname,M_.fname,...
-                                          vartan,nvar,arg1,arg2);
+          case 'variance'    
+            oo_ = covariance_posterior_analysis(SampleSize,M_.dname,M_.fname,...
+                                          vartan,nvar,arg1,arg2,options_.mh_conf_sig,oo_);
           case 'decomposition'
-            variance_decomposition_posterior_analysis(NumberOfFiles,SampleSize,M_.dname,M_.fname,...
-                                                      M_.exo_names,arg2,vartan,nvar,arg1);
+            oo_ = variance_decomposition_posterior_analysis(SampleSize,M_.dname,M_.fname,...
+                                                      M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_);
           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
+    end
\ No newline at end of file
diff --git a/matlab/variance_decomposition_posterior_analysis.m b/matlab/variance_decomposition_posterior_analysis.m
new file mode 100644
index 0000000000000000000000000000000000000000..94426e64dc49e8ea04b28e073837f2e06d0a86f9
--- /dev/null
+++ b/matlab/variance_decomposition_posterior_analysis.m
@@ -0,0 +1,59 @@
+function oo_ = variance_decomposition_posterior_analysis(NumberOfSimulations,dname,fname, ...
+                                          exonames,exo,vartan,var,mh_conf_sig,oo_)
+    indx = 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
+    tmp = dir([ dname '/metropolis/'  fname '_PosteriorVarianceDecomposition*.mat']);
+    NumberOfFiles = length(tmp);
+    i1 = 1; tmp = zeros(NumberOfSimulations,1);
+    indice = (indx-1)*rows(exonames)+jndx;
+    for file = 1:NumberOfFiles
+        load([dname '/metropolis/' fname '_PosteriorVarianceDecomposition' int2str(file) '.mat']);
+        i2 = i1 + rows(Decomposition_array) - 1;
+        tmp(i1:i2) = Decomposition_array(:,indice);
+        i1 = i2+1;
+    end
+    name = [ var '.' exo ];
+    if isfield(oo_,'PosteriorTheoreticalMoments')
+        if isfield(oo_.PosteriorTheoreticalMoments,'dsge')
+            if isfield(oo_.PosteriorTheoreticalMoments.dsge,'VarianceDecomposition')
+                if isfield(oo_.PosteriorTheoreticalMoments.dsge.VarianceDecomposition.mean,name)
+                    % Nothing to do.
+                    return
+                end
+            end
+        end
+    end
+    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,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;']);
\ No newline at end of file