diff --git a/doc/dynare.texi b/doc/dynare.texi
index 3647a824e2f24609f606f946ad8cc1cc2c33ef82..279a512371196ef99f3b733a825dccc077884ba4 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -3161,7 +3161,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. Currently, variance decompositions are only implemented for @code{order=1}. In case of @code{order=2}, Dynare will thus not display the variance decomposition based on a second order approximation, but on a first order approximation. 
+(@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)}).  
 
 @item pruning
 Discard higher order terms when iteratively computing simulations of
@@ -3292,7 +3292,7 @@ in declaration order.
 @defvr {MATLAB/Octave variable} oo_.var
 After a run of @code{stoch_simul}, contains the variance-covariance of
 the endogenous variables. Contains theoretical variance if the
-@code{periods} option is not present, and empirical variance
+@code{periods} option is not present (or an approximation thereof for @code{order=2}), and empirical variance
 otherwise. The variables are arranged in declaration order.
 @end defvr
 
@@ -3303,7 +3303,7 @@ autocorrelation matrices of the endogenous variables. The element
 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, and
+autocorrelations if the @code{periods} option is not present (or an approximation thereof for @code{order=2}), and
 empirical autocorrelations otherwise.
 
 The element @code{oo_.autocorr@{i@}(k,l)} is equal to the correlation
@@ -3338,6 +3338,8 @@ If a second order approximation has been requested, contains the
 vector of the mean correction terms.
 @end table
 
+In case of @code{order=2}, the theoretical second moments are a second order accurate approximation of the true second moments, see @code{conditional_variance_decomposition}. 
+
 @end defvr
 
 @defvr {MATLAB/Octave variable} oo_.irfs
diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m
index f61f1e17aaf8aff70609470d1c3250c91d34f71e..cf3397a646bd5406ceb1be376f00a55746ad0f92 100644
--- a/matlab/disp_th_moments.m
+++ b/matlab/disp_th_moments.m
@@ -51,7 +51,11 @@ oo_.mean = m;
 oo_.var = oo_.gamma_y{1};
 
 if ~options_.noprint %options_.nomoments == 0
-    title='THEORETICAL MOMENTS';
+    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
@@ -62,7 +66,7 @@ if ~options_.noprint %options_.nomoments == 0
     if M_.exo_nbr > 1 && size(stationary_vars, 1) > 0
         disp(' ')
         if options_.order == 2
-            title='VARIANCE DECOMPOSITION (in percent), based on first order approximation';            
+            title='APPROXIMATED VARIANCE DECOMPOSITION (in percent)';            
         else
             title='VARIANCE DECOMPOSITION (in percent)';
         end
