diff --git a/doc/dynare.texi b/doc/dynare.texi
index 40278494b9bcb5a8fcca56243596ca6be57b0247..67863a2766b4f2c1943ee16ee873c4cecace9459 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -3703,7 +3703,7 @@ period(s). The periods must be strictly positive. Conditional variances are give
 decomposition provides the decomposition of the effects of shocks upon
 impact. The results are stored in
 @code{oo_.conditional_variance_decomposition}
-(@pxref{oo_.conditional_variance_decomposition}). The variance decomposition is only conducted, if theoretical moments are requested, i.e. using the @code{periods=0}-option. In case of @code{order=2}, Dynare provides a second-order accurate approximation to the true second moments based on the linear terms of the second-order solution (see @cite{Kim, Kim, Schaumburg and Sims (2008)}).
+(@pxref{oo_.conditional_variance_decomposition}). The variance decomposition is only conducted, if theoretical moments are requested, i.e. using the @code{periods=0}-option. In case of @code{order=2}, Dynare provides a second-order accurate approximation to the true second moments based on the linear terms of the second-order solution (see @cite{Kim, Kim, Schaumburg and Sims (2008)}). Note that the unconditional variance decomposition (i.e. at horizon infinity) is automatically conducted if theoretical moments are requested (@pxref {oo_.variance_decomposition})
 
 @item pruning
 Discard higher order terms when iteratively computing simulations of
@@ -3852,7 +3852,7 @@ number of the matrix in the cell array corresponds to the order of
 autocorrelation. The option @code{ar} specifies the number of
 autocorrelation matrices available. Contains theoretical
 autocorrelations if the @code{periods} option is not present (or an approximation thereof for @code{order=2}), and
-empirical autocorrelations otherwise.
+empirical autocorrelations otherwise. The field is only created if stationary variables are present.
 
 The element @code{oo_.autocorr@{i@}(k,l)} is equal to the correlation
 between @math{y^k_t} and @math{y^l_{t-i}}, where @math{y^k}
@@ -3879,7 +3879,7 @@ details. Beware, this is the @i{autocorrelation} function, not the
 @i{autocovariance} function.
 
 @item oo_.gamma@{nar+2@}
-Variance decomposition.
+Unconditional variance decomposition @pxref{oo_.variance_decomposition}
 
 @item oo_.gamma@{nar+3@}
 If a second order approximation has been requested, contains the
@@ -3890,6 +3890,22 @@ In case of @code{order=2}, the theoretical second moments are a second order acc
 
 @end defvr
 
+@anchor{oo_.variance_decomposition}
+@defvr {MATLAB/Octave variable} oo_.variance_decomposition
+After a run of @code{stoch_simul} when requesting theoretical moments (@code{periods=0}), contains a matrix with the result of the unconditional variance decomposition (i.e. at horizon infinity). The first dimension corresponds to the endogenous variables (in the order of declaration) and the second dimension corresponds to exogenous variables (in the order of declaration). Numbers are in percent and sum up to 100 across columns.
+@end defvr
+
+@anchor{oo_.conditional_variance_decomposition}
+@defvr {MATLAB/Octave variable} oo_.conditional_variance_decomposition
+After a run of @code{stoch_simul} with the
+@code{conditional_variance_decomposition} option, contains a
+three-dimensional array with the result of the decomposition. The
+first dimension corresponds to forecast horizons (as declared with the
+option), the second dimension corresponds to endogenous variables (in
+the order of declaration), the third dimension corresponds to
+exogenous variables (in the order of declaration).
+@end defvr
+
 @defvr {MATLAB/Octave variable} oo_.irfs
 After a run of @code{stoch_simul} with option @code{irf} different
 from zero, contains the impulse responses, with the following naming
@@ -4183,16 +4199,6 @@ three times in the unfolded @math{G_3} matrix, they must be multiplied
 by 3 when computing the decision rules.
 @end itemize
 
