diff --git a/matlab/epilogue_shock_decomposition.m b/matlab/epilogue_shock_decomposition.m
index b610c5f35e305079ccbab7560cefed881a2bd13f..4a76cd2defcf2ac9e832d3225dabfab7dd1d440d 100644
--- a/matlab/epilogue_shock_decomposition.m
+++ b/matlab/epilogue_shock_decomposition.m
@@ -27,7 +27,11 @@ nterms = size(z,2);
 for k=1:size(y,1)
     ytmp  = squeeze(y(k,:,:));
     yres = ytmp(end,:) - sum(ytmp(1:end-1,:));
-    w = abs(ytmp(1:end-1,:))./sum(abs(ytmp(1:end-1,:)));
+    if ~isoctave && matlab_ver_less_than('9.1') % Automatic broadcasting was introduced in MATLAB R2016b
+        w = bsxfun(@rdivide, abs(ytmp(1:end-1,:)), sum(abs(ytmp(1:end-1,:))))
+    else
+        w = abs(ytmp(1:end-1,:))./sum(abs(ytmp(1:end-1,:)));
+    end
     %     ytmp(1:end-1,:) = ytmp(1:end-1,:) + repmat(yres,[nterms-1 1])/(nterms-1);
     ytmp(1:end-1,:) = ytmp(1:end-1,:) + repmat(yres,[nterms-1 1]).*w;
     y(k,:,:) = ytmp;
diff --git a/matlab/missing/stats/+gamrnd/best_1978.m b/matlab/missing/stats/+gamrnd/best_1978.m
index e2820a4d914b3fdef7527e9eae52abd016e83ef6..1934c71cb37338379e45808d4a7985e332320aea 100644
--- a/matlab/missing/stats/+gamrnd/best_1978.m
+++ b/matlab/missing/stats/+gamrnd/best_1978.m
@@ -9,7 +9,7 @@ function  g = best_1978(a ,b)
 % OUTPUTS
 % - g    [double]     n*1 vector, gamma variates.
 
-% Copyright (C) 2006-2018 Dynare Team
+% Copyright (C) 2006-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -44,7 +44,12 @@ while mm
     X(index) = bb(index)+Y(index); % x
     id1 = index(X(index)<0); % Reject.
     id2 = setdiff(index, id1);
-    Z(id2) = 64.0*(W(id2).^3).*(rand(length(id2),1).^2); % d
+    if numel(id2) ~= 0 || isoctave || ~matlab_ver_less_than('9.1')
+        % If id2=[], LHS of the .* has size [0,0], while RHS is [0,1].
+        % Since there is no automatic broadcast in MATLAB < R2016b, skip the
+        % statement in that case.
+        Z(id2) = 64.0*(W(id2).^3).*(rand(length(id2),1).^2); % d
+    end
     id3 = id2(Z(id2)>1.0-2.0*Y(id2).*Y(id2)./X(id2)); % Reject.
     id4 = id3(log(Z(id3))>2.0*(bb(id3).*log(X(id3)./bb(id3))-Y(id3))); % Reject.
     index = [id1, id4];
diff --git a/matlab/missing/stats/+gamrnd/knuth.m b/matlab/missing/stats/+gamrnd/knuth.m
index 02a3e25592d587cdf7f168598c72b0c92839e2e5..ed4ca4b518525cf07583650a8a6a91d12c6fcd9a 100644
--- a/matlab/missing/stats/+gamrnd/knuth.m
+++ b/matlab/missing/stats/+gamrnd/knuth.m
@@ -9,7 +9,7 @@ function  g = knuth(a, b)
 % OUTPUTS
 % - g    [double]     n*1 vector, gamma variates.
 
-% Copyright (C) 2006-2018 Dynare Team
+% Copyright (C) 2006-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -38,9 +38,17 @@ while mm
     X(index) = Y(index).*bb(index) + a(index) - 1;
     id1 = index(X(index)<=0); % Rejected draws.
     id2 = setdiff(index, id1);
-    id3 = id2(rand(length(id2), 1)>(1+Y(id2).*Y(id2)).*exp((a(id2)-1).*(log(X(id2))-log(a(id2)-1))-bb(id2).*Y(id2))); % Rejected draws.
+    if numel(id2) == 0 && ~isoctave && matlab_ver_less_than('9.1')
+        % The LHS of the > comparison in the "else" branch has size = [0, 1]
+        % while the RHS has size = [0, 0].
+        % Automatic broadcasting was only introduced in MATLAB R2016b, so we
+        % must handle this case separately to avoid an error.
+        id3 = [];
+    else
+        id3 = id2(rand(length(id2), 1)>(1+Y(id2).*Y(id2)).*exp((a(id2)-1).*(log(X(id2))-log(a(id2)-1))-bb(id2).*Y(id2))); % Rejected draws.
+    end
     index = [id1, id3];
     mm = length(index);
 end
 
-g = X.*b;
\ No newline at end of file
+g = X.*b;