diff --git a/src/@dseries/nanmean.m b/src/@dseries/nanmean.m
index 3a0ca58aa7249498d9c3f304941af85be0949ef8..a77e2e7b4f41408a69f3683f0d58fdf2a2448edb 100644
--- a/src/@dseries/nanmean.m
+++ b/src/@dseries/nanmean.m
@@ -1,15 +1,16 @@
-function m = nanmean(o) % --*-- Unitary tests --*--
+function m = nanmean(o, geometric) % --*-- Unitary tests --*--
 
-% Returns the mean of the variables in a @dseries object o.
+% Returns the mean of the variables in a @dseries object o (robust
+% to the presence of NaNs).
 %
 % INPUTS
-%  o o             dseries object [mandatory].
-%  o geometric     logical [default is false], if true returns the geometric mean.
+% - o             dseries object [mandatory].
+% - geometric     logical [default is false], if true returns the geometric mean.
 %
 % OUTPUTS
-%  o m             1*vobs(o) vector of doubles.
+% - m             1*vobs(o) vector of doubles.
 
-% Copyright (C) 2019 Dynare Team
+% Copyright © 2019-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -26,26 +27,59 @@ function m = nanmean(o) % --*-- Unitary tests --*--
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-m = nanmean(o.data);
+if nargin<2
+    geometric = false;
+end
 
+if geometric
+    m = NaN(1,o.vobs());
+    for i = 1:o.vobs()
+        tmp = o.data(~isnan(o.data(:,i)),i);
+        m(i) = prod(tmp)^(1.0/length(tmp));
+    end
+else
+    m = nanmean(o.data);
+end
+
+return
 %@test:1
-%$ % Define a dataset.
-%$ A = repmat([1.005, 1.05], 10, 1);
-%$ A(3,1) = NaN;
-%$ A(5,2) = NaN;
-%$
-%$ % Instantiate a time series object and compute the mean.
-%$ try
-%$    ts = dseries(A);
-%$    m = nanmean(ts);
-%$    t(1) = 1;
-%$ catch
-%$    t = 0;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = dassert(isequal(size(m),[1, 2]), true);
-%$    t(3) = dassert(m, [1.005, 1.05]);
-%$ end
-%$ T = all(t);
-%@eof:1
\ No newline at end of file
+% Define a dataset.
+A = repmat([1.005, 1.05], 10, 1);
+A(3,1) = NaN;
+A(5,2) = NaN;
+
+% Instantiate a time series object and compute the mean.
+try
+   ts = dseries(A);
+   m = nanmean(ts);
+   t(1) = 1;
+catch
+   t = 0;
+end
+
+if t(1)
+   t(2) = dassert(isequal(size(m),[1, 2]), true);
+   t(3) = dassert(m, [1.005, 1.05]);
+end
+T = all(t);
+%@eof:1
+
+%@test:2
+% Define a dataset.
+a = [1 0; NaN 2; 3 4];
+
+% Instantiate a time series object and compute the geometric mean.
+try
+   ts = dseries(a);
+   m = nanmean(ts, true);
+   t(1) = 1;
+catch
+   t = 0;
+end
+
+if t(1)
+   t(2) = dassert(isequal(size(m),[1, 2]), true);
+   t(3) = dassert(m, [sqrt(3), 0]);
+end
+T = all(t);
+%@eof:2
\ No newline at end of file