-@anchor{oo_.conditional_variance_decomposition}
-@defvr {MATLAB/Octave variable} oo_.conditional_variance_decomposition
-After a run of @code{stoch_simul} with the
-@code{conditional_variance_decomposition} option, contains a
-three-dimensional array with the result of the decomposition. The
-first dimension corresponds to forecast horizons (as declared with the
-option), the second dimension corresponds to endogenous variables (in
-the order of declaration), the third dimension corresponds to
-exogenous variables (in the order of declaration).
-@end defvr
 
 @node Estimation
 @section Estimation
@@ -5044,8 +5050,7 @@ model. Default: @code{4}.
 @anchor{moments_varendo} Triggers the computation of the posterior
 distribution of the theoretical moments of the endogenous
 variables. Results are stored in
-@code{oo_.PosteriorTheoreticalMoments} (see below for a description of
-this variable). The number of lags in the autocorrelation function is
+@code{oo_.PosteriorTheoreticalMoments} (@pxref{oo_.PosteriorTheoreticalMoments}). The number of lags in the autocorrelation function is
 controlled by the @code{ar} option.
 
 @item conditional_variance_decomposition = @var{INTEGER}
@@ -5062,7 +5067,7 @@ period 1, the conditional variance decomposition provides the
 decomposition of the effects of shocks upon impact. The results are
 stored in
 @code{oo_.PosteriorTheoreticalMoments.dsge.ConditionalVarianceDecomposition},
-but currently there is no output. Note that this option requires the
+but currently there is no displayed output. Note that this option requires the
 option @code{moments_varendo} to be specified.
 
 @item filtered_vars
@@ -5489,6 +5494,7 @@ After an estimation with Metropolis, fields are of the form:
 @end defvr
 
 @defvr {MATLAB/Octave variable} oo_.PosteriorTheoreticalMoments
+@anchor{oo_.PosteriorTheoreticalMoments}
 Variable set by the @code{estimation} command, if it is used with the
 @code{moments_varendo} option. Fields are of the form:
 @example
@@ -5506,7 +5512,7 @@ Auto- and cross-correlation of endogenous variables. Fields are vectors with cor
 
 
 @item VarianceDecomposition
-Decomposition of variance@footnote{When the shocks are correlated, it
+Decomposition of variance (unconditional variance, i.e. at horizon infinity)@footnote{When the shocks are correlated, it
 is the decomposition of orthogonalized shocks via Cholesky
 decomposition according to the order of declaration of shocks
 (@pxref{Variable declarations})}
diff --git a/matlab/disp_moments.m b/matlab/disp_moments.m
index 058bebc97c7c248816c25728bc77da5547ec9cb5..478936e9659af94efb2337e66e94b8ec3859e1d0 100644
--- a/matlab/disp_moments.m
+++ b/matlab/disp_moments.m
@@ -80,6 +80,10 @@ if options_.nocorr == 0
     end
 end
 
+if options_.noprint == 0 && length(options_.conditional_variance_decomposition)
+   fprintf('\nSTOCH_SIMUL: conditional_variance_decomposition requires theoretical moments, i.e. periods=0.\n') 
+end
+
 ar = options_.ar;
 if ar > 0
     autocorr = [];
diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m
index d1a4876a296cd575dbf51016093277acdcb45767..0c1a245a6db5874bd0d500ef703aa1fa72a2993f 100644
--- a/matlab/disp_th_moments.m
+++ b/matlab/disp_th_moments.m
@@ -50,43 +50,54 @@ z = [ m sd s2 ];
 oo_.mean = m;
 oo_.var = oo_.gamma_y{1};
 
-if ~options_.noprint %options_.nomoments == 0
-    if options_.order == 2
-        title='APROXIMATED THEORETICAL MOMENTS';
-    else
-        title='THEORETICAL MOMENTS';
-    end
-    if options_.hp_filter
-        title = [title ' (HP filter, lambda = ' num2str(options_.hp_filter) ')'];
-    end
-    headers=char('VARIABLE','MEAN','STD. DEV.','VARIANCE');
-    labels = deblank(M_.endo_names(ivar,:));
-    lh = size(labels,2)+2;
-    dyntable(title,headers,labels,z,lh,11,4);
-    if M_.exo_nbr > 1 && size(stationary_vars, 1) > 0
-        skipline()
+if M_.exo_nbr > 1 && size(stationary_vars, 1) > 0
+    oo_.variance_decomposition=100*oo_.gamma_y{options_.ar+2};
+    if ~options_.noprint %options_.nomoments == 0
         if options_.order == 2
-            title='APPROXIMATED VARIANCE DECOMPOSITION (in percent)';            
+            title='APROXIMATED THEORETICAL MOMENTS';
         else
-            title='VARIANCE DECOMPOSITION (in percent)';
+            title='THEORETICAL MOMENTS';
         end
         if options_.hp_filter
-            title = [title ' (HP filter, lambda = ' ...
-                     num2str(options_.hp_filter) ')'];
+            title = [title ' (HP filter, lambda = ' num2str(options_.hp_filter) ')'];
         end
-        headers = M_.exo_names;
-        headers(M_.exo_names_orig_ord,:) = headers;
-        headers = char(' ',headers);
-        lh = size(deblank(M_.endo_names(ivar(stationary_vars),:)),2)+2;
-        dyntable(title,headers,deblank(M_.endo_names(ivar(stationary_vars), ...
-                                                     :)),100*oo_.gamma_y{options_.ar+2}(stationary_vars,:),lh,8,2);
+        headers=char('VARIABLE','MEAN','STD. DEV.','VARIANCE');
+        labels = deblank(M_.endo_names(ivar,:));
+        lh = size(labels,2)+2;
+        dyntable(title,headers,labels,z,lh,11,4);
+
+        skipline()
+            if options_.order == 2
+                title='APPROXIMATED VARIANCE DECOMPOSITION (in percent)';            
+            else
+                title='VARIANCE DECOMPOSITION (in percent)';
+            end
+            if options_.hp_filter
+                title = [title ' (HP filter, lambda = ' ...
+                         num2str(options_.hp_filter) ')'];
+            end
+            headers = M_.exo_names;
+            headers(M_.exo_names_orig_ord,:) = headers;
+            headers = char(' ',headers);
+            lh = size(deblank(M_.endo_names(ivar(stationary_vars),:)),2)+2;
+            dyntable(title,headers,deblank(M_.endo_names(ivar(stationary_vars), ...
+                                                         :)),100*oo_.gamma_y{options_.ar+2}(stationary_vars,:),lh,8,2);
     end
     
     conditional_variance_steps = options_.conditional_variance_decomposition;
     if length(conditional_variance_steps)
-        oo_ = display_conditional_variance_decomposition(conditional_variance_steps,...
-                                                         ivar,dr,M_, ...
-                                                         options_,oo_);
+        StateSpaceModel.number_of_state_equations = M_.endo_nbr;
+        StateSpaceModel.number_of_state_innovations = M_.exo_nbr;
+        StateSpaceModel.sigma_e_is_diagonal = M_.sigma_e_is_diagonal;
+        [StateSpaceModel.transition_matrix,StateSpaceModel.impulse_matrix] = kalman_transition_matrix(dr,(1:M_.endo_nbr)',M_.nstatic+(1:M_.nspred)',M_.exo_nbr);
+        StateSpaceModel.state_innovations_covariance_matrix = M_.Sigma_e;
+        StateSpaceModel.order_var = dr.order_var;
+        oo_.conditional_variance_decomposition = conditional_variance_decomposition(StateSpaceModel,conditional_variance_steps,ivar);
+        
+        if options_.noprint == 0
+            display_conditional_variance_decomposition(oo_.conditional_variance_decomposition,conditional_variance_steps,...
+                                                         ivar,M_,options_);
+        end
     end
 end
 
diff --git a/matlab/display_conditional_variance_decomposition.m b/matlab/display_conditional_variance_decomposition.m
index 7ce71478c2fceb1a342e3fd41c84c1b860d1fdd6..2e15a1d9e08021fd14759b8586ea5193097c2220 100644
--- a/matlab/display_conditional_variance_decomposition.m
+++ b/matlab/display_conditional_variance_decomposition.m
@@ -1,22 +1,19 @@
-function oo_ = display_conditional_variance_decomposition(Steps, SubsetOfVariables, dr,M_,options_,oo_)
-% This function computes the conditional variance decomposition of a given state space model
+function display_conditional_variance_decomposition(conditional_decomposition_array,Steps,SubsetOfVariables,M_,options_)
+% This function displays the conditional variance decomposition of a given state space model
 % for a subset of endogenous variables.
 % 
 % INPUTS 
-%   StateSpaceModel     [structure]   Specification of the state space model.
+%   conditional_decomposition_array     [matrix]   Output matrix from compute_conditional_variance_decomposition
 %   Steps               [integer]     1*h vector of dates.
 %   SubsetOfVariables   [integer]     1*q vector of indices.
-%    
+%   M_                  [structure]   Dynare structure containing the
+%                                     Model description
+%   options_            [structure]   Dynare structure containing the
+%                                     options
 % OUTPUTS 
-%   PackedConditionalVarianceDecomposition  [double] n(n+1)/2*p matrix, where p is the number of state innovations and
-%                                                    n is equal to length(SubsetOfVariables).    
+%   none
 %
-% SPECIAL REQUIREMENTS
-%
-% [1] The covariance matrix of the state innovations needs to be diagonal.
-% [2] In this version, absence of measurement errors is assumed...
-
-% Copyright (C) 2010-2013 Dynare Team
+% Copyright (C) 2010-2014 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -33,49 +30,28 @@ function oo_ = display_conditional_variance_decomposition(Steps, SubsetOfVariabl
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-endo_nbr = M_.endo_nbr;
-exo_nbr = M_.exo_nbr;
-StateSpaceModel.number_of_state_equations = M_.endo_nbr;
-StateSpaceModel.number_of_state_innovations = exo_nbr;
-StateSpaceModel.sigma_e_is_diagonal = M_.sigma_e_is_diagonal;
-
-iv = (1:endo_nbr)';
-ic = M_.nstatic+(1:M_.nspred)';
-
-[StateSpaceModel.transition_matrix,StateSpaceModel.impulse_matrix] = kalman_transition_matrix(dr,iv,ic,exo_nbr);
-StateSpaceModel.state_innovations_covariance_matrix = M_.Sigma_e;
-StateSpaceModel.order_var = dr.order_var;
-
-conditional_decomposition_array = conditional_variance_decomposition(StateSpaceModel,Steps,SubsetOfVariables );
-
-if options_.noprint == 0
-  if options_.order == 2
-    skipline()                
-    disp('APPROXIMATED CONDITIONAL VARIANCE DECOMPOSITION (in percent)')
-  else
-    skipline()               
-    disp('CONDITIONAL VARIANCE DECOMPOSITION (in percent)')
-  end
+if options_.order == 2
+skipline()                
+disp('APPROXIMATED CONDITIONAL VARIANCE DECOMPOSITION (in percent)')
+else
+skipline()               
+disp('CONDITIONAL VARIANCE DECOMPOSITION (in percent)')
 end
 
-vardec_i = zeros(length(SubsetOfVariables),exo_nbr);
+vardec_i = zeros(length(SubsetOfVariables),M_.exo_nbr);
 
 for i=1:length(Steps)
     disp(['Period ' int2str(Steps(i)) ':'])
-    
-    for j=1:exo_nbr
+
+    for j=1:M_.exo_nbr
         vardec_i(:,j) = 100*conditional_decomposition_array(:, ...
                                                           i,j);
     end
-    if options_.noprint == 0
-        headers = M_.exo_names;
-        headers(M_.exo_names_orig_ord,:) = headers;
-        headers = char(' ',headers);
-        lh = size(deblank(M_.endo_names(SubsetOfVariables,:)),2)+2;
-        dyntable('',headers,...
-                 deblank(M_.endo_names(SubsetOfVariables,:)),...
-                 vardec_i,lh,8,2);
-    end
-end
-
-oo_.conditional_variance_decomposition = conditional_decomposition_array;
+    headers = M_.exo_names;
+    headers(M_.exo_names_orig_ord,:) = headers;
+    headers = char(' ',headers);
+    lh = size(deblank(M_.endo_names(SubsetOfVariables,:)),2)+2;
+    dyntable('',headers,...
+             deblank(M_.endo_names(SubsetOfVariables,:)),...
+             vardec_i,lh,8,2);
+end
\ No newline at end of file