diff --git a/doc/bvar-a-la-sims.tex b/doc/bvar-a-la-sims.tex
index 637bc5cad..9296485b8 100644
--- a/doc/bvar-a-la-sims.tex
+++ b/doc/bvar-a-la-sims.tex
@@ -169,7 +169,7 @@ $$\left[
 $$\left[
 \begin{array}{cc}
 0 & 0 \\
-0 & 0 
+0 & 0
 \end{array}
 \right]
 =
@@ -499,7 +499,7 @@ The syntax for computing the marginal density is:
 
 The options are those described above.
 
-The command will actually compute the marginal density for several models: first for the model with one lag, then with two lags, and so on up to \textit{max\_number\_of\_lags} lags.
+The command will actually compute the marginal density for several models: first for the model with one lag, then with two lags, and so on up to \textit{max\_number\_of\_lags} lags. Results will be stored in a \textit{max\_number\_of\_lags} by 1 vector \texttt{oo\_.bvar.log\_marginal\_data\_density}.
 
 \subsection{Forecasting}
 
diff --git a/matlab/bvar_density.m b/matlab/bvar_density.m
index 665a9dae5..db6271082 100644
--- a/matlab/bvar_density.m
+++ b/matlab/bvar_density.m
@@ -29,6 +29,10 @@ function bvar_density(maxnlags)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
+global oo_
+
+oo_.bvar.log_marginal_data_density=NaN(maxnlags,1);
+
 for nlags = 1:maxnlags
     [ny, nx, posterior, prior] = bvar_toolbox(nlags);
     
@@ -39,6 +43,8 @@ for nlags = 1:maxnlags
     
     log_dnsty = posterior_int - prior_int - 0.5*ny*lik_nobs*log(2*pi);
     
+    oo_.bvar.log_marginal_data_density(nlags)=log_dnsty;
+    
     skipline()
     fprintf('The marginal log density of the BVAR(%g) model is equal to %10.4f\n', ...
             nlags, log_dnsty);