@@ -97,7 +101,11 @@ if options_.nocorr == 0 && size(stationary_vars, 1) > 0
     corr = oo_.gamma_y{1}(i1,i1)./(sd(i1)*sd(i1)');
     if ~options_.noprint,
         disp(' ')
-        title='MATRIX OF CORRELATIONS';
+        if options_.order == 2
+            title='APPROXIMATED MATRIX OF CORRELATIONS';            
+        else
+            title='MATRIX OF CORRELATIONS';
+        end
         if options_.hp_filter
             title = [title ' (HP filter, lambda = ' num2str(options_.hp_filter) ')'];
         end
@@ -115,7 +123,11 @@ if options_.ar > 0 && size(stationary_vars, 1) > 0
     end
     if ~options_.noprint,      
         disp(' ')    
-        title='COEFFICIENTS OF AUTOCORRELATION';      
+        if options_.order == 2
+            title='APPROXIMATED COEFFICIENTS OF AUTOCORRELATION';            
+        else
+            title='COEFFICIENTS OF AUTOCORRELATION';
+        end
         if options_.hp_filter        
             title = [title ' (HP filter, lambda = ' num2str(options_.hp_filter) ')'];      
         end      
diff --git a/matlab/display_conditional_variance_decomposition.m b/matlab/display_conditional_variance_decomposition.m
index 97811025a610e87f3141686562bdfb674620faad..ed4f334fcea4de3d708f10b288bb8fec6e46871c 100644
--- a/matlab/display_conditional_variance_decomposition.m
+++ b/matlab/display_conditional_variance_decomposition.m
@@ -51,7 +51,7 @@ conditional_decomposition_array = conditional_variance_decomposition(StateSpaceM
 if options_.noprint == 0
   if options_.order == 2
     disp(' ')                
-    disp('CONDITIONAL VARIANCE DECOMPOSITION (in percent), based on first order approximation')
+    disp('APPROXIMATED CONDITIONAL VARIANCE DECOMPOSITION (in percent)')
   else
     disp(' ')                
     disp('CONDITIONAL VARIANCE DECOMPOSITION (in percent)')
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index 23bfdcbd2b82a55b7ef66f0f4998e11a4b8608cf..694c1655e934c2262c6beedcff9eee3841dfba85 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -761,8 +761,12 @@ else
     lnprior = priordens(xparam1,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
 end
 if DynareOptions.endogenous_prior==1
-  [lnpriormom]  = endogenous_prior(Y,Pstar,BayesInfo);
-  fval    = (likelihood-lnprior-lnpriormom);
+  if DynareOptions.lik_init==2 || DynareOptions.lik_init==3
+    error('Endogenous prior not supported with non-stationary models')
+  else
+    [lnpriormom]  = endogenous_prior(Y,Pstar,BayesInfo,H);
+    fval    = (likelihood-lnprior-lnpriormom);
+  end
 else
   fval    = (likelihood-lnprior);
 end
diff --git a/matlab/endogenous_prior.m b/matlab/endogenous_prior.m
index b9b7f54aba61aa65c946837fc2bc0196bbbd5f8f..c2930abad8a9b82e316df4ce2c4e26483bb0c85f 100644
--- a/matlab/endogenous_prior.m
+++ b/matlab/endogenous_prior.m
@@ -1,4 +1,4 @@
-function [lnpriormom] = endogenous_prior(data,Pstar,BayesInfo)
+function [lnpriormom] = endogenous_prior(data,Pstar,BayesInfo,H)
 % Computes the endogenous log prior addition to the initial prior
 %
 % INPUTS
@@ -82,8 +82,7 @@ mf=BayesInfo.mf1;
 II=eye(size(Pstar,2));
 Z=II(mf,:);
 % This is Ftheta, variance of model variables, given param vector theta:
-Ftheta=diag(Z*Pstar(:,mf));  % +H;
-
+Ftheta=diag(Z*Pstar(:,mf)+H);
 % below commented out line is for Del Negro Schorfheide style priors:
 %     lnpriormom=-.5*n*TT*log(2*pi)-.5*TT*log(det(sigma))-.5*TT*trace(inv(sigma)*(gamyy-2*phi'*gamxy+phi'*gamxx*phi));
 lnpriormom=.5*n*log(Tsamp/(2*pi))-.5*log(det(Shat))-.5*Tsamp*(Fhat-Ftheta)'/Shat*(Fhat-Ftheta);
diff --git a/matlab/rplot.m b/matlab/rplot.m
index 101a188a5144d7b3ce7998186fe4018ed05224e8..164d26d172512ed2fd4d1776ca3b2d2a40f150b7 100644
--- a/matlab/rplot.m
+++ b/matlab/rplot.m
@@ -58,7 +58,7 @@ if rplottype == 0
     for j = 1:size(y,1)
         t = [t s1(j,:) ' '] ;
     end
-    figure ;
+    figure('Name',['Simulated Trajectory']) ;
     plot(ix(i),y(:,i)) ;
     title (t,'Interpreter','none') ;
     xlabel('Periods') ;
@@ -72,13 +72,13 @@ if rplottype == 0
     end
 elseif rplottype == 1
     for j = 1:size(y,1)
-        figure ;
+        figure('Name',['Simulated Trajectory']) ;
         plot(ix(i),y(j,i)) ;
         title(['Plot of ' s1(j,:)],'Interpreter','none') ;
         xlabel('Periods') ;
     end
 elseif rplottype == 2
-    figure ;
+    figure('Name',['Simulated Trajectory']) ;
     nl = max(1,fix(size(y,1)/4)) ;
     nc = ceil(size(y,1)/nl) ;
     for j = 1:size(y,1)
@@ -96,6 +96,7 @@ end
 % 02/28/01 MJ replaced bseastr by MATLAB's strmatch
 % 06/19/01 MJ added 'exact' to strmatch calls
 % 06/25/03 MJ correction when options_.smpl ~= 0
+% 03/18/13 JP bugfix for rplottype>0; added figure names
 
 
 
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index ec48719e0c436f9e80d8ad7c20ab1fe6bb9727ce..e06b413bb99fa78ce94e7d320a24d97ef5e54ac2 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -122,7 +122,7 @@ end
 if options_.periods > 0 && ~PI_PCL_solver
     if options_.periods <= options_.drop
         disp(['STOCH_SIMUL error: The horizon of simulation is shorter' ...
-              ' than the number of observations to be DROPed'])
+              ' than the number of observations to be dropped'])
         options_ =options_old;
         return
     end
@@ -142,9 +142,6 @@ if options_.nomoments == 0
     elseif options_.periods == 0
         % There is no code for theoretical moments at 3rd order
         if options_.order <= 2
-        if options_.order == 2
-            warning('You have requested a second order approximation, but variance decompositions currently only allow for first order. Displaying decompositions at order=1 instead.')
-        end
             disp_th_moments(oo_.dr,var_list);
         end
     else