diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index db05b8538238d1335741bb7a6b2f32008862969f..5da7340488986e31ea4465bf3d13f9bdd2b9f738 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -101,10 +101,10 @@ It's useful to contribute `.mod` files that test some aspect of Dynare that is n
 
 ### Unit tests
 
-So-called unit tests allow the test suite to check the correct functioning of the MATLAB/Octave functions contained in Dynare. To add a unit test you need to 
-1. add the keyword ` % --*-- Unitary tests --*--` at the end of the `function` header to tell the testsuite that the file contains unit tests. 
-1. Add the particular tests at the end of the file after a `return` by
-   1. Starting a test with `%@test:INTEGER` 
+So-called unit tests allow the test suite to check the correct functioning of the MATLAB/Octave functions contained in Dynare. To add a unit test you need to
+1. add the `return % --*-- Unit tests --*--` at the end of the `function` to tell the testsuite that the file contains unit tests.
+1. Add the particular tests at the end of the file after the `return` statement by
+   1. Starting a test with `%@test:INTEGER`
    2. Adding a MATLAB/Octave test code that provides a pass/fail indicator `T` that takes on `true` if the test passed.
    3. Closing the test with `%@eof:INTEGER`
    where `INTEGER` denotes the number of the test.
@@ -112,9 +112,9 @@ So-called unit tests allow the test suite to check the correct functioning of th
 An example testing the correct functionality of mode-computations for a normal distribution is
 
 ```
-function m = compute_prior_mode(hyperparameters,shape) % --*-- Unitary tests --*--
+function m = compute_prior_mode(hyperparameters,shape)
 
-return
+return  % --*-- Unit tests --*--
 
 %@test:1
 % Normal density
diff --git a/matlab/cubature_with_gaussian_weight.m b/matlab/cubature_with_gaussian_weight.m
index 48b970ee75a544049e47cd3787ac6f2b3ec43898..20642b0cd8ca8426b0ac85a0010400947b776c3c 100644
--- a/matlab/cubature_with_gaussian_weight.m
+++ b/matlab/cubature_with_gaussian_weight.m
@@ -1,4 +1,4 @@
-function [nodes, weights] = cubature_with_gaussian_weight(d,n,method)  % --*-- Unitary tests --*--
+function [nodes, weights] = cubature_with_gaussian_weight(d, n, method)
 
 % Computes nodes and weights for a n-order cubature with gaussian weight.
 %
@@ -17,7 +17,7 @@ function [nodes, weights] = cubature_with_gaussian_weight(d,n,method)  % --*-- U
 %     ∫   f(x) × e        dx
 %      -∞
 
-% Copyright © 2012-2019 Dynare Team
+% Copyright © 2012-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -108,7 +108,7 @@ m(:,2) =  e(n,i)-e(n,j);
 m(:,3) = -m(:,2);
 m(:,4) = -m(:,1);
 
-return
+return % --*-- Unit tests --*--
 
 %@test:1
 d = 4;
@@ -267,4 +267,4 @@ if t(1)
     t(4) = dassert(m3,zeros(d,1),1e-12);
     T = all(t);
 end
-%@eof:6
\ No newline at end of file
+%@eof:6
diff --git a/matlab/distributions/beta_specification.m b/matlab/distributions/beta_specification.m
index 097880ec178f7e48530a5811ea4a1a162451463b..e6dc77d9031821b018bb577dcf12e44f04b81b4b 100644
--- a/matlab/distributions/beta_specification.m
+++ b/matlab/distributions/beta_specification.m
@@ -1,4 +1,4 @@
-function [a, b] = beta_specification(mu, sigma2, lb, ub, name)   % --*-- Unitary tests --*--
+function [a, b] = beta_specification(mu, sigma2, lb, ub, name)
 
 % Returns the hyperparameters of the beta distribution given the expectation and variance.
 %
@@ -12,7 +12,7 @@ function [a, b] = beta_specification(mu, sigma2, lb, ub, name)   % --*-- Unitary
 % - a      [double]   First hyperparameter of the Beta density.
 % - b      [double]   Second hyperparameter of the Beta density.
 
-% Copyright © 2015-2017 Dynare Team
+% Copyright © 2015-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -66,67 +66,69 @@ end
 a = (1-mu)*mu*mu/sigma2-mu;
 b = a*(1/mu-1);
 
+return % --*-- Unit tests --*--
+
 %@test:1
-%$ try
-%$    [a, b] = beta_specification(.5, .05);
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = abs(0.5-a/(a+b))<1e-12;
-%$    t(3) = abs(0.05-a*b/((a+b)^2*(a+b+1)))<1e-12;
-%$ end
-%$ T = all(t);
+try
+   [a, b] = beta_specification(.5, .05);
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = abs(0.5-a/(a+b))<1e-12;
+    t(3) = abs(0.05-a*b/((a+b)^2*(a+b+1)))<1e-12;
+end
+T = all(t);
 %@eof:1
 
 %@test:2
-%$ try
-%$    [a, b] = beta_specification(0.5, .05, 1, 2);
-%$    t(1) = false;
-%$ catch
-%$    t(1) = true;
-%$ end
-%$
-%$ T = all(t);
+try
+    [a, b] = beta_specification(0.5, .05, 1, 2);
+    t(1) = false;
+catch
+    t(1) = true;
+end
+
+T = all(t);
 %@eof:2
 
 %@test:3
-%$ try
-%$    [a, b] = beta_specification(2.5, .05, 1, 2);
-%$    t(1) = false;
-%$ catch
-%$    t(1) = true;
-%$ end
-%$
-%$ T = all(t);
+try
+    [a, b] = beta_specification(2.5, .05, 1, 2);
+    t(1) = false;
+catch
+    t(1) = true;
+end
+
+T = all(t);
 %@eof:3
 
 %@test:4
-%$ try
-%$    [a, b] = beta_specification(.4, .6*.4+1e-4);
-%$    t(1) = false;
-%$ catch
-%$    t(1) = true;
-%$ end
-%$
-%$ T = all(t);
+try
+    [a, b] = beta_specification(.4, .6*.4+1e-4);
+    t(1) = false;
+catch
+    t(1) = true;
+end
+
+T = all(t);
 %@eof:4
 
 %@test:5
-%$ try
-%$    [a, b] = beta_specification(.4, .6*.4+1e-6);
-%$    t(1) = false;
-%$ catch
-%$    t(1) = true;
-%$ end
-%$ try
-%$    [a, b] = beta_specification(.4, .6*.4-1e-6);
-%$    t(2) = true;
-%$ catch
-%$    t(2) = false;
-%$ end
-%$
-%$ T = all(t);
+try
+   [a, b] = beta_specification(.4, .6*.4+1e-6);
+   t(1) = false;
+catch
+   t(1) = true;
+end
+try
+    [a, b] = beta_specification(.4, .6*.4-1e-6);
+    t(2) = true;
+catch
+    t(2) = false;
+end
+
+T = all(t);
 %@eof:5
\ No newline at end of file
diff --git a/matlab/distributions/compute_prior_mode.m b/matlab/distributions/compute_prior_mode.m
index eed1f802bd3095d084a18c8731f2c08a177b9530..6a9ea7c301d1266b733ccdbf9d4e270a39b08457 100644
--- a/matlab/distributions/compute_prior_mode.m
+++ b/matlab/distributions/compute_prior_mode.m
@@ -1,4 +1,4 @@
-function m = compute_prior_mode(hyperparameters,shape) % --*-- Unitary tests --*--
+function m = compute_prior_mode(hyperparameters,shape)
 % This function computes the mode of the prior distribution given the (two, three or four) hyperparameters
 % of the prior distribution.
 %
@@ -23,7 +23,7 @@ function m = compute_prior_mode(hyperparameters,shape) % --*-- Unitary tests --*
 % [3] The uniform distribution has an infinity of modes. In this case the function returns the prior mean.
 % [4] For the beta distribution we can have 1, 2 or an infinity of modes.
 
-% Copyright © 2009-2017 Dynare Team
+% Copyright © 2009-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -102,140 +102,142 @@ switch shape
     error('Unknown prior shape!')
 end
 
+return % --*-- Unit tests --*--
+
 %@test:1
-%$ % Beta density
-%$ try
-%$     m1 = compute_prior_mode([2 1],1);
-%$     m2 = compute_prior_mode([2 5 1 4],1); % Wolfram Alpha: BetaDistribution[2,5]
-%$     t(1) = true;
-%$ catch
-%$     t(1) = false;
-%$ end
-%$
-%$ % Check the results
-%$ if t(1)
-%$     t(2) = dassert(m1,1,1e-6);
-%$     t(3) = dassert(m2,1/5*3+1,1e-6);
-%$ end
-%$ T = all(t);
+% Beta density
+try
+    m1 = compute_prior_mode([2 1],1);
+    m2 = compute_prior_mode([2 5 1 4],1); % Wolfram Alpha: BetaDistribution[2,5]
+    t(1) = true;
+catch
+    t(1) = false;
+end
+
+% Check the results
+if t(1)
+    t(2) = dassert(m1,1,1e-6);
+    t(3) = dassert(m2,1/5*3+1,1e-6);
+end
+T = all(t);
 %@eof:1
 
 %@test:2
-%$ % Gamma density
-%$ try
-%$     m1 = compute_prior_mode([1 2],2);
-%$     m2 = compute_prior_mode([9 0.5 1],2);  % Wolfram Alpha: GammaDistribution[9,0.5]
-%$     t(1) = true;
-%$ catch
-%$     t(1) = false;
-%$ end
-%$
-%$ % Check the results
-%$ if t(1)
-%$     t(2) = dassert(m1,0,1e-6);
-%$     t(3) = dassert(m2,4+1,1e-6);
-%$ end
-%$ T = all(t);
+% Gamma density
+try
+    m1 = compute_prior_mode([1 2],2);
+    m2 = compute_prior_mode([9 0.5 1],2);  % Wolfram Alpha: GammaDistribution[9,0.5]
+    t(1) = true;
+catch
+    t(1) = false;
+end
+
+% Check the results
+if t(1)
+    t(2) = dassert(m1,0,1e-6);
+    t(3) = dassert(m2,4+1,1e-6);
+end
+T = all(t);
 %@eof:2
 
 %@test:3
-%$ % Normal density
-%$ try
-%$     m1 = compute_prior_mode([1 1],3);
-%$     m2 = compute_prior_mode([2 5],3);
-%$     t(1) = true;
-%$ catch
-%$     t(1) = false;
-%$ end
-%$
-%$ % Check the results
-%$ if t(1)
-%$     t(2) = dassert(m1,1,1e-6);
-%$     t(3) = dassert(m2,2,1e-6);
-%$ end
-%$ T = all(t);
+% Normal density
+try
+    m1 = compute_prior_mode([1 1],3);
+    m2 = compute_prior_mode([2 5],3);
+    t(1) = true;
+catch
+    t(1) = false;
+end
+
+% Check the results
+if t(1)
+    t(2) = dassert(m1,1,1e-6);
+    t(3) = dassert(m2,2,1e-6);
+end
+T = all(t);
 %@eof:3
 
 %@test:4
-%$ % Inverse Gamma I density
-%$ try
-%$     m1 = compute_prior_mode([8 2],4);
-%$     m2 = compute_prior_mode([8 2 1],4);
-%$     t(1) = true;
-%$ catch
-%$     t(1) = false;
-%$ end
-%$
-%$ % Check the results
-%$ if t(1)
-%$     t(2) = dassert(m1,1.632993161855452,1e-6);
-%$     t(3) = dassert(m2,1.632993161855452+1,1e-6);
-%$ end
-%$ T = all(t);
+% Inverse Gamma I density
+try
+    m1 = compute_prior_mode([8 2],4);
+    m2 = compute_prior_mode([8 2 1],4);
+    t(1) = true;
+catch
+    t(1) = false;
+end
+
+% Check the results
+if t(1)
+    t(2) = dassert(m1,1.632993161855452,1e-6);
+    t(3) = dassert(m2,1.632993161855452+1,1e-6);
+end
+T = all(t);
 %@eof:4
 
 %@test:5
-%$ % Uniform density
-%$ try
-%$     m1 = compute_prior_mode([0.5 1/sqrt(12)],5);
-%$     m2 = compute_prior_mode([2 5 1 2],5);
-%$     t(1) = true;
-%$ catch
-%$     t(1) = false;
-%$ end
-%$
-%$ % Check the results
-%$ if t(1)
-%$     t(2) = dassert(m1,0.5,1e-6);
-%$     t(3) = dassert(m2,2,1e-6);
-%$ end
-%$ T = all(t);
+% Uniform density
+try
+    m1 = compute_prior_mode([0.5 1/sqrt(12)],5);
+    m2 = compute_prior_mode([2 5 1 2],5);
+    t(1) = true;
+catch
+    t(1) = false;
+end
+
+% Check the results
+if t(1)
+    t(2) = dassert(m1,0.5,1e-6);
+    t(3) = dassert(m2,2,1e-6);
+end
+T = all(t);
 %@eof:5
 
 %@test:6
-%$ % Inverse Gamma II density, parameterized with s and nu where  s=2*beta and nu=2*alpha
-%$ try
-%$     m1 = compute_prior_mode([8 2],6);  % Wolfram Alpha, parameterized with alpha beta: InversegammaDistribution[1,4]
-%$     m2 = compute_prior_mode([8 4 1],6); % Wolfram Alpha, parameterized with alpha beta: InversegammaDistribution[2,4]
-%$     t(1) = true;
-%$ catch
-%$     t(1) = false;
-%$ end
-%$
-%$ % Check the results
-%$ if t(1)
-%$     t(2) = dassert(m1,2,1e-6);
-%$     t(3) = dassert(m2,1+4/3,1e-6);
-%$ end
-%$ T = all(t);
+% Inverse Gamma II density, parameterized with s and nu where  s=2*beta and nu=2*alpha
+try
+    m1 = compute_prior_mode([8 2],6);  % Wolfram Alpha, parameterized with alpha beta: InversegammaDistribution[1,4]
+    m2 = compute_prior_mode([8 4 1],6); % Wolfram Alpha, parameterized with alpha beta: InversegammaDistribution[2,4]
+    t(1) = true;
+catch
+    t(1) = false;
+end
+
+% Check the results
+if t(1)
+    t(2) = dassert(m1,2,1e-6);
+    t(3) = dassert(m2,1+4/3,1e-6);
+end
+T = all(t);
 %@eof:6
 
 %@test:7
-%$ % Weibull density
-%$ try
-%$     m1 = compute_prior_mode([1 1],8);
-%$     m2 = compute_prior_mode([2 1 1],8); % Wolfram Alpha: WeibullDistribution[2,1]
-%$     t(1) = true;
-%$ catch
-%$     t(1) = false;
-%$ end
-%$
-%$ % Check the results
-%$ if t(1)
-%$     t(2) = dassert(m1,0,1e-6);
-%$     t(3) = dassert(m2,1+1/sqrt(2),1e-6);
-%$ end
-%$ T = all(t);
+% Weibull density
+try
+    m1 = compute_prior_mode([1 1],8);
+    m2 = compute_prior_mode([2 1 1],8); % Wolfram Alpha: WeibullDistribution[2,1]
+    t(1) = true;
+catch
+    t(1) = false;
+end
+
+% Check the results
+if t(1)
+    t(2) = dassert(m1,0,1e-6);
+    t(3) = dassert(m2,1+1/sqrt(2),1e-6);
+end
+T = all(t);
 %@eof:7
 
 %@test:8
-%$ % Unknown density
-%$ try
-%$     m1 = compute_prior_mode([1 1],7);
-%$     t(1) = false;
-%$ catch
-%$     t(1) = true;
-%$ end
-%$
-%$ T = all(t);
-%@eof:8
+% Unknown density
+try
+    m1 = compute_prior_mode([1 1],7);
+    t(1) = false;
+catch
+    t(1) = true;
+end
+
+T = all(t);
+%@eof:8
\ No newline at end of file
diff --git a/matlab/distributions/gamma_specification.m b/matlab/distributions/gamma_specification.m
index dbe1b0657e98f1c0525b8472f9c5cce33fa3dbe8..7019acf50db5fd0319960de544f32507c973b059 100644
--- a/matlab/distributions/gamma_specification.m
+++ b/matlab/distributions/gamma_specification.m
@@ -1,4 +1,4 @@
-function [a, b] = gamma_specification(mu, sigma2, lb, name)   % --*-- Unitary tests --*--
+function [a, b] = gamma_specification(mu, sigma2, lb, name)
 
 % Returns the hyperparameters of the gamma distribution given the expectation and variance.
 %
@@ -12,7 +12,7 @@ function [a, b] = gamma_specification(mu, sigma2, lb, name)   % --*-- Unitary te
 % - a      [double]   First hyperparameter of the Gamma density (shape).
 % - b      [double]   Second hyperparameter of the Gamma density (scale).
 
-% Copyright © 2015-2017 Dynare Team
+% Copyright © 2015-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -53,54 +53,56 @@ mu = mu-lb;
 b  = sigma2/mu;
 a  = mu/b;
 
+return % --*-- Unit tests --*--
+
 %@test:1
-%$ try
-%$    [a, b] = gamma_specification(.5, 1);
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = abs(0.5-a*b)<1e-12;
-%$    t(3) = abs(1.0-a*b*b)<1e-12;
-%$ end
-%$ T = all(t);
+try
+   [a, b] = gamma_specification(.5, 1);
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = abs(0.5-a*b)<1e-12;
+    t(3) = abs(1.0-a*b*b)<1e-12;
+end
+T = all(t);
 %@eof:1
 
 %@test:2
-%$ try
-%$    [a, b] = gamma_specification(2, 1);
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = abs(2.0-a*b)<1e-12;
-%$    t(3) = abs(1.0-a*b*b)<1e-12;
-%$ end
-%$ T = all(t);
+try
+   [a, b] = gamma_specification(2, 1);
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = abs(2.0-a*b)<1e-12;
+    t(3) = abs(1.0-a*b*b)<1e-12;
+end
+T = all(t);
 %@eof:2
 
 %@test:3
-%$ try
-%$    [a, b] = gamma_specification(2, 1, 3);
-%$    t(1) = false;
-%$ catch
-%$    t(1) = true;
-%$ end
-%$
-%$ T = all(t);
+try
+    [a, b] = gamma_specification(2, 1, 3);
+    t(1) = false;
+catch
+    t(1) = true;
+end
+
+T = all(t);
 %@eof:3
 
 %@test:4
-%$ try
-%$    [a, b] = gamma_specification(2, Inf, 3);
-%$    t(1) = false;
-%$ catch
-%$    t(1) = true;
-%$ end
-%$
-%$ T = all(t);
+try
+    [a, b] = gamma_specification(2, Inf, 3);
+    t(1) = false;
+catch
+    t(1) = true;
+end
+
+T = all(t);
 %@eof:4
\ No newline at end of file
diff --git a/matlab/distributions/inverse_gamma_specification.m b/matlab/distributions/inverse_gamma_specification.m
index cedca4d6700b036dbdf1deaf38aac0de893533b6..e5136381ad7591d55a029d6a25b924eb6b08a74f 100644
--- a/matlab/distributions/inverse_gamma_specification.m
+++ b/matlab/distributions/inverse_gamma_specification.m
@@ -1,4 +1,4 @@
-function [s,nu] = inverse_gamma_specification(mu, sigma2, lb, type, use_fzero_flag, name) % --*-- Unitary tests --*--
+function [s,nu] = inverse_gamma_specification(mu, sigma2, lb, type, use_fzero_flag, name)
 
 % Computes the inverse Gamma hyperparameters from the prior mean and variance.
 %
@@ -21,7 +21,7 @@ function [s,nu] = inverse_gamma_specification(mu, sigma2, lb, type, use_fzero_fl
 % more often in finding an interval for nu containing a signe change because it expands the interval on both sides and eventually
 % violates  the condition nu>2.
 
-% Copyright © 2003-2017 Dynare Team
+% Copyright © 2003-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -147,62 +147,64 @@ else
     error('inverse_gamma_specification: unkown type')
 end
 
+return % --*-- Unit tests --*--
+
 %@test:1
-%$ try
-%$    [s, nu] = inverse_gamma_specification(.5, .05, 0, 1);
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = abs(0.5-sqrt(.5*s)*gamma(.5*(nu-1))/gamma(.5*nu))<1e-12;
-%$    t(3) = abs(0.05-s/(nu-2)+.5^2)<1e-12;
-%$ end
-%$ T = all(t);
+try
+   [s, nu] = inverse_gamma_specification(.5, .05, 0, 1);
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = abs(0.5-sqrt(.5*s)*gamma(.5*(nu-1))/gamma(.5*nu))<1e-12;
+    t(3) = abs(0.05-s/(nu-2)+.5^2)<1e-12;
+end
+T = all(t);
 %@eof:1
 
 %@test:2
-%$ try
-%$    [s, nu] = inverse_gamma_specification(.5, .05, 0, 2);
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = abs(0.5-s/(nu-2))<1e-12;
-%$    t(3) = abs(0.05-2*.5^2/(nu-4))<1e-12;
-%$ end
-%$ T = all(t);
+try
+   [s, nu] = inverse_gamma_specification(.5, .05, 0, 2);
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = abs(0.5-s/(nu-2))<1e-12;
+    t(3) = abs(0.05-2*.5^2/(nu-4))<1e-12;
+end
+T = all(t);
 %@eof:2
 
 %@test:3
-%$ try
-%$    [s, nu] = inverse_gamma_specification(.5, Inf, 0, 1);
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = abs(0.5-sqrt(.5*s)*gamma(.5*(nu-1))/gamma(.5*nu))<1e-12;
-%$    t(3) = isequal(nu, 2); %abs(0.05-2*.5^2/(nu-4))<1e-12;
-%$ end
-%$ T = all(t);
+try
+   [s, nu] = inverse_gamma_specification(.5, Inf, 0, 1);
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = abs(0.5-sqrt(.5*s)*gamma(.5*(nu-1))/gamma(.5*nu))<1e-12;
+    t(3) = isequal(nu, 2); %abs(0.05-2*.5^2/(nu-4))<1e-12;
+end
+T = all(t);
 %@eof:3
 
 %@test:4
-%$ try
-%$    [s, nu] = inverse_gamma_specification(.5, Inf, 0, 2);
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = abs(0.5-s/(nu-2))<1e-12;
-%$    t(3) = isequal(nu, 4);
-%$ end
-%$ T = all(t);
+try
+   [s, nu] = inverse_gamma_specification(.5, Inf, 0, 2);
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = abs(0.5-s/(nu-2))<1e-12;
+    t(3) = isequal(nu, 4);
+end
+T = all(t);
 %@eof:4
\ No newline at end of file
diff --git a/matlab/distributions/lpdfgweibull.m b/matlab/distributions/lpdfgweibull.m
index 29e0212d165976b0be265ca098a9525f57e94f2f..b6c88da94b3bf92199ccdcd8b183063a73810ca6 100644
--- a/matlab/distributions/lpdfgweibull.m
+++ b/matlab/distributions/lpdfgweibull.m
@@ -1,4 +1,4 @@
-function [ldens,Dldens,D2ldens] = lpdfgweibull(x,a,b,c)  % --*-- Unitary tests --*--
+function [ldens,Dldens,D2ldens] = lpdfgweibull(x,a,b,c)
 
 % Evaluates the logged Weibull PDF at x.
 %
@@ -16,7 +16,7 @@ function [ldens,Dldens,D2ldens] = lpdfgweibull(x,a,b,c)  % --*-- Unitary tests -
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2015-2020 Dynare Team
+% Copyright © 2015-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -96,344 +96,346 @@ if nargout>1
     end
 end
 
+return % --*-- Unit tests --*--
+
 %@test:1
-%$ try
-%$    lpdfgweibull(1.0);
-%$    t(1) = false;
-%$ catch
-%$    t(1) = true;
-%$ end
-%$
-%$ T = all(t);
+try
+    lpdfgweibull(1.0);
+    t(1) = false;
+catch
+    t(1) = true;
+end
+
+T = all(t);
 %@eof:1
 
 %@test:2
-%$ try
-%$    lpdfgweibull(1.0, .5);
-%$    t(1) = false;
-%$ catch
-%$    t(1) = true;
-%$ end
-%$
-%$ T = all(t);
+try
+    lpdfgweibull(1.0, .5);
+    t(1) = false;
+catch
+    t(1) = true;
+end
+
+T = all(t);
 %@eof:2
 
 %@test:3
-%$ try
-%$    lpdfgweibull(ones(2,2), .5*ones(2,1), .1*ones(3,1));
-%$    t(1) = false;
-%$ catch
-%$    t(1) = true;
-%$ end
-%$
-%$ T = all(t);
+try
+    lpdfgweibull(ones(2,2), .5*ones(2,1), .1*ones(3,1));
+    t(1) = false;
+catch
+    t(1) = true;
+end
+
+T = all(t);
 %@eof:3
 
 %@test:4
-%$ try
-%$    a = lpdfgweibull(-1, .5, .1);
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = isinf(a);
-%$ end
-%$
-%$ T = all(t);
+try
+   a = lpdfgweibull(-1, .5, .1);
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = isinf(a);
+end
+
+T = all(t);
 %@eof:4
 
 %@test:5
-%$ try
-%$    a = lpdfgweibull(0, 1, 1);
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = abs(a-1.0)<1e-10;
-%$ end
-%$
-%$ T = all(t);
+try
+   a = lpdfgweibull(0, 1, 1);
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = abs(a-1.0)<1e-10;
+end
+
+T = all(t);
 %@eof:5
 
 %@test:6
-%$ try
-%$    a = lpdfgweibull([0, -1], [1 1], [1 1]);
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = abs(a(1)-1.0)<1e-10;
-%$    t(3) = isinf(a(2));
-%$ end
-%$
-%$ T = all(t);
+try
+   a = lpdfgweibull([0, -1], [1 1], [1 1]);
+   t(1) = true;
+catch
+   t(1) = false;
+end
+
+if t(1)
+    t(2) = abs(a(1)-1.0)<1e-10;
+    t(3) = isinf(a(2));
+end
+
+T = all(t);
 %@eof:6
 
 %@test:7
-%$ scale = 1;
-%$ shape = 2;
-%$ mode  = scale*((shape-1)/shape)^(1/shape);
-%$
-%$ try
-%$    [a, b, c] = lpdfgweibull(mode, shape, scale);
-%$    p = rand(1000,1)*4;
-%$    am = lpdfgweibull(p, shape*ones(size(p)), scale*ones(size(p)));
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$
-%$    t(2) = abs(b)<1e-8;
-%$    t(3) = c<0;
-%$    t(4) = all(am<a);
-%$ end
-%$
-%$ T = all(t);
+scale = 1;
+shape = 2;
+mode  = scale*((shape-1)/shape)^(1/shape);
+
+try
+   [a, b, c] = lpdfgweibull(mode, shape, scale);
+   p = rand(1000,1)*4;
+   am = lpdfgweibull(p, shape*ones(size(p)), scale*ones(size(p)));
+   t(1) = true;
+catch
+   t(1) = false;
+end
+
+if t(1)
+
+    t(2) = abs(b)<1e-8;
+    t(3) = c<0;
+    t(4) = all(am<a);
+end
+
+T = all(t);
 %@eof:7
 
 %@test:8
-%$ scale = 1;
-%$ shape = 2;
-%$ density  = @(x) exp(lpdfgweibull(x,shape,scale));
-%$
-%$ try
-%$    if isoctave
-%$        s = quadl(density, .0000000001, 100000, 1e-10);
-%$    else
-%$        s = integral(density, 0, 100000);
-%$    end
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = abs(s-1)<1e-6;
-%$ end
-%$
-%$ T = all(t);
+scale = 1;
+shape = 2;
+density  = @(x) exp(lpdfgweibull(x,shape,scale));
+
+try
+   if isoctave
+       s = quadl(density, .0000000001, 100000, 1e-10);
+   else
+       s = integral(density, 0, 100000);
+   end
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = abs(s-1)<1e-6;
+end
+
+T = all(t);
 %@eof:8
 
 %@test:9
-%$ scale = 1;
-%$ shape = 1;
-%$ density  = @(x) exp(lpdfgweibull(x,shape,scale));
-%$
-%$ try
-%$    if isoctave
-%$        s = quadl(density, .0000000001, 100000, 1e-10);
-%$    else
-%$        s = integral(density, 0, 100000);
-%$    end
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = abs(s-1)<1e-6;
-%$ end
-%$
-%$ T = all(t);
+scale = 1;
+shape = 1;
+density  = @(x) exp(lpdfgweibull(x,shape,scale));
+
+try
+   if isoctave
+       s = quadl(density, .0000000001, 100000, 1e-10);
+   else
+       s = integral(density, 0, 100000);
+   end
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = abs(s-1)<1e-6;
+end
+
+T = all(t);
 %@eof:9
 
 %@test:10
-%$ scale = 1;
-%$ shape = .5;
-%$ density  = @(x) exp(lpdfgweibull(x,shape,scale));
-%$
-%$ try
-%$    if isoctave
-%$        s = quadl(density, .0000000001, 100000, 1e-10);
-%$    else
-%$        s = integral(density, 0, 100000);
-%$    end
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$        if isoctave
-%$            t(2) = abs(s-1)<5e-5;
-%$        else
-%$            t(2) = abs(s-1)<1e-6;
-%$        end
-%$ end
-%$
-%$ T = all(t);
+scale = 1;
+shape = .5;
+density  = @(x) exp(lpdfgweibull(x,shape,scale));
+
+try
+   if isoctave
+       s = quadl(density, .0000000001, 100000, 1e-10);
+   else
+       s = integral(density, 0, 100000);
+   end
+   t(1) = true;
+catch
+   t(1) = false;
+end
+
+if t(1)
+    if isoctave
+        t(2) = abs(s-1)<5e-5;
+    else
+        t(2) = abs(s-1)<1e-6;
+    end
+end
+
+T = all(t);
 %@eof:10
 
 %@test:11
-%$ scale = 1;
-%$ shape = 2;
-%$ xdens = @(x) x.*exp(lpdfgweibull(x,shape,scale));
-%$
-%$ try
-%$    if isoctave
-%$        s = quadgk(xdens, .0000000001, 100000, 1e-10);
-%$    else
-%$        s = integral(xdens, 0, 100000);
-%$    end
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = abs(s-scale*gamma(1+1/shape))<1e-6;
-%$ end
-%$
-%$ T = all(t);
+scale = 1;
+shape = 2;
+xdens = @(x) x.*exp(lpdfgweibull(x,shape,scale));
+
+try
+   if isoctave
+       s = quadgk(xdens, .0000000001, 100000, 1e-10);
+   else
+       s = integral(xdens, 0, 100000);
+   end
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = abs(s-scale*gamma(1+1/shape))<1e-6;
+end
+
+T = all(t);
 %@eof:11
 
 %@test:12
-%$ scale = 1;
-%$ shape = 1;
-%$ xdens = @(x) x.*exp(lpdfgweibull(x,shape,scale));
-%$
-%$ try
-%$    if isoctave
-%$        s = quadl(xdens, .0000000001, 100000, 1e-10);
-%$    else
-%$        s = integral(xdens, 0, 100000);
-%$    end
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = abs(s-scale*gamma(1+1/shape))<1e-6;
-%$ end
-%$
-%$ T = all(t);
+scale = 1;
+shape = 1;
+xdens = @(x) x.*exp(lpdfgweibull(x,shape,scale));
+
+try
+   if isoctave
+       s = quadl(xdens, .0000000001, 100000, 1e-10);
+   else
+       s = integral(xdens, 0, 100000);
+   end
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = abs(s-scale*gamma(1+1/shape))<1e-6;
+end
+
+T = all(t);
 %@eof:12
 
 %@test:13
-%$ scale = 1;
-%$ shape = .5;
-%$ xdens = @(x) x.*exp(lpdfgweibull(x,shape,scale));
-%$
-%$ try
-%$    if isoctave
-%$        s = quadl(xdens, .0000000001, 100000, 1e-10);
-%$    else
-%$        s = integral(xdens, 0, 100000);
-%$    end
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = abs(s-scale*gamma(1+1/shape))<1e-6;
-%$ end
-%$
-%$ T = all(t);
+scale = 1;
+shape = .5;
+xdens = @(x) x.*exp(lpdfgweibull(x,shape,scale));
+
+try
+   if isoctave
+       s = quadl(xdens, .0000000001, 100000, 1e-10);
+   else
+       s = integral(xdens, 0, 100000);
+   end
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = abs(s-scale*gamma(1+1/shape))<1e-6;
+end
+
+T = all(t);
 %@eof:13
 
 %@test:14
-%$ scale = 1;
-%$ shape = 2;
-%$ density = @(x) exp(lpdfgweibull(x,shape,scale));
-%$ n = 200;
-%$
-%$ try
-%$    s = NaN(n, 1);
-%$    for i=1:n
-%$        if isoctave
-%$             s(i) = quadl(density, .0000000001, .1*i, 1e-10);
-%$        else
-%$            s(i) = integral(density, 0, .1*i);
-%$        end
-%$    end
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    for i=1:n
-%$        x = .1*i;
-%$        q = 1-exp(-(x/scale)^shape);
-%$        t(i+1) = abs(s(i)-q)<1e-6;
-%$    end
-%$ end
-%$
-%$ T = all(t);
+scale = 1;
+shape = 2;
+density = @(x) exp(lpdfgweibull(x,shape,scale));
+n = 200;
+
+try
+   s = NaN(n, 1);
+   for i=1:n
+       if isoctave
+            s(i) = quadl(density, .0000000001, .1*i, 1e-10);
+       else
+           s(i) = integral(density, 0, .1*i);
+       end
+   end
+   t(1) = true;
+catch
+   t(1) = false;
+end
+
+if t(1)
+    for i=1:n
+        x = .1*i;
+        q = 1-exp(-(x/scale)^shape);
+        t(i+1) = abs(s(i)-q)<1e-6;
+    end
+end
+
+T = all(t);
 %@eof:14
 
 %@test:15
-%$ scale = 1;
-%$ shape = 1;
-%$ density = @(x) exp(lpdfgweibull(x,shape,scale));
-%$ n = 200;
-%$
-%$ try
-%$    s = NaN(n, 1);
-%$    for i=1:n
-%$        if isoctave
-%$            s(i) = quadl(density, .0000000001, .1*i, 1e-10);
-%$        else
-%$            s(i) = integral(density, 0, .1*i);
-%$        end
-%$    end
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    for i=1:n
-%$        x = .1*i;
-%$        q = 1-exp(-(x/scale)^shape);
-%$        t(i+1) = abs(s(i)-q)<1e-6;
-%$    end
-%$ end
-%$
-%$ T = all(t);
+scale = 1;
+shape = 1;
+density = @(x) exp(lpdfgweibull(x,shape,scale));
+n = 200;
+
+try
+   s = NaN(n, 1);
+   for i=1:n
+       if isoctave
+           s(i) = quadl(density, .0000000001, .1*i, 1e-10);
+       else
+           s(i) = integral(density, 0, .1*i);
+       end
+   end
+   t(1) = true;
+catch
+   t(1) = false;
+end
+
+if t(1)
+    for i=1:n
+        x = .1*i;
+        q = 1-exp(-(x/scale)^shape);
+        t(i+1) = abs(s(i)-q)<1e-6;
+    end
+end
+
+T = all(t);
 %@eof:15
 
 %@test:16
-%$ scale = 1;
-%$ shape = .5;
-%$ density = @(x) exp(lpdfgweibull(x,shape,scale));
-%$ n = 200;
-%$
-%$ try
-%$    s = NaN(n, 1);
-%$    for i=1:n
-%$        if isoctave
-%$            s(i) = quadl(density, .0000000001, .1*i, 1e-10);
-%$        else
-%$            s(i) = integral(density, 0, .1*i);
-%$        end
-%$    end
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    for i=1:n
-%$        x = .1*i;
-%$        q = 1-exp(-(x/scale)^shape);
-%$        if isoctave
-%$            t(i+1) = abs(s(i)-q)<5e-5;
-%$        else
-%$            t(i+1) = abs(s(i)-q)<1e-6;
-%$        end
-%$    end
-%$ end
-%$
-%$ T = all(t);
+scale = 1;
+shape = .5;
+density = @(x) exp(lpdfgweibull(x,shape,scale));
+n = 200;
+
+try
+   s = NaN(n, 1);
+   for i=1:n
+       if isoctave
+           s(i) = quadl(density, .0000000001, .1*i, 1e-10);
+       else
+           s(i) = integral(density, 0, .1*i);
+       end
+   end
+   t(1) = true;
+catch
+   t(1) = false;
+end
+
+if t(1)
+   for i=1:n
+       x = .1*i;
+       q = 1-exp(-(x/scale)^shape);
+       if isoctave
+           t(i+1) = abs(s(i)-q)<5e-5;
+       else
+           t(i+1) = abs(s(i)-q)<1e-6;
+       end
+   end
+end
+
+T = all(t);
 %@eof:16
diff --git a/matlab/distributions/weibull_specification.m b/matlab/distributions/weibull_specification.m
index 56becdef17af9b4f34d37be45eb3664540ca0d97..e9808f1d0ec1c31419726ecc4a965720f2997abd 100644
--- a/matlab/distributions/weibull_specification.m
+++ b/matlab/distributions/weibull_specification.m
@@ -1,4 +1,4 @@
-function [shape, scale] = weibull_specification(mu, sigma2, lb, name)   % --*-- Unitary tests --*--
+function [shape, scale] = weibull_specification(mu, sigma2, lb, name)
 
 % Returns the hyperparameters of the Weibull distribution given the expectation and variance.
 %
@@ -9,7 +9,7 @@ function [shape, scale] = weibull_specification(mu, sigma2, lb, name)   % --*--
 %
 %
 
-% Copyright © 2015-2017 Dynare Team
+% Copyright © 2015-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -65,78 +65,80 @@ end
 
 scale = mu/gamma(1+1/shape);
 
+return % --*-- Unit tests --*--
+
 %@test:1
-%$ debug = false;
-%$ scales = 1:.01:5;
-%$ shapes = .5:.01:2;
-%$ n_scales = length(scales);
-%$ n_shapes = length(shapes);
-%$ mu = NaN(n_scales, n_shapes);
-%$ s2 = NaN(n_scales, n_shapes);
-%$ for i=1:n_shapes
-%$     g1 = gamma(1+1/shapes(i));
-%$     g2 = gamma(1+2/shapes(i));
-%$     g3 = g1*g1;
-%$     for j=1:n_scales
-%$         mu(j, i) = scales(j)*g1;
-%$         s2(j, i) = scales(j)*scales(j)*(g2-g3);
-%$     end
-%$ end
-%$ if debug
-%$    success = [];
-%$    failed1 = [];
-%$    failed1_ = [];
-%$    failed2 = [];
-%$ end
-%$ try
-%$    for i=1:n_shapes
-%$        for j=1:n_scales
-%$           if debug
-%$               disp(sprintf('... mu=%s and s2=%s', num2str(mu(j,i)),num2str(s2(j,i))))
-%$           end
-%$           if ~isnan(mu(j,i)) && ~isnan(s2(j,i)) && ~isinf(mu(j,i)) && ~isinf(s2(j,i))
-%$               [shape, scale] = weibull_specification(mu(j,i), s2(j,i));
-%$               if isnan(scale)
-%$                  t = false;
-%$               else
-%$                  if abs(scales(j)-scale)<1e-9 && abs(shapes(i)-shape)<1e-9
-%$                      t = true;
-%$                  else
-%$                      t = false;
-%$                  end
-%$               end
-%$               if ~t && debug
-%$                  failed1 = [failed1; mu(j,i) s2(j,i)];
-%$                  failed1_ = [failed1_; shapes(i) scales(j)];
-%$                  error('UnitTest','Cannot compute scale and shape hyperparameters for mu=%s and s2=%s', num2str(mu(j,i)), num2str(s2(j,i)))
-%$               end
-%$               if debug
-%$                   success = [success; mu(j,i) s2(j,i)];
-%$               end
-%$           else
-%$               failed2 = [failed2; shapes(i) scales(j)];
-%$               continue % Pass this test
-%$           end
-%$       end
-%$   end
-%$ catch
-%$    t = false;
-%$ end
-%$
-%$ if debug
-%$     figure(1)
-%$     plot(success(:,1),success(:,2),'ok');
-%$     if ~isempty(failed1)
-%$         hold on
-%$         plot(failed1(:,1),failed1(:,2),'or');
-%$         hold off
-%$         figure(2)
-%$         plot(failed1_(:,1),failed1_(:,2),'or')
-%$     end
-%$     if ~isempty(failed2)
-%$         figure(2)
-%$         plot(failed2(:,1),failed2(:,2),'or');
-%$     end
-%$ end
-%$ T = all(t);
+debug = false;
+scales = 1:.01:5;
+shapes = .5:.01:2;
+n_scales = length(scales);
+n_shapes = length(shapes);
+mu = NaN(n_scales, n_shapes);
+s2 = NaN(n_scales, n_shapes);
+for i=1:n_shapes
+    g1 = gamma(1+1/shapes(i));
+    g2 = gamma(1+2/shapes(i));
+    g3 = g1*g1;
+    for j=1:n_scales
+        mu(j, i) = scales(j)*g1;
+        s2(j, i) = scales(j)*scales(j)*(g2-g3);
+    end
+end
+if debug
+   success = [];
+   failed1 = [];
+   failed1_ = [];
+   failed2 = [];
+end
+try
+   for i=1:n_shapes
+       for j=1:n_scales
+          if debug
+              disp(sprintf('... mu=%s and s2=%s', num2str(mu(j,i)),num2str(s2(j,i))))
+          end
+          if ~isnan(mu(j,i)) && ~isnan(s2(j,i)) && ~isinf(mu(j,i)) && ~isinf(s2(j,i))
+              [shape, scale] = weibull_specification(mu(j,i), s2(j,i));
+              if isnan(scale)
+                 t = false;
+              else
+                 if abs(scales(j)-scale)<1e-9 && abs(shapes(i)-shape)<1e-9
+                     t = true;
+                 else
+                     t = false;
+                 end
+              end
+              if ~t && debug
+                 failed1 = [failed1; mu(j,i) s2(j,i)];
+                 failed1_ = [failed1_; shapes(i) scales(j)];
+                 error('UnitTest','Cannot compute scale and shape hyperparameters for mu=%s and s2=%s', num2str(mu(j,i)), num2str(s2(j,i)))
+              end
+              if debug
+                  success = [success; mu(j,i) s2(j,i)];
+              end
+          else
+              failed2 = [failed2; shapes(i) scales(j)];
+              continue % Pass this test
+          end
+      end
+  end
+catch
+   t = false;
+end
+
+if debug
+    figure(1)
+    plot(success(:,1),success(:,2),'ok');
+    if ~isempty(failed1)
+        hold on
+        plot(failed1(:,1),failed1(:,2),'or');
+        hold off
+        figure(2)
+        plot(failed1_(:,1),failed1_(:,2),'or')
+    end
+    if ~isempty(failed2)
+        figure(2)
+        plot(failed2(:,1),failed2(:,2),'or');
+    end
+end
+T = all(t);
 %@eof:1
diff --git a/matlab/hessian.m b/matlab/hessian.m
index a0aafc9b1c19e191cdd8a153e0b9f5b3a026297d..2df67f431665eb15c8c720aa7a4abfb9826c1992 100644
--- a/matlab/hessian.m
+++ b/matlab/hessian.m
@@ -1,4 +1,4 @@
-function hessian_mat = hessian(func,x, gstep, varargin) % --*-- Unitary tests --*--
+function hessian_mat = hessian(func,x, gstep, varargin)
 
 % Computes second order partial derivatives
 %
@@ -26,7 +26,7 @@ function hessian_mat = hessian(func,x, gstep, varargin) % --*-- Unitary tests --
 %    none
 %
 
-% Copyright © 2001-2017 Dynare Team
+% Copyright © 2001-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -98,53 +98,54 @@ for i=1:n
     end
 end
 
+return % --*-- Unit tests --*--
 
 %@test:1
-%$ % Create a function.
-%$ fid = fopen('exfun.m','w+');
-%$ fprintf(fid,'function [f,g,H] = exfun(xvar)\\n');
-%$ fprintf(fid,'x = xvar(1);\\n');
-%$ fprintf(fid,'y = xvar(2);\\n');
-%$ fprintf(fid,'f = x^2* log(y);\\n');
-%$ fprintf(fid,'if nargout>1\\n');
-%$ fprintf(fid,'    g = zeros(2,1);\\n');
-%$ fprintf(fid,'    g(1) = 2*x*log(y);\\n');
-%$ fprintf(fid,'    g(2) = x*x/y;\\n');
-%$ fprintf(fid,'end\\n');
-%$ fprintf(fid,'if nargout>2\\n');
-%$ fprintf(fid,'    H = zeros(2,2);\\n');
-%$ fprintf(fid,'    H(1,1) = 2*log(y);\\n');
-%$ fprintf(fid,'    H(1,2) = 2*x/y;\\n');
-%$ fprintf(fid,'    H(2,1) = H(1,2);\\n');
-%$ fprintf(fid,'    H(2,2) = -x*x/(y*y);\\n');
-%$ fprintf(fid,'    H = H(:);\\n');
-%$ fprintf(fid,'end\\n');
-%$ fclose(fid);
-%$
-%$ rehash;
-%$
-%$ t = zeros(5,1);
-%$
-%$ % Evaluate the Hessian at (1,e)
-%$ try
-%$    H = hessian('exfun',[1; exp(1)],[1e-2; 1]);
-%$    t(1) = 1;
-%$ catch
-%$    t(1) = 0;
-%$ end
-%$
-%$ % Compute the true Hessian matrix
-%$ [f, g, Htrue] = exfun([1 exp(1)]);
-%$
-%$ % Delete exfun routine from disk.
-%$ delete('exfun.m');
-%$
-%$ % Compare the values in H and Htrue
-%$ if t(1)
-%$    t(2) = dassert(abs(H(1)-Htrue(1))<1e-6,true);
-%$    t(3) = dassert(abs(H(2)-Htrue(2))<1e-6,true);
-%$    t(4) = dassert(abs(H(3)-Htrue(3))<1e-6,true);
-%$    t(5) = dassert(abs(H(4)-Htrue(4))<1e-6,true);
-%$ end
-%$ T = all(t);
+% Create a function.
+fid = fopen('exfun.m','w+');
+fprintf(fid,'function [f,g,H] = exfun(xvar)\\n');
+fprintf(fid,'x = xvar(1);\\n');
+fprintf(fid,'y = xvar(2);\\n');
+fprintf(fid,'f = x^2* log(y);\\n');
+fprintf(fid,'if nargout>1\\n');
+fprintf(fid,'    g = zeros(2,1);\\n');
+fprintf(fid,'    g(1) = 2*x*log(y);\\n');
+fprintf(fid,'    g(2) = x*x/y;\\n');
+fprintf(fid,'end\\n');
+fprintf(fid,'if nargout>2\\n');
+fprintf(fid,'    H = zeros(2,2);\\n');
+fprintf(fid,'    H(1,1) = 2*log(y);\\n');
+fprintf(fid,'    H(1,2) = 2*x/y;\\n');
+fprintf(fid,'    H(2,1) = H(1,2);\\n');
+fprintf(fid,'    H(2,2) = -x*x/(y*y);\\n');
+fprintf(fid,'    H = H(:);\\n');
+fprintf(fid,'end\\n');
+fclose(fid);
+
+rehash;
+
+t = zeros(5,1);
+
+% Evaluate the Hessian at (1,e)
+try
+   H = hessian('exfun',[1; exp(1)],[1e-2; 1]);
+   t(1) = 1;
+catch
+   t(1) = 0;
+end
+
+% Compute the true Hessian matrix
+[f, g, Htrue] = exfun([1 exp(1)]);
+
+% Delete exfun routine from disk.
+delete('exfun.m');
+
+% Compare the values in H and Htrue
+if t(1)
+    t(2) = dassert(abs(H(1)-Htrue(1))<1e-6,true);
+    t(3) = dassert(abs(H(2)-Htrue(2))<1e-6,true);
+    t(4) = dassert(abs(H(3)-Htrue(3))<1e-6,true);
+    t(5) = dassert(abs(H(4)-Htrue(4))<1e-6,true);
+end
+T = all(t);
 %@eof:1
diff --git a/matlab/internals.m b/matlab/internals.m
index 3b09b842814d8d5895d37e02cb36fae1e484d57b..670e5b8c8c9c882fd47c6242a8073ade749936ed 100644
--- a/matlab/internals.m
+++ b/matlab/internals.m
@@ -4,13 +4,13 @@ function internals(flag, varargin)
 %! @deftypefn {Function File} internals (@var{flag},@var{a},@var{b}, ...)
 %! @anchor{internals}
 %! @sp 1
-%! This command provides internal documentation and unitary tests for the matlab routines.
+%! This command provides internal documentation and unit tests for the matlab routines.
 %! @sp 2
 %! @strong{Inputs}
 %! @sp 1
 %! @table @ @var
 %! @item flag
-%! Mandatory argument: --doc (for displaying internal documentation) or --test (for performing unitary tests).
+%! Mandatory argument: --doc (for displaying internal documentation) or --test (for performing unit tests).
 %! @item b
 %! Name of the routine to be tested or for which internal documentation is needed.
 %! @item c
@@ -36,13 +36,13 @@ function internals(flag, varargin)
 %! @example
 %! internals --test particle/local_state_iteration
 %! @end example
-%! will execute the unitary tests associated the routine local_state_iteration.
+%! will execute the unit tests associated the routine local_state_iteration.
 %! @sp 2
 %! @strong{Remarks}
 %! @sp 1
 %! [1] It is not possible to display the internal documentation of more than one routine.
 %! @sp 1
-%! [2] It is possible to perform unitary tests on a list of routines.
+%! [2] It is possible to perform unit tests on a list of routines.
 %! @sp 1
 %! [3] For displaying the internal documentation, matlab calls texinfo which has to be installed.
 %! @sp 2
@@ -56,7 +56,7 @@ function internals(flag, varargin)
 %! @end deftypefn
 %@eod:
 
-% Copyright © 2011-2014 Dynare Team
+% Copyright © 2011-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
diff --git a/matlab/isolder.m b/matlab/isolder.m
index c3f8e1d3a3d769e0019a8b67aedb7d80164d9eca..237395436ff5ec52e476d05d87fdeb2f37328f50 100644
--- a/matlab/isolder.m
+++ b/matlab/isolder.m
@@ -1,4 +1,4 @@
-function b = isolder(f, F) % --*-- Unitary tests --*--
+function b = isolder(f, F)
 
 % Returns true if f is older than any file in folder (and subfolders) F.
 %
@@ -9,7 +9,7 @@ function b = isolder(f, F) % --*-- Unitary tests --*--
 % OUTPUT
 % - b   [logical]
 
-% Copyright © 2015-2017 Dynare Team
+% Copyright © 2015-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -56,38 +56,40 @@ for i=1:length(files)
     end
 end
 
+return % --*-- Unit tests --*--
+
 %@test:1
-%$ t = false(3,1);
-%$ mkdir toto;
-%$ cd toto
-%$ mkdir titi;
-%$ mkdir tata;
-%$ cd tata
-%$ mkdir tutu;
-%$ cd tutu
-%$ system('touch a.m');
-%$ cd ..
-%$ system('touch b.m');
-%$ system('touch c.m');
-%$ cd ../titi
-%$ system('touch d.m');
-%$ cd ..
-%$ pause(1)
-%$ system('touch e.m');
-%$ t(1) = isequal(isolder('e.m'), false);
-%$ pause(1)
-%$ system('touch tata/tutu/a.m');
-%$ system('touch tata/b.m');
-%$ system('touch tata/c.m');
-%$ system('touch titi/d.m');
-%$ t(2) = isequal(isolder('e.m'), true);
-%$ pause(1)
-%$ system('touch e.m');
-%$ t(3) = isequal(isolder('e.m'), false);
-%$ cd ..
-%$ if isoctave()
-%$   confirm_recursive_rmdir(false, 'local');
-%$ end
-%$ rmdir('toto','s');
-%$ T = all(t);
+t = false(3,1);
+mkdir toto;
+cd toto
+mkdir titi;
+mkdir tata;
+cd tata
+mkdir tutu;
+cd tutu
+system('touch a.m');
+cd ..
+system('touch b.m');
+system('touch c.m');
+cd ../titi
+system('touch d.m');
+cd ..
+pause(1)
+system('touch e.m');
+t(1) = isequal(isolder('e.m'), false);
+pause(1)
+system('touch tata/tutu/a.m');
+system('touch tata/b.m');
+system('touch tata/c.m');
+system('touch titi/d.m');
+t(2) = isequal(isolder('e.m'), true);
+pause(1)
+system('touch e.m');
+t(3) = isequal(isolder('e.m'), false);
+cd ..
+if isoctave()
+    confirm_recursive_rmdir(false, 'local');
+end
+rmdir('toto','s');
+T = all(t);
 %@eof:1
diff --git a/matlab/load_m_file_data_legacy.m b/matlab/load_m_file_data_legacy.m
index f3782cbb5e3cf3364c6328deb01f4eab3eba3ee9..8b2afe4faa2a552416c1e4d8729c9b5ba49a0acd 100644
--- a/matlab/load_m_file_data_legacy.m
+++ b/matlab/load_m_file_data_legacy.m
@@ -1,6 +1,6 @@
-function o2WysrOISH  = load_m_file_data_legacy(datafile, U7ORsJ0vy3) % --*-- Unitary tests --*--
+function o2WysrOISH  = load_m_file_data_legacy(datafile, U7ORsJ0vy3)
 
-% Copyright © 2014-2017 Dynare Team
+% Copyright © 2014-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -78,7 +78,7 @@ JSmvfqTSXT = repmat(' vec(%s) ', 1, length(U7ORsJ0vy3));
 VbO4y7zOlh = sprintf('[%s]', JSmvfqTSXT);
 o2WysrOISH = dseries(eval(sprintf(VbO4y7zOlh, U7ORsJ0vy3{:})), [], U7ORsJ0vy3);
 
-return
+return % --*-- Unit tests --*--
 
 %@test:1
 % Write a data file
diff --git a/matlab/lyapunov_solver.m b/matlab/lyapunov_solver.m
index c467e6aabfbfaa1f7ef00b1786ba294cff74368e..8dcb3b9e66de1e438b4df89703560d7c89c9135f 100644
--- a/matlab/lyapunov_solver.m
+++ b/matlab/lyapunov_solver.m
@@ -1,4 +1,4 @@
-function P=lyapunov_solver(T,R,Q,DynareOptions) % --*-- Unitary tests --*--
+function P=lyapunov_solver(T,R,Q,DynareOptions)
 % function P=lyapunov_solver(T,R,Q,DynareOptions)
 % Solves the Lyapunov equation P-T*P*T' = R*Q*R' arising in a state-space
 % system, where P is the variance of the states
@@ -23,7 +23,7 @@ function P=lyapunov_solver(T,R,Q,DynareOptions) % --*-- Unitary tests --*--
 %       Square-root solver for discrete-time Lyapunov equations (requires Matlab System Control toolbox
 %       or Octave control package)
 
-% Copyright © 2016-2022 Dynare Team
+% Copyright © 2016-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -65,122 +65,124 @@ else
     P = lyapunov_symm(T,R*Q*R',DynareOptions.lyapunov_fixed_point_tol,DynareOptions.qz_criterium,DynareOptions.lyapunov_complex_threshold, [], DynareOptions.debug);
 end
 
+return % --*-- Unit tests --*--
+
 %@test:1
-%$ t = NaN(13,1);
-%$ options_.lyapunov_complex_threshold = 1e-15;
-%$ options_.qz_zero_threshold = 1e-6;
-%$ options_.qz_criterium=1-options_.qz_zero_threshold;
-%$ options_.lyapunov_fixed_point_tol = 1e-10;
-%$ options_.lyapunov_doubling_tol = 1e-16;
-%$ options_.debug=false;
-%$
-%$ n_small=8;
-%$ m_small=10;
-%$ T_small=randn(n_small,n_small);
-%$ T_small=0.99*T_small/max(abs(eigs(T_small)));
-%$ tmp2=randn(m_small,m_small);
-%$ Q_small=tmp2*tmp2';
-%$ R_small=randn(n_small,m_small);
-%$
-%$ n_large=9;
-%$ m_large=11;
-%$ T_large=randn(n_large,n_large);
-%$ T_large=0.99*T_large/max(abs(eigs(T_large)));
-%$ tmp2=randn(m_large,m_large);
-%$ Q_large=tmp2*tmp2';
-%$ R_large=randn(n_large,m_large);
-%$
-%$ % DynareOptions.lyapunov_fp == 1
-%$ options_.lyapunov_fp = true;
-%$ try
-%$    Pstar1_small = lyapunov_solver(T_small,R_small,Q_small,options_);
-%$    Pstar1_large = lyapunov_solver(T_large,R_large,Q_large,options_);
-%$    t(1) = 1;
-%$ catch
-%$    t(1) = 0;
-%$ end
-%$
-%$ % Dynareoptions.lyapunov_db == 1
-%$ options_.lyapunov_fp = false;
-%$ options_.lyapunov_db = true;
-%$ try
-%$    Pstar2_small = lyapunov_solver(T_small,R_small,Q_small,options_);
-%$    Pstar2_large = lyapunov_solver(T_large,R_large,Q_large,options_);
-%$    t(2) = 1;
-%$ catch
-%$    t(2) = 0;
-%$ end
-%$
-%$ % Dynareoptions.lyapunov_srs == 1
-%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox'))
-%$     options_.lyapunov_db = false;
-%$     options_.lyapunov_srs = true;
-%$     try
-%$        Pstar3_small = lyapunov_solver(T_small,R_small,Q_small,options_);
-%$        Pstar3_large = lyapunov_solver(T_large,R_large,Q_large,options_);
-%$        t(3) = 1;
-%$     catch
-%$        t(3) = 0;
-%$     end
-%$ else
-%$     t(3) = 1;
-%$ end
-%$
-%$ % Standard
-%$     options_.lyapunov_srs = false;
-%$ try
-%$    Pstar4_small = lyapunov_solver(T_small,R_small,Q_small,options_);
-%$    Pstar4_large = lyapunov_solver(T_large,R_large,Q_large,options_);
-%$    t(4) = 1;
-%$ catch
-%$    t(4) = 0;
-%$ end
-%$
-%$ % Test the results.
-%$
-%$ if max(max(abs(Pstar1_small-Pstar2_small)))>1e-8
-%$    t(5) = 0;
-%$ else
-%$    t(5) = 1;
-%$ end
-%$
-%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox'))
-%$    if max(max(abs(Pstar1_small-Pstar3_small)))>1e-8
-%$       t(6) = 0;
-%$    else
-%$       t(6) = 1;
-%$    end
-%$ else
-%$    t(6) = 1;
-%$ end
-%$
-%$ if max(max(abs(Pstar1_small-Pstar4_small)))>1e-8
-%$    t(7) = 0;
-%$ else
-%$    t(7) = 1;
-%$ end
-%$
-%$ if max(max(abs(Pstar1_large-Pstar2_large)))>2e-8
-%$    t(8) = 0;
-%$ else
-%$    t(8) = 1;
-%$ end
-%$
-%$ if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox'))
-%$    if max(max(abs(Pstar1_large-Pstar3_large)))>1e-8
-%$       t(9) = 0;
-%$    else
-%$       t(9) = 1;
-%$    end
-%$ else
-%$    t(9) = 1;
-%$ end
-%$
-%$ if max(max(abs(Pstar1_large-Pstar4_large)))>2e-8
-%$    t(10) = 0;
-%$ else
-%$    t(10) = 1;
-%$ end
-%$ 
-%$ T = all(t);
+t = NaN(13,1);
+options_.lyapunov_complex_threshold = 1e-15;
+options_.qz_zero_threshold = 1e-6;
+options_.qz_criterium=1-options_.qz_zero_threshold;
+options_.lyapunov_fixed_point_tol = 1e-10;
+options_.lyapunov_doubling_tol = 1e-16;
+options_.debug=false;
+
+n_small=8;
+m_small=10;
+T_small=randn(n_small,n_small);
+T_small=0.99*T_small/max(abs(eigs(T_small)));
+tmp2=randn(m_small,m_small);
+Q_small=tmp2*tmp2';
+R_small=randn(n_small,m_small);
+
+n_large=9;
+m_large=11;
+T_large=randn(n_large,n_large);
+T_large=0.99*T_large/max(abs(eigs(T_large)));
+tmp2=randn(m_large,m_large);
+Q_large=tmp2*tmp2';
+R_large=randn(n_large,m_large);
+
+% DynareOptions.lyapunov_fp == 1
+options_.lyapunov_fp = true;
+try
+   Pstar1_small = lyapunov_solver(T_small,R_small,Q_small,options_);
+   Pstar1_large = lyapunov_solver(T_large,R_large,Q_large,options_);
+   t(1) = 1;
+catch
+   t(1) = 0;
+end
+
+% Dynareoptions.lyapunov_db == 1
+options_.lyapunov_fp = false;
+options_.lyapunov_db = true;
+try
+   Pstar2_small = lyapunov_solver(T_small,R_small,Q_small,options_);
+   Pstar2_large = lyapunov_solver(T_large,R_large,Q_large,options_);
+   t(2) = 1;
+catch
+   t(2) = 0;
+end
+
+% Dynareoptions.lyapunov_srs == 1
+if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox'))
+    options_.lyapunov_db = false;
+    options_.lyapunov_srs = true;
+    try
+       Pstar3_small = lyapunov_solver(T_small,R_small,Q_small,options_);
+       Pstar3_large = lyapunov_solver(T_large,R_large,Q_large,options_);
+       t(3) = 1;
+    catch
+       t(3) = 0;
+    end
+else
+    t(3) = 1;
+end
+
+% Standard
+    options_.lyapunov_srs = false;
+try
+   Pstar4_small = lyapunov_solver(T_small,R_small,Q_small,options_);
+   Pstar4_large = lyapunov_solver(T_large,R_large,Q_large,options_);
+   t(4) = 1;
+catch
+   t(4) = 0;
+end
+
+% Test the results.
+
+if max(max(abs(Pstar1_small-Pstar2_small)))>1e-8
+   t(5) = 0;
+else
+   t(5) = 1;
+end
+
+if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox'))
+   if max(max(abs(Pstar1_small-Pstar3_small)))>1e-8
+      t(6) = 0;
+   else
+      t(6) = 1;
+   end
+else
+   t(6) = 1;
+end
+
+if max(max(abs(Pstar1_small-Pstar4_small)))>1e-8
+   t(7) = 0;
+else
+   t(7) = 1;
+end
+
+if max(max(abs(Pstar1_large-Pstar2_large)))>2e-8
+   t(8) = 0;
+else
+   t(8) = 1;
+end
+
+if (isoctave && user_has_octave_forge_package('control')) || (~isoctave && user_has_matlab_license('control_toolbox'))
+   if max(max(abs(Pstar1_large-Pstar3_large)))>1e-8
+      t(9) = 0;
+   else
+      t(9) = 1;
+   end
+else
+   t(9) = 1;
+end
+
+if max(max(abs(Pstar1_large-Pstar4_large)))>2e-8
+    t(10) = 0;
+else
+    t(10) = 1;
+end
+
+T = all(t);
 %@eof:1
diff --git a/matlab/lyapunov_symm.m b/matlab/lyapunov_symm.m
index cae88bd252a0f6f9780ef067ec5aebb19526a945..e7c4f153cb53f112edc4ee87b218d5d701c45907 100644
--- a/matlab/lyapunov_symm.m
+++ b/matlab/lyapunov_symm.m
@@ -1,4 +1,4 @@
-function [x,u] = lyapunov_symm(a,b,lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,method,debug)  % --*-- Unitary tests --*--
+function [x,u] = lyapunov_symm(a,b,lyapunov_fixed_point_tol,qz_criterium,lyapunov_complex_threshold,method,debug)
 % Solves the Lyapunov equation x-a*x*a' = b, for b and x symmetric matrices.
 % If a has some unit roots, the function computes only the solution of the stable subsystem.
 %
@@ -25,7 +25,7 @@ function [x,u] = lyapunov_symm(a,b,lyapunov_fixed_point_tol,qz_criterium,lyapuno
 % SPECIAL REQUIREMENTS
 %   None
 
-% Copyright © 2006-2017 Dynare Team
+% Copyright © 2006-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -163,4 +163,4 @@ if i == 1
     x(1,1) = (B(1,1)+c)/(1-T(1,1)*T(1,1));
 end
 x = U(:,k+1:end)*x*U(:,k+1:end)';
-u = U(:,1:k);
\ No newline at end of file
+u = U(:,1:k);
diff --git a/matlab/missing/mex/cycle_reduction/cycle_reduction.m b/matlab/missing/mex/cycle_reduction/cycle_reduction.m
index 8f33c04c42031ccd52b2fe95b9ff987dd5f8c603..d3e126de2bd39b02dfad7e59a089c02b88873ad6 100644
--- a/matlab/missing/mex/cycle_reduction/cycle_reduction.m
+++ b/matlab/missing/mex/cycle_reduction/cycle_reduction.m
@@ -1,4 +1,4 @@
-function [X, info] = cycle_reduction(A0, A1, A2, cvg_tol, ch) % --*-- Unitary tests --*--
+function [X, info] = cycle_reduction(A0, A1, A2, cvg_tol, ch)
 
 %@info:
 %! @deftypefn {Function File} {[@var{X}, @var{info}] =} cycle_reduction (@var{A0},@var{A1},@var{A2},@var{cvg_tol},@var{ch})
@@ -43,7 +43,7 @@ function [X, info] = cycle_reduction(A0, A1, A2, cvg_tol, ch) % --*-- Unitary te
 %! @end deftypefn
 %@eod:
 
-% Copyright © 2013-2017 Dynare Team
+% Copyright © 2013-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -112,40 +112,42 @@ if (nargin == 5 && ~isempty(ch) )
     end
 end
 
+return % --*-- Unit tests --*--
+
 %@test:1
-%$
-%$ t = zeros(3,1);
-%$
-%$ % Set the dimension of the problem to be solved.
-%$ n = 500;
-%$
-%$ % Set the equation to be solved
-%$ A = eye(n);
-%$ B = diag(30*ones(n,1)); B(1,1) = 20; B(end,end) = 20; B = B - diag(10*ones(n-1,1),-1); B = B - diag(10*ones(n-1,1),1);
-%$ C = diag(15*ones(n,1)); C = C - diag(5*ones(n-1,1),-1); C = C - diag(5*ones(n-1,1),1);
-%$
-%$ % Solve the equation with the cycle reduction algorithm
-%$ try
-%$     tic; X1 = cycle_reduction(C,B,A,1e-7); elapsedtime = toc;
-%$     disp(['Elapsed time for cycle reduction algorithm is: ' num2str(elapsedtime) ' (n=' int2str(n) ').'])
-%$     t(1) = 1;
-%$ catch
-%$     % nothing to do here.
-%$ end
-%$
-%$ % Solve the equation with the logarithmic reduction algorithm
-%$ try
-%$     tic; X2 = logarithmic_reduction(A,B,C,1e-16,100); elapsedtime = toc;
-%$     disp(['Elapsed time for logarithmic reduction algorithm is: ' num2str(elapsedtime) ' (n=' int2str(n) ').'])
-%$     t(2) = 1;
-%$ catch
-%$     % nothing to do here.
-%$ end
-%$
-%$ % Check the results.
-%$ if t(1) && t(2)
-%$     t(3) = dassert(X1,X2,1e-12);
-%$ end
-%$
-%$ T = all(t);
+
+t = zeros(3,1);
+
+% Set the dimension of the problem to be solved.
+n = 500;
+
+% Set the equation to be solved
+A = eye(n);
+B = diag(30*ones(n,1)); B(1,1) = 20; B(end,end) = 20; B = B - diag(10*ones(n-1,1),-1); B = B - diag(10*ones(n-1,1),1);
+C = diag(15*ones(n,1)); C = C - diag(5*ones(n-1,1),-1); C = C - diag(5*ones(n-1,1),1);
+
+% Solve the equation with the cycle reduction algorithm
+try
+    tic; X1 = cycle_reduction(C,B,A,1e-7); elapsedtime = toc;
+    disp(['Elapsed time for cycle reduction algorithm is: ' num2str(elapsedtime) ' (n=' int2str(n) ').'])
+    t(1) = 1;
+catch
+    % nothing to do here.
+end
+
+% Solve the equation with the logarithmic reduction algorithm
+try
+    tic; X2 = logarithmic_reduction(A,B,C,1e-16,100); elapsedtime = toc;
+    disp(['Elapsed time for logarithmic reduction algorithm is: ' num2str(elapsedtime) ' (n=' int2str(n) ').'])
+    t(2) = 1;
+catch
+    % nothing to do here.
+end
+
+% Check the results.
+if t(1) && t(2)
+    t(3) = dassert(X1,X2,1e-12);
+end
+
+T = all(t);
 %@eof:1
\ No newline at end of file
diff --git a/matlab/missing/mex/local_state_space_iterations/local_state_space_iteration_2.m b/matlab/missing/mex/local_state_space_iterations/local_state_space_iteration_2.m
index 1328cb9ca4434d9139a12435794e7e84e86f0140..ef678fb8c31cf47e2c744d59d02d8768ee1f1f9f 100644
--- a/matlab/missing/mex/local_state_space_iterations/local_state_space_iteration_2.m
+++ b/matlab/missing/mex/local_state_space_iterations/local_state_space_iteration_2.m
@@ -1,4 +1,4 @@
-function [y, y_] = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, a, b, c) % --*-- Unitary tests --*--
+function [y, y_] = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, constant, ghxx, ghuu, ghxu, a, b, c)
 
 % Given the demeaned states (yhat) and structural innovations (epsilon), this routine computes the level of selected endogenous variables when the
 % model is approximated by an order two taylor expansion around the deterministic steady state. Depending on the number of input/output
@@ -26,7 +26,7 @@ function [y, y_] = local_state_space_iteration_2(yhat, epsilon, ghx, ghu, consta
 % 1. If the function has more than nine input arguments (pruning) then it must have two output arguments (otherwise only one input).
 % 2. Ninth input argument is not used, it is here only to have a common interface with the mex version.
 
-% Copyright © 2011-2022 Dynare Team
+% Copyright © 2011-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -77,7 +77,7 @@ else
     end
 end
 
-return
+return % --*-- Unit tests --*--
 
 %@test:1
 n = 2;
diff --git a/matlab/missing/mex/mjdgges/mjdgges.m b/matlab/missing/mex/mjdgges/mjdgges.m
index e163e6caca47873eda2f919a2340cc19e3a7c2a1..85f2b78a44282f74d928a607187ed6fb32d18a6b 100644
--- a/matlab/missing/mex/mjdgges/mjdgges.m
+++ b/matlab/missing/mex/mjdgges/mjdgges.m
@@ -1,4 +1,4 @@
-function [ss, tt, zz, sdim, eigval, info] = mjdgges(e, d, qz_criterium, zhreshold) % --*-- Unitary tests --*--
+function [ss, tt, zz, sdim, eigval, info] = mjdgges(e, d, qz_criterium, zhreshold)
 %
 % INPUTS
 %   e            [double] real square (n*n) matrix.
@@ -18,7 +18,7 @@ function [ss, tt, zz, sdim, eigval, info] = mjdgges(e, d, qz_criterium, zhreshol
 % SPECIAL REQUIREMENTS
 %   none.
 
-% Copyright © 1996-2022 Dynare Team
+% Copyright © 1996-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -66,18 +66,20 @@ catch
     info = 1; % Not as precise as lapack's info!
 end
 
+return % --*-- Unit tests --*--
+
 %@test:1
-%$ try
-%$     E =[0,0,0,0,0,0,0,-1,0,0;0,0,0,0,0,0,0,0,-1,0;0,0,0,0,0,0,0,0,0,-0.990099009900990;0,0,0,0,0,0,0,0,0,0.0990099009900990;0,0,0,-1.01010101010101,0,0.0427672955974843,0,0,0,0;0,0,0,0,0,0.128301886792453,-1,0,0,0;0.800000000000000,0,0,0,0,0,0,0,0,0;0,1,0,0,0,0,1,0,0,0;0,0,0.900000000000000,0,0,0,0,0,0,0;0,0,0,-1.01010101010101,-1,0,2,0,-1,0];
-%$     D=[0,0,0,0,-1,0,0,-0.792000000000000,0,0;0,0,0,0,0,0,0,0,-0.990000000000000,0;0,0,0,-0.000493818030899887,0,0,0,0,0,-0.882178217821782;0,0,0,-1.00493818030900,0,0,0,0,0,0.0882178217821782;0,0,0,-1,0.128301886792453,0,0,0,0,0;-1,0,0,0,0,0,-0.990000000000000,0,0,0;1,0,0,0,0,0,0,0,0,0;0,1,0,0,0,0,0,0,0,0;0,0,1,0,0,0,0,0,0,0;0,0,0,0,-1,0,0,0,0,0];
-%$     [ss, tt, w, sdim, dr.eigval, info1]=mjdgges(E, D, 1.000001, 1e-06);
-%$     if sdim==5
-%$         t(1) = 1;
-%$     else 
-%$         t(1) = 0;
-%$     end
-%$ catch
-%$     t(1) = 0;
-%$ end
-%$ T = all(t);
-%@eof:1
+try
+    E =[0,0,0,0,0,0,0,-1,0,0;0,0,0,0,0,0,0,0,-1,0;0,0,0,0,0,0,0,0,0,-0.990099009900990;0,0,0,0,0,0,0,0,0,0.0990099009900990;0,0,0,-1.01010101010101,0,0.0427672955974843,0,0,0,0;0,0,0,0,0,0.128301886792453,-1,0,0,0;0.800000000000000,0,0,0,0,0,0,0,0,0;0,1,0,0,0,0,1,0,0,0;0,0,0.900000000000000,0,0,0,0,0,0,0;0,0,0,-1.01010101010101,-1,0,2,0,-1,0];
+    D=[0,0,0,0,-1,0,0,-0.792000000000000,0,0;0,0,0,0,0,0,0,0,-0.990000000000000,0;0,0,0,-0.000493818030899887,0,0,0,0,0,-0.882178217821782;0,0,0,-1.00493818030900,0,0,0,0,0,0.0882178217821782;0,0,0,-1,0.128301886792453,0,0,0,0,0;-1,0,0,0,0,0,-0.990000000000000,0,0,0;1,0,0,0,0,0,0,0,0,0;0,1,0,0,0,0,0,0,0,0;0,0,1,0,0,0,0,0,0,0;0,0,0,0,-1,0,0,0,0,0];
+    [ss, tt, w, sdim, dr.eigval, info1]=mjdgges(E, D, 1.000001, 1e-06);
+    if sdim==5
+        t(1) = 1;
+    else 
+        t(1) = 0;
+    end
+catch
+    t(1) = 0;
+end
+T = all(t);
+%@eof:1
\ No newline at end of file
diff --git a/matlab/missing/stats/corr.m b/matlab/missing/stats/corr.m
index 26c27fe3bd91ef10752fcebb7396fb2e36121eb6..6f8cd0b8cb6f2c17809c652db69f33fa290751e9 100644
--- a/matlab/missing/stats/corr.m
+++ b/matlab/missing/stats/corr.m
@@ -1,4 +1,4 @@
-function retval = corr(x, y) % --*-- Unitary tests --*--
+function retval = corr(x, y)
 %@info:
 %! @deftypefn  {Function File} {} corr (@var{x})
 %! @deftypefnx {Function File} {} corr (@var{x}, @var{y})
@@ -39,7 +39,7 @@ function retval = corr(x, y) % --*-- Unitary tests --*--
 % Copyright © 1993-1996 Kurt Hornik
 % Copyright © 1996-2015 John W. Eaton
 % Copyright © 2013-2015 Julien Bect
-% Copyright © 2016-2017 Dynare Team
+% Copyright © 2016-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -110,79 +110,77 @@ c = x'*y;
 s = sx'*sy;
 retval = c./(n * s);
 
-
-end
-
+return % --*-- Unit tests --*--
 
 %@test:1
-%$ x = rand (10);
-%$ cc1 = corr(x);
-%$ cc2 = corr(x, x);
-%$ t(1) = dassert(size(cc1),[10, 10]);
-%$ t(2) = dassert(size (cc2),[10, 10]);
-%$ t(3) = dassert(cc1, cc2, sqrt (eps));
-%$ T = all(t);
+x = rand (10);
+cc1 = corr(x);
+cc2 = corr(x, x);
+t(1) = dassert(size(cc1),[10, 10]);
+t(2) = dassert(size (cc2),[10, 10]);
+t(3) = dassert(cc1, cc2, sqrt (eps));
+T = all(t);
 %@eof:1
 
 %@test:2
-%$ x = [1:3]';
-%$ y = [3:-1:1]';
-%$ t = zeros(3,1);
-%$ t(1) = dassert(corr(x, y), -1, 5*eps);
-%$ t(2) = dassert(corr(x, flipud (y)), 1, 5*eps);
-%$ t(3) = dassert(corr([x, y]), [1 -1; -1 1], 5*eps);
-%$ T = all(t);
+x = [1:3]';
+y = [3:-1:1]';
+t = zeros(3,1);
+t(1) = dassert(corr(x, y), -1, 5*eps);
+t(2) = dassert(corr(x, flipud (y)), 1, 5*eps);
+t(3) = dassert(corr([x, y]), [1 -1; -1 1], 5*eps);
+T = all(t);
 %@eof:2
 
 %@test:3
-%$ if ~isoctave()
-%$   t = zeros(3,1);
-%$   t(1) = dassert(corr(5), NaN);
-%$   t(2) = dassert(corr([1 2 3],5),NaN(3,1));
-%$   t(3) = dassert(corr(5,[1 2 3]),NaN(1,3));
-%$ else
-%$   t = 1;
-%$ end
-%$ T = all(t);
+if ~isoctave()
+    t = zeros(3,1);
+    t(1) = dassert(corr(5), NaN);
+    t(2) = dassert(corr([1 2 3],5),NaN(3,1));
+    t(3) = dassert(corr(5,[1 2 3]),NaN(1,3));
+else
+    t = 1;
+end
+T = all(t);
 %@eof:3
 
 %@test:4
-%$ t = zeros(6,1);
-%$ try
-%$     corr()
-%$     t(1) = false;
-%$ catch
-%$     t(1) = true;
-%$ end
-%$ try
-%$     corr(1, 2, 3)
-%$     t(2) = false;
-%$ catch
-%$     t(2) = true;
-%$ end
-%$ try
-%$     corr([1; 2], ['A', 'B'])
-%$     t(3) = false;
-%$ catch
-%$     t(3) = true;
-%$ end
-%$ try
-%$     corr([1; 2], ['A', 'B'])
-%$     t(4) = false;
-%$ catch
-%$     t(4) = true;
-%$ end
-%$ try
-%$     corr(ones (2,2,2))
-%$     t(5) = false;
-%$ catch
-%$     t(5) = true;
-%$ end
-%$ try
-%$     corr(ones (2,2), ones (2,2,2))
-%$     t(6) = false;
-%$ catch
-%$     t(6) = true;
-%$ end
-%$ T = all(t);
+t = zeros(6,1);
+try
+    corr()
+    t(1) = false;
+catch
+    t(1) = true;
+end
+try
+    corr(1, 2, 3)
+    t(2) = false;
+catch
+    t(2) = true;
+end
+try
+    corr([1; 2], ['A', 'B'])
+    t(3) = false;
+catch
+    t(3) = true;
+end
+try
+    corr([1; 2], ['A', 'B'])
+    t(4) = false;
+catch
+    t(4) = true;
+end
+try
+    corr(ones (2,2,2))
+    t(5) = false;
+catch
+    t(5) = true;
+end
+try
+    corr(ones (2,2), ones (2,2,2))
+    t(6) = false;
+catch
+    t(6) = true;
+end
+T = all(t);
 %@eof:4
\ No newline at end of file
diff --git a/matlab/missing/stats/gamrnd.m b/matlab/missing/stats/gamrnd.m
index da24cc2df8036e5d94293833d8e6aeae1d17b1b6..a791332e052f8802757d91a48ae17dd450bad3be 100644
--- a/matlab/missing/stats/gamrnd.m
+++ b/matlab/missing/stats/gamrnd.m
@@ -1,4 +1,4 @@
-function rnd = gamrnd(a, b, method) % --*-- Unitary tests --*--
+function rnd = gamrnd(a, b, method)
 
 % This function produces independent random variates from the Gamma distribution.
 %
@@ -126,7 +126,7 @@ if ~isempty(ddx)
     end
 end
 
-return
+return % --*-- Unit tests --*--
 
 %@test:1
 if ~isoctave && ~user_has_matlab_license('statistics_toolbox')
diff --git a/matlab/missing/stats/quantile.m b/matlab/missing/stats/quantile.m
index 28738c6faa9385395d4281b2ebf64d8c20f0faf2..54c9832d3e73612b0ad7a09dbff47c2b873e3437 100644
--- a/matlab/missing/stats/quantile.m
+++ b/matlab/missing/stats/quantile.m
@@ -1,4 +1,4 @@
-function [q,N] = quantile(X, p, dim, method, weights) % --*-- Unitary tests --*--
+function [q,N] = quantile(X, p, dim, method, weights)
 
 % Quantiles of a sample via various methods.
 %
@@ -49,7 +49,7 @@ function [q,N] = quantile(X, p, dim, method, weights) % --*-- Unitary tests --*-
 % http://fr.mathworks.com/matlabcentral/fileexchange/46555-quantile-calculation
 %
 % Copyright © 2014-2016 University of Surrey (Christopher Hummersone)
-% Copyright © 2016-2017 Dynare Team
+% Copyright © 2016-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -234,24 +234,25 @@ function y = irearrange(x,order,shape)
 y = reshape(x,shape);
 y = ipermute(y,order);
 
+return % --*-- Unit tests --*--
 
 %@test:1
-%$ X = randn(10000000, 1);
-%$
-%$ try
-%$   q = quantile(X, [.25, .5, .75, .95 ]);
-%$   t(1) = true;
-%$ catch
-%$   t(1) = false;
-%$ end
-%$
-%$ e = [-0.674489750196082, 0, 0.674489750196082, 1.644853626951472];
-%$
-%$ if t(1)
-%$    for i=1:4
-%$        t(i+1) = abs(q(i)-e(i))<2e-3;
-%$    end
-%$ end
-%$
-%$ T = all(t);
-%@eof:1
\ No newline at end of file
+X = randn(10000000, 1);
+
+try
+  q = quantile(X, [.25, .5, .75, .95 ]);
+  t(1) = true;
+catch
+  t(1) = false;
+end
+
+e = [-0.674489750196082, 0, 0.674489750196082, 1.644853626951472];
+
+if t(1)
+    for i=1:4
+        t(i+1) = abs(q(i)-e(i))<2e-3;
+    end
+end
+
+T = all(t);
+%@eof:1
diff --git a/matlab/missing/stats/wblcdf.m b/matlab/missing/stats/wblcdf.m
index 7fe25e36c6657bec328fdfa300e0db09786eb6c9..20d0aa5d77737729dcd4f4542b111e3e6d4d73c3 100644
--- a/matlab/missing/stats/wblcdf.m
+++ b/matlab/missing/stats/wblcdf.m
@@ -1,4 +1,4 @@
-function p = wblcdf(x, scale, shape)   % --*-- Unitary tests --*--
+function p = wblcdf(x, scale, shape)
 
 % Cumulative distribution function for the Weibull distribution.
 %
@@ -10,7 +10,7 @@ function p = wblcdf(x, scale, shape)   % --*-- Unitary tests --*--
 % OUTPUTS
 % - p     [double] Positive scalar between
 
-% Copyright © 2015-2017 Dynare Team
+% Copyright © 2015-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -63,86 +63,86 @@ end
 p = 1-exp(-(x/scale)^shape);
 
 %@test:1
-%$ try
-%$     p = wblcdf(-1, .5, .1);
-%$     t(1) = true;
-%$ catch
-%$     t(1) = false;
-%$ end
-%$
-%$ % Check the results
-%$ if t(1)
-%$     t(2) = isequal(p, 0);
-%$ end
-%$ T = all(t);
+try
+    p = wblcdf(-1, .5, .1);
+    t(1) = true;
+catch
+    t(1) = false;
+end
+
+% Check the results
+if t(1)
+    t(2) = isequal(p, 0);
+end
+T = all(t);
 %@eof:1
 
 %@test:2
-%$ try
-%$     p = wblcdf(Inf, .5, .1);
-%$     t(1) = true;
-%$ catch
-%$     t(1) = false;
-%$ end
-%$
-%$ % Check the results
-%$ if t(1)
-%$     t(2) = isequal(p, 1);
-%$ end
-%$ T = all(t);
+try
+    p = wblcdf(Inf, .5, .1);
+    t(1) = true;
+catch
+    t(1) = false;
+end
+
+% Check the results
+if t(1)
+    t(2) = isequal(p, 1);
+end
+T = all(t);
 %@eof:2
 
 %@test:3
-%$ % Set the hyperparameters of a Weibull definition.
-%$ scale = .5;
-%$ shape = 1.5;
-%$
-%$ % Compute the median of the weibull distribution.
-%$ m = scale*log(2)^(1/shape);
-%$
-%$ try
-%$     p = wblcdf(m, scale, shape);
-%$     t(1) = true;
-%$ catch
-%$     t(1) = false;
-%$ end
-%$
-%$ % Check the results
-%$ if t(1)
-%$     t(2) = abs(p-.5)<1e-12;
-%$ end
-%$ T = all(t);
+% Set the hyperparameters of a Weibull definition.
+scale = .5;
+shape = 1.5;
+
+% Compute the median of the weibull distribution.
+m = scale*log(2)^(1/shape);
+
+try
+    p = wblcdf(m, scale, shape);
+    t(1) = true;
+catch
+    t(1) = false;
+end
+
+% Check the results
+if t(1)
+    t(2) = abs(p-.5)<1e-12;
+end
+T = all(t);
 %@eof:3
 
 %@test:4
-%$ % Consistency check between wblinv and wblcdf.
-%$
-%$ % Set the hyperparameters of a Weibull definition.
-%$ scale = .5;
-%$ shape = 1.5;
-%$
-%$ % Compute quatiles of the weibull distribution.
-%$ q = 0:.05:1;
-%$ m = zeros(size(q));
-%$ p = zeros(size(q));
-%$ for i=1:length(q)
-%$     m(i) = wblinv(q(i), scale, shape);
-%$ end
-%$
-%$ try
-%$     for i=1:length(q)
-%$         p(i) = wblcdf(m(i), scale, shape);
-%$     end
-%$     t(1) = true;
-%$ catch
-%$     t(1) = false;
-%$ end
-%$
-%$ % Check the results
-%$ if t(1)
-%$     for i=1:length(q)
-%$         t(i+1) = abs(p(i)-q(i))<1e-12;
-%$     end
-%$ end
-%$ T = all(t);
+% Consistency check between wblinv and wblcdf.
+
+% Set the hyperparameters of a Weibull definition.
+scale = .5;
+shape = 1.5;
+
+% Compute quatiles of the weibull distribution.
+q = 0:.05:1;
+m = zeros(size(q));
+p = zeros(size(q));
+for i=1:length(q)
+    m(i) = wblinv(q(i), scale, shape);
+end
+
+try
+    for i=1:length(q)
+        p(i) = wblcdf(m(i), scale, shape);
+    end
+    t(1) = true;
+catch
+    t(1) = false;
+end
+
+% Check the results
+if t(1)
+    for i=1:length(q)
+        t(i+1) = abs(p(i)-q(i))<1e-12;
+    end
+end
+T = all(t);
 %@eof:4
diff --git a/matlab/missing/stats/wblinv.m b/matlab/missing/stats/wblinv.m
index 8742865298d33635d3978b55659efaa20e31da8d..aee1f1706ea5f505a68323bd1c2dcddee0fff30c 100644
--- a/matlab/missing/stats/wblinv.m
+++ b/matlab/missing/stats/wblinv.m
@@ -1,4 +1,4 @@
-function t = wblinv(proba, scale, shape)   % --*-- Unitary tests --*--
+function t = wblinv(proba, scale, shape)
 
 % Inverse cumulative distribution function.
 %
@@ -10,7 +10,7 @@ function t = wblinv(proba, scale, shape)   % --*-- Unitary tests --*--
 % OUTPUTS
 % - t     [double] scalar such that P(X<=t)=proba
 
-% Copyright © 2015-2020 Dynare Team
+% Copyright © 2015-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -58,110 +58,112 @@ end
 
 t = exp(log(scale)+log(-log(1-proba))/shape);
 
+return % --*-- Unit tests --*--
+
 %@test:1
-%$ try
-%$    x = wblinv(0, 1, 2);
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = isequal(x, 0);
-%$ end
-%$ T = all(t);
+try
+   x = wblinv(0, 1, 2);
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = isequal(x, 0);
+end
+T = all(t);
 %@eof:1
 
 %@test:2
-%$ try
-%$    x = wblinv(1, 1, 2);
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    t(2) = isinf(x);
-%$ end
-%$ T = all(t);
+try
+   x = wblinv(1, 1, 2);
+   t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = isinf(x);
+end
+T = all(t);
 %@eof:2
 
 %@test:3
-%$ scales = [.5, 1, 5];
-%$ shapes = [.1, 1, 2];
-%$ x = NaN(9,1);
-%$
-%$ try
-%$    k = 0;
-%$    for i=1:3
-%$       for j=1:3
-%$           k = k+1;
-%$           x(k) = wblinv(.5, scales(i), shapes(j));
-%$       end
-%$    end
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    k = 1;
-%$    for i=1:3
-%$       for j=1:3
-%$           k = k+1;
-%$           t(k) = abs(x(k-1)-scales(i)*log(2)^(1/shapes(j)))<1e-12;
-%$       end
-%$    end
-%$ end
-%$ T = all(t);
+scales = [.5, 1, 5];
+shapes = [.1, 1, 2];
+x = NaN(9,1);
+
+try
+   k = 0;
+   for i=1:3
+      for j=1:3
+          k = k+1;
+          x(k) = wblinv(.5, scales(i), shapes(j));
+      end
+   end
+   t(1) = true;
+catch
+   t(1) = false;
+end
+
+if t(1)
+    k = 1;
+    for i=1:3
+        for j=1:3
+            k = k+1;
+            t(k) = abs(x(k-1)-scales(i)*log(2)^(1/shapes(j)))<1e-12;
+        end
+    end
+end
+T = all(t);
 %@eof:3
 
 %@test:4
-%$ debug = false;
-%$ scales = [ .5, 1, 5];
-%$ shapes = [ 1, 2, 3];
-%$ x = NaN(9,1);
-%$ p = 1e-1;
-%$
-%$ try
-%$    k = 0;
-%$    for i=1:3
-%$       for j=1:3
-%$           k = k+1;
-%$           x(k) = wblinv(p, scales(i), shapes(j));
-%$       end
-%$    end
-%$    t(1) = true;
-%$ catch
-%$    t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$    k = 1;
-%$    for i=1:3
-%$       for j=1:3
-%$           k = k+1;
-%$           shape = shapes(j);
-%$           scale = scales(i);
-%$           density = @(z) exp(lpdfgweibull(z,shape,scale));
-%$           if debug
-%$               [shape, scale, x(k-1)]
-%$           end
-%$           if isoctave
-%$               s = quadv(density, 0, x(k-1),1e-10);
-%$           else
-%$               s = integral(density, 0, x(k-1));
-%$           end
-%$           if debug
-%$               [s, abs(p-s)]
-%$           end
-%$         if isoctave
-%$           t(k) = abs(p-s)<1e-9;
-%$         else
-%$           t(k) = abs(p-s)<1e-12;
-%$         end
-%$       end
-%$    end
-%$ end
-%$ T = all(t);
+debug = false;
+scales = [ .5, 1, 5];
+shapes = [ 1, 2, 3];
+x = NaN(9,1);
+p = 1e-1;
+
+try
+   k = 0;
+   for i=1:3
+      for j=1:3
+          k = k+1;
+          x(k) = wblinv(p, scales(i), shapes(j));
+      end
+   end
+   t(1) = true;
+catch
+   t(1) = false;
+end
+
+if t(1)
+   k = 1;
+   for i=1:3
+      for j=1:3
+          k = k+1;
+          shape = shapes(j);
+          scale = scales(i);
+          density = @(z) exp(lpdfgweibull(z,shape,scale));
+          if debug
+              [shape, scale, x(k-1)]
+          end
+          if isoctave
+              s = quadv(density, 0, x(k-1),1e-10);
+          else
+              s = integral(density, 0, x(k-1));
+          end
+          if debug
+              [s, abs(p-s)]
+          end
+          if isoctave
+              t(k) = abs(p-s)<1e-9;
+          else
+              t(k) = abs(p-s)<1e-12;
+          end
+      end
+   end
+end
+T = all(t);
 %@eof:4
diff --git a/matlab/modules/dseries b/matlab/modules/dseries
index a519e8a68177e2a5beccff3d04ddcb9e4f64030d..f2960bf774113d278eeb43faf9a30c31495e6843 160000
--- a/matlab/modules/dseries
+++ b/matlab/modules/dseries
@@ -1 +1 @@
-Subproject commit a519e8a68177e2a5beccff3d04ddcb9e4f64030d
+Subproject commit f2960bf774113d278eeb43faf9a30c31495e6843
diff --git a/matlab/prior_draw.m b/matlab/prior_draw.m
index e5d244b388a90475754e87064b4d5c789e7436ce..d5e0f436b107ded33bf1d6776405a2bedccae31f 100644
--- a/matlab/prior_draw.m
+++ b/matlab/prior_draw.m
@@ -1,4 +1,4 @@
-function pdraw = prior_draw(BayesInfo, prior_trunc, uniform) % --*-- Unitary tests --*--
+function pdraw = prior_draw(BayesInfo, prior_trunc, uniform)
 
 % This function generate one draw from the joint prior distribution and
 % allows sampling uniformly from the prior support (uniform==1 when called with init==1)
@@ -26,7 +26,7 @@ function pdraw = prior_draw(BayesInfo, prior_trunc, uniform) % --*-- Unitary tes
 % NOTE 3. This code relies on bayestopt_ as created in the base workspace
 %           by the preprocessor (or as updated in subsequent pieces of code and handed to the base workspace)
 %
-% Copyright © 2006-2017 Dynare Team
+% Copyright © 2006-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -174,101 +174,103 @@ if weibull_draws
     end
 end
 
+return % --*-- Unit tests --*--
+
 %@test:1
-%$ % Fill global structures with required fields...
-%$ prior_trunc = 1e-10;
-%$ p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1);    % Prior shape
-%$ p1 = .4*ones(14,1);                          % Prior mean
-%$ p2 = .2*ones(14,1);                          % Prior std.
-%$ p3 = NaN(14,1);
-%$ p4 = NaN(14,1);
-%$ p5 = NaN(14,1);
-%$ p6 = NaN(14,1);
-%$ p7 = NaN(14,1);
-%$
-%$ for i=1:14
-%$     switch p0(i)
-%$       case 1
-%$         % Beta distribution
-%$         p3(i) = 0;
-%$         p4(i) = 1;
-%$         [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
-%$         p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
-%$       case 2
-%$         % Gamma distribution
-%$         p3(i) = 0;
-%$         p4(i) = Inf;
-%$         [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
-%$         p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
-%$       case 3
-%$         % Normal distribution
-%$         p3(i) = -Inf;
-%$         p4(i) = Inf;
-%$         p6(i) = p1(i);
-%$         p7(i) = p2(i);
-%$         p5(i) = p1(i);
-%$       case 4
-%$         % Inverse Gamma (type I) distribution
-%$         p3(i) = 0;
-%$         p4(i) = Inf;
-%$         [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
-%$         p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
-%$       case 5
-%$         % Uniform distribution
-%$         [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
-%$         p3(i) = p6(i);
-%$         p4(i) = p7(i);
-%$         p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
-%$       case 6
-%$         % Inverse Gamma (type II) distribution
-%$         p3(i) = 0;
-%$         p4(i) = Inf;
-%$         [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
-%$         p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
-%$       case 8
-%$         % Weibull distribution
-%$         p3(i) = 0;
-%$         p4(i) = Inf;
-%$         [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
-%$         p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
-%$       otherwise
-%$         error('This density is not implemented!')
-%$     end
-%$ end
-%$
-%$ BayesInfo.pshape = p0;
-%$ BayesInfo.p1 = p1;
-%$ BayesInfo.p2 = p2;
-%$ BayesInfo.p3 = p3;
-%$ BayesInfo.p4 = p4;
-%$ BayesInfo.p5 = p5;
-%$ BayesInfo.p6 = p6;
-%$ BayesInfo.p7 = p7;
-%$
-%$ ndraws = 1e5;
-%$ m0 = BayesInfo.p1; %zeros(14,1);
-%$ v0 = diag(BayesInfo.p2.^2); %zeros(14);
-%$
-%$ % Call the tested routine
-%$ try
-%$     % Initialization of prior_draws.
-%$     prior_draw(BayesInfo, prior_trunc, false);
-%$     % Do simulations in a loop and estimate recursively the mean and teh variance.
-%$     for i = 1:ndraws
-%$          draw = transpose(prior_draw());
-%$          m1 = m0 + (draw-m0)/i;
-%$          m2 = m1*m1';
-%$          v0 = v0 + ((draw*draw'-m2-v0) + (i-1)*(m0*m0'-m2'))/i;
-%$          m0 = m1;
-%$     end
-%$     t(1) = true;
-%$ catch
-%$     t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$     t(2) = all(abs(m0-BayesInfo.p1)<3e-3);
-%$     t(3) = all(all(abs(v0-diag(BayesInfo.p2.^2))<5e-3));
-%$ end
-%$ T = all(t);
+% Fill global structures with required fields...
+prior_trunc = 1e-10;
+p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1);    % Prior shape
+p1 = .4*ones(14,1);                          % Prior mean
+p2 = .2*ones(14,1);                          % Prior std.
+p3 = NaN(14,1);
+p4 = NaN(14,1);
+p5 = NaN(14,1);
+p6 = NaN(14,1);
+p7 = NaN(14,1);
+
+for i=1:14
+    switch p0(i)
+      case 1
+        % Beta distribution
+        p3(i) = 0;
+        p4(i) = 1;
+        [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
+        p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
+      case 2
+        % Gamma distribution
+        p3(i) = 0;
+        p4(i) = Inf;
+        [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
+        p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
+      case 3
+        % Normal distribution
+        p3(i) = -Inf;
+        p4(i) = Inf;
+        p6(i) = p1(i);
+        p7(i) = p2(i);
+        p5(i) = p1(i);
+      case 4
+        % Inverse Gamma (type I) distribution
+        p3(i) = 0;
+        p4(i) = Inf;
+        [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
+        p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
+      case 5
+        % Uniform distribution
+        [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
+        p3(i) = p6(i);
+        p4(i) = p7(i);
+        p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
+      case 6
+        % Inverse Gamma (type II) distribution
+        p3(i) = 0;
+        p4(i) = Inf;
+        [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
+        p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
+      case 8
+        % Weibull distribution
+        p3(i) = 0;
+        p4(i) = Inf;
+        [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
+        p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
+      otherwise
+        error('This density is not implemented!')
+    end
+end
+
+BayesInfo.pshape = p0;
+BayesInfo.p1 = p1;
+BayesInfo.p2 = p2;
+BayesInfo.p3 = p3;
+BayesInfo.p4 = p4;
+BayesInfo.p5 = p5;
+BayesInfo.p6 = p6;
+BayesInfo.p7 = p7;
+
+ndraws = 1e5;
+m0 = BayesInfo.p1; %zeros(14,1);
+v0 = diag(BayesInfo.p2.^2); %zeros(14);
+
+% Call the tested routine
+try
+    % Initialization of prior_draws.
+    prior_draw(BayesInfo, prior_trunc, false);
+    % Do simulations in a loop and estimate recursively the mean and teh variance.
+    for i = 1:ndraws
+         draw = transpose(prior_draw());
+         m1 = m0 + (draw-m0)/i;
+         m2 = m1*m1';
+         v0 = v0 + ((draw*draw'-m2-v0) + (i-1)*(m0*m0'-m2'))/i;
+         m0 = m1;
+    end
+    t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = all(abs(m0-BayesInfo.p1)<3e-3);
+    t(3) = all(all(abs(v0-diag(BayesInfo.p2.^2))<5e-3));
+end
+T = all(t);
 %@eof:1
diff --git a/matlab/priordens.m b/matlab/priordens.m
index e432403006daf5b82ceef1f23b7d89bc1b3980b6..f68376794b55ef4e49dcad1e24fa2c43ba59bbcc 100644
--- a/matlab/priordens.m
+++ b/matlab/priordens.m
@@ -1,4 +1,4 @@
-function [logged_prior_density, dlprior, d2lprior, info] = priordens(x, pshape, p6, p7, p3, p4, initialization) % --*-- Unitary tests --*--
+function [logged_prior_density, dlprior, d2lprior, info] = priordens(x, pshape, p6, p7, p3, p4, initialization)
 % Computes a prior density for the structural parameters of DSGE models
 %
 % INPUTS
@@ -15,7 +15,7 @@ function [logged_prior_density, dlprior, d2lprior, info] = priordens(x, pshape,
 %    info                  [double]  error code for index of Inf-prior parameter
 %
 
-% Copyright © 2003-2017 Dynare Team
+% Copyright © 2003-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -191,85 +191,87 @@ if nargout==3
     d2lprior = diag(d2lprior);
 end
 
+return % --*-- Unit tests --*--
+
 %@test:1
-%$ % Fill global structures with required fields...
-%$ prior_trunc = 1e-10;
-%$ p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1);    % Prior shape
-%$ p1 = .4*ones(14,1);                          % Prior mean
-%$ p2 = .2*ones(14,1);                          % Prior std.
-%$ p3 = NaN(14,1);
-%$ p4 = NaN(14,1);
-%$ p5 = NaN(14,1);
-%$ p6 = NaN(14,1);
-%$ p7 = NaN(14,1);
-%$
-%$ for i=1:14
-%$     switch p0(i)
-%$       case 1
-%$         % Beta distribution
-%$         p3(i) = 0;
-%$         p4(i) = 1;
-%$         [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
-%$         p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
-%$       case 2
-%$         % Gamma distribution
-%$         p3(i) = 0;
-%$         p4(i) = Inf;
-%$         [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
-%$         p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
-%$       case 3
-%$         % Normal distribution
-%$         p3(i) = -Inf;
-%$         p4(i) = Inf;
-%$         p6(i) = p1(i);
-%$         p7(i) = p2(i);
-%$         p5(i) = p1(i);
-%$       case 4
-%$         % Inverse Gamma (type I) distribution
-%$         p3(i) = 0;
-%$         p4(i) = Inf;
-%$         [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
-%$         p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
-%$       case 5
-%$         % Uniform distribution
-%$         [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
-%$         p3(i) = p6(i);
-%$         p4(i) = p7(i);
-%$         p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
-%$       case 6
-%$         % Inverse Gamma (type II) distribution
-%$         p3(i) = 0;
-%$         p4(i) = Inf;
-%$         [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
-%$         p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
-%$       case 8
-%$         % Weibull distribution
-%$         p3(i) = 0;
-%$         p4(i) = Inf;
-%$         [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
-%$         p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
-%$       otherwise
-%$         error('This density is not implemented!')
-%$     end
-%$ end
-%$
-%$ % Call the tested routine
-%$ try
-%$     % Initialization of priordens.
-%$     lpdstar = priordens(p5, p0, p6, p7, p3, p4, true);
-%$     % Do simulations in a loop and estimate recursively the mean and teh variance.
-%$     LPD = NaN(10000,1);
-%$     for i = 1:10000
-%$          draw = p5+randn(size(p5))*.02;
-%$          LPD(i) = priordens(p5, p0, p6, p7, p3, p4);
-%$     end
-%$     t(1) = true;
-%$ catch
-%$     t(1) = false;
-%$ end
-%$
-%$ if t(1)
-%$     t(2) = all(LPD<=lpdstar);
-%$ end
-%$ T = all(t);
-%@eof:1
\ No newline at end of file
+% Fill global structures with required fields...
+prior_trunc = 1e-10;
+p0 = repmat([1; 2; 3; 4; 5; 6; 8], 2, 1);    % Prior shape
+p1 = .4*ones(14,1);                          % Prior mean
+p2 = .2*ones(14,1);                          % Prior std.
+p3 = NaN(14,1);
+p4 = NaN(14,1);
+p5 = NaN(14,1);
+p6 = NaN(14,1);
+p7 = NaN(14,1);
+
+for i=1:14
+    switch p0(i)
+      case 1
+        % Beta distribution
+        p3(i) = 0;
+        p4(i) = 1;
+        [p6(i), p7(i)] = beta_specification(p1(i), p2(i)^2, p3(i), p4(i));
+        p5(i) = compute_prior_mode([p6(i) p7(i)], 1);
+      case 2
+        % Gamma distribution
+        p3(i) = 0;
+        p4(i) = Inf;
+        [p6(i), p7(i)] = gamma_specification(p1(i), p2(i)^2, p3(i), p4(i));
+        p5(i) = compute_prior_mode([p6(i) p7(i)], 2);
+      case 3
+        % Normal distribution
+        p3(i) = -Inf;
+        p4(i) = Inf;
+        p6(i) = p1(i);
+        p7(i) = p2(i);
+        p5(i) = p1(i);
+      case 4
+        % Inverse Gamma (type I) distribution
+        p3(i) = 0;
+        p4(i) = Inf;
+        [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 1, false);
+        p5(i) = compute_prior_mode([p6(i) p7(i)], 4);
+      case 5
+        % Uniform distribution
+        [p1(i), p2(i), p6(i), p7(i)] = uniform_specification(p1(i), p2(i), p3(i), p4(i));
+        p3(i) = p6(i);
+        p4(i) = p7(i);
+        p5(i) = compute_prior_mode([p6(i) p7(i)], 5);
+      case 6
+        % Inverse Gamma (type II) distribution
+        p3(i) = 0;
+        p4(i) = Inf;
+        [p6(i), p7(i)] = inverse_gamma_specification(p1(i), p2(i)^2, p3(i), 2, false);
+        p5(i) = compute_prior_mode([p6(i) p7(i)], 6);
+      case 8
+        % Weibull distribution
+        p3(i) = 0;
+        p4(i) = Inf;
+        [p6(i), p7(i)] = weibull_specification(p1(i), p2(i)^2, p3(i));
+        p5(i) = compute_prior_mode([p6(i) p7(i)], 8);
+      otherwise
+        error('This density is not implemented!')
+    end
+end
+
+% Call the tested routine
+try
+    % Initialization of priordens.
+    lpdstar = priordens(p5, p0, p6, p7, p3, p4, true);
+    % Do simulations in a loop and estimate recursively the mean and teh variance.
+    LPD = NaN(10000,1);
+    for i = 1:10000
+         draw = p5+randn(size(p5))*.02;
+         LPD(i) = priordens(p5, p0, p6, p7, p3, p4);
+    end
+    t(1) = true;
+catch
+    t(1) = false;
+end
+
+if t(1)
+    t(2) = all(LPD<=lpdstar);
+end
+T = all(t);
+%@eof:1
diff --git a/matlab/utilities/general/compare_vectors.m b/matlab/utilities/general/compare_vectors.m
index 1e548c854373526d29d6f51f2538e20532fd8697..0e9b9a4694850fb66bd63068291aed4088b93171 100644
--- a/matlab/utilities/general/compare_vectors.m
+++ b/matlab/utilities/general/compare_vectors.m
@@ -1,4 +1,4 @@
-function C = compare_vectors(f, A, B)  % --*-- Unitary tests --*--
+function C = compare_vectors(f, A, B)
 
 % Performs lexicographical comparison of vectors.
 %
@@ -13,7 +13,7 @@ function C = compare_vectors(f, A, B)  % --*-- Unitary tests --*--
 % REMARKS
 %  o It is assumed that vectors A and B have the same number of elements.
 
-% Copyright © 2013-2017 Dynare Team
+% Copyright © 2013-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -76,19 +76,21 @@ else
     end
 end
 
+return % --*-- Unit tests --*--
+
 %@test:1
-%$ t(1) = dassert(compare_vectors(@lt, [1990 3], [1991 3]), 1);
-%$ t(2) = dassert(compare_vectors(@lt, [1990 3], [1990 3]), 0);
-%$ t(3) = dassert(compare_vectors(@le, [1990 3], [1990 3]), 1);
-%$ t(4) = dassert(compare_vectors(@lt, [1990 3], [1990 4]), 1);
-%$ t(5) = dassert(compare_vectors(@le, [1990 3], [1990 4]), 1);
-%$ t(6) = dassert(compare_vectors(@gt, [1990 3], [1991 3]), 0);
-%$ t(7) = dassert(compare_vectors(@gt, [1990 3], [1990 3]), 0);
-%$ t(8) = dassert(compare_vectors(@ge, [1990 3], [1990 3]), 1);
-%$ t(9) = dassert(compare_vectors(@gt, [1990 3], [1990 4]), 0);
-%$ t(10) = dassert(compare_vectors(@ge, [1990 3], [1990 4]), 0);
-%$ t(11) = dassert(compare_vectors(@le, [1991 3], [1990 4]), 0);
-%$ t(12) = dassert(compare_vectors(@le, [1991 3], [1990 2]), 0);
-%$ t(13) = dassert(compare_vectors(@le, [1945 2], [1950, 1]),1);
-%$ T = all(t);
-%@eof:1
\ No newline at end of file
+t(1) = dassert(compare_vectors(@lt, [1990 3], [1991 3]), 1);
+t(2) = dassert(compare_vectors(@lt, [1990 3], [1990 3]), 0);
+t(3) = dassert(compare_vectors(@le, [1990 3], [1990 3]), 1);
+t(4) = dassert(compare_vectors(@lt, [1990 3], [1990 4]), 1);
+t(5) = dassert(compare_vectors(@le, [1990 3], [1990 4]), 1);
+t(6) = dassert(compare_vectors(@gt, [1990 3], [1991 3]), 0);
+t(7) = dassert(compare_vectors(@gt, [1990 3], [1990 3]), 0);
+t(8) = dassert(compare_vectors(@ge, [1990 3], [1990 3]), 1);
+t(9) = dassert(compare_vectors(@gt, [1990 3], [1990 4]), 0);
+t(10) = dassert(compare_vectors(@ge, [1990 3], [1990 4]), 0);
+t(11) = dassert(compare_vectors(@le, [1991 3], [1990 4]), 0);
+t(12) = dassert(compare_vectors(@le, [1991 3], [1990 2]), 0);
+t(13) = dassert(compare_vectors(@le, [1945 2], [1950, 1]),1);
+T = all(t);
+%@eof:1
diff --git a/matlab/utilities/tests b/matlab/utilities/tests
index 343d9e13654a07c348fad863f7f79a03896c6bc9..3c4079a8e1e495be50be3cd81855eb37189fceab 160000
--- a/matlab/utilities/tests
+++ b/matlab/utilities/tests
@@ -1 +1 @@
-Subproject commit 343d9e13654a07c348fad863f7f79a03896c6bc9
+Subproject commit 3c4079a8e1e495be50be3cd81855eb37189fceab
diff --git a/matlab/utilities/version/ver_greater_than.m b/matlab/utilities/version/ver_greater_than.m
index 8130daf8fa7b9f9db62567dba6750d35838e6825..447b4b8a624d1d045419a16a11b9ae78e1cff8c1 100644
--- a/matlab/utilities/version/ver_greater_than.m
+++ b/matlab/utilities/version/ver_greater_than.m
@@ -1,4 +1,4 @@
-function tf = ver_greater_than(ver1, ver2)  % --*-- Unitary tests --*--
+function tf = ver_greater_than(ver1, ver2)
 %function tf = ver_greater_than(ver1, ver2)
 % ver1 > ver2 ? 1 : 0;
 %
@@ -12,7 +12,7 @@ function tf = ver_greater_than(ver1, ver2)  % --*-- Unitary tests --*--
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright © 2015-2021 Dynare Team
+% Copyright © 2015-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -30,7 +30,8 @@ function tf = ver_greater_than(ver1, ver2)  % --*-- Unitary tests --*--
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 tf = ~ver_less_than(ver1, ver2) && ~strcmp(ver1, ver2);
-return
+
+return % --*-- Unit tests --*--
 
 %@test:1
 ver2='4.4';
@@ -38,18 +39,21 @@ ver1='4.5.2';
 t(1)=dassert(ver_greater_than(ver1,ver2),true);
 T = all(t);
 %@eof:1
+
 %@test:2
 ver1='6-unstable-2021-12-15-1737-21a8a579';
 ver2='4.4';
 t(1)=dassert(ver_greater_than(ver1,ver2),true);
 T = all(t);
 %@eof:2
+
 %@test:3
 ver2='5.0';
 ver1='5.1';
 t(1)=dassert(ver_greater_than(ver1,ver2),true);
 T = all(t);
 %@eof:3
+
 %@test:4
 ver2='6-unstable-2021-12-18-1227-c43777f6';
 ver1='6-unstable-2021-12-19-1953-d841fc7c';
diff --git a/matlab/utilities/version/ver_less_than.m b/matlab/utilities/version/ver_less_than.m
index 16676c76df0a47ea5ac11a2e618cfa236820b914..76a3dcef2f056403fada1ecdad15a5489b18819f 100644
--- a/matlab/utilities/version/ver_less_than.m
+++ b/matlab/utilities/version/ver_less_than.m
@@ -1,4 +1,4 @@
-function tf = ver_less_than(ver1, ver2) % --*-- Unitary tests --*--
+function tf = ver_less_than(ver1, ver2)
 %function tf = ver_less_than(ver1, ver2)
 % ver1 < ver2 ? 1 : 0;
 %
@@ -54,7 +54,7 @@ elseif (maj_ver1 == maj_ver2) && (min_ver1 > min_ver2)
 end
 
 %deal with revision in Dynare 4 and unstable versions involved
-if min(maj_ver1,maj_ver2)<5 %old versioning scheme with three digits 
+if min(maj_ver1,maj_ver2)<5 %old versioning scheme with three digits
     if (length(ver1) == length(ver2) && length(ver1) == 3)
         %check if master branch (unstable) or stable
         ismaster1 = isnan(str2double(ver1{3}));
@@ -63,7 +63,7 @@ if min(maj_ver1,maj_ver2)<5 %old versioning scheme with three digits
             %ver2 is the unstable
             return
         end
-        
+
         if ~ismaster1 && ~ismaster2 %both are stable versions
             rev_ver1 = str2double(ver1{3});
             rev_ver2 = str2double(ver2{3});
@@ -76,7 +76,7 @@ if min(maj_ver1,maj_ver2)<5 %old versioning scheme with three digits
         %ver1 is an unstable version
         error('Case is undefined, please contact the developers')
     end
-elseif min(maj_ver1,maj_ver2)>=5 %new versioning scheme with three digits 
+elseif min(maj_ver1,maj_ver2)>=5 %new versioning scheme with three digits
     if strcmp(ver1{2},'x') || strcmp(ver1{2},'unstable')
         date_number_1=datenum([ver1{3} ver1{4} ver1{5}],'YYYYMMDD');
         stable_version_indicator_1=0;
@@ -102,7 +102,7 @@ elseif min(maj_ver1,maj_ver2)>=5 %new versioning scheme with three digits
 end
 tf = false;
 
-return
+return % --*-- Unit tests --*--
 
 %@test:1
 ver1='4.4';
@@ -110,26 +110,24 @@ ver2='4.5.2';
 t(1)=dassert(ver_less_than(ver1,ver2),true);
 T = all(t);
 %@eof:1
+
 %@test:2
 ver1='4.4';
 ver2='6-unstable-2021-12-15-1737-21a8a579';
 t(1)=dassert(ver_less_than(ver1,ver2),true);
 T = all(t);
 %@eof:2
+
 %@test:3
 ver1='5.0';
 ver2='5.1';
 t(1)=dassert(ver_less_than(ver1,ver2),true);
 T = all(t);
 %@eof:3
+
 %@test:4
 ver1='6-unstable-2021-12-18-1227-c43777f6';
 ver2='6-unstable-2021-12-19-1953-d841fc7c';
 t(1)=dassert(ver_less_than(ver1,ver2),true);
 T = all(t);
 %@eof:4
-% %@test:5
-% ver1='5.0';
-% ver2='5.x-2021-12-14-1101-25c1e0c0';
-% t(5)=dassert(ver_less_than(ver1,ver2),true)
-
diff --git a/tests/Makefile.am b/tests/Makefile.am
index e0f8aead52e2a8b3573ee6a1da6bf0a48fc43d16..fef87573981524dd4388040be853162042a07674 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -1259,7 +1259,7 @@ M_TRS_FILES = $(patsubst %.mod, %.m.trs, $(MODFILES))
 M_TRS_FILES += $(patsubst %.mod, %.m.trs, $(PARTICLEFILES))
 M_TRS_FILES += 	run_block_byte_tests_matlab.m.trs \
 		run_reporting_test_matlab.m.trs \
-		run_all_unitary_tests.m.trs \
+		run_all_unit_tests.m.trs \
 		histval_initval_file_unit_tests.m.trs \
 		test_aggregate_routine_1_2.m.trs \
 		test_aggregate_routine_1_2_3.m.trs \
@@ -1276,7 +1276,7 @@ O_TRS_FILES = $(patsubst %.mod, %.o.trs, $(MODFILES))
 O_TRS_FILES += $(patsubst %.mod, %.o.trs, $(PARTICLEFILES))
 O_TRS_FILES += 	run_block_byte_tests_octave.o.trs \
 		run_reporting_test_octave.o.trs \
-		run_all_unitary_tests.o.trs \
+		run_all_unit_tests.o.trs \
 		histval_initval_file_unit_tests.o.trs \
 		run_kronecker_tests.o.trs \
 		nonlinearsolvers.o.trs \
@@ -1306,7 +1306,7 @@ EXTRA_DIST = \
 	run_block_byte_tests_octave.m \
 	run_reporting_test_matlab.m \
 	run_reporting_test_octave.m \
-	run_all_unitary_tests.m \
+	run_all_unit_tests.m \
 	+matlab/+namespace/y_k.m \
 	reporting/AnnualTable.m \
 	reporting/CommResidTablePage.m \
diff --git a/tests/run_all_unitary_tests.m b/tests/run_all_unit_tests.m
similarity index 92%
rename from tests/run_all_unitary_tests.m
rename to tests/run_all_unit_tests.m
index c8d65b4fc07382b78d52e384f4c5a1a2ab499964..302653311f98bfed36f795188d5d8a67219fbf6f 100644
--- a/tests/run_all_unitary_tests.m
+++ b/tests/run_all_unit_tests.m
@@ -1,4 +1,4 @@
-% Copyright © 2013-2022 Dynare Team
+% Copyright © 2013-2023 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -52,7 +52,7 @@ counter = 0;
 
 for i = 1:length(mlist)
     f = [top_test_dir filesep mlist{i} ];
-    if is_unitary_test_available(f)
+    if is_unit_test_available(f)
         [check, info] = mtest(f);
         for j = 1:size(info, 1)
             counter = counter + 1;
@@ -65,9 +65,9 @@ end
 
 cd(getenv('TOP_TEST_DIR'));
 if isoctave
-    fid = fopen('run_all_unitary_tests.o.trs', 'w+');
+    fid = fopen('run_all_unit_tests.o.trs', 'w+');
 else
-    fid = fopen('run_all_unitary_tests.m.trs', 'w+');
+    fid = fopen('run_all_unit_tests.m.trs', 'w+');
 end
 if length(failedtests) > 0
   fprintf(fid,':test-result: FAIL\n');