diff --git a/tests/Makefile.am b/tests/Makefile.am
index b278925dc6a4d5884391593c3f09a70ec24d6b8c..543753d99757fb48215fb822c88aa1b3c151ce30 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -371,9 +371,6 @@ MODFILES = \
 	observation_trends_and_prefiltering/calib_smoother/Tr_prefilter_loglin_calib_smoother.mod \
 	observation_trends_and_prefiltering/calib_smoother/Tr_prefil_f_obs_loglin_cal_smoother.mod \
 	observation_trends_and_prefiltering/ML/Trend_no_prefilter_selected_var.mod \
-	bgp/solow-1/solow.mod \
-	bgp/nk-1/nk.mod \
-	bgp/ramsey-1/ramsey.mod \
 	dynare-command-options/ramst.mod \
 	particle/local_state_space_iteration_k_test.mod
 
diff --git a/tests/bgp/fs2000/fs2000.mod b/tests/bgp/fs2000/fs2000.mod
deleted file mode 100644
index e5ef2fc4994806565f29838bd906f637db6a9d84..0000000000000000000000000000000000000000
--- a/tests/bgp/fs2000/fs2000.mod
+++ /dev/null
@@ -1,72 +0,0 @@
-/*
- * This file is a modified version of 'fs2000.mod'.
- *
- * The difference is that, here, the equations are written in non-stationary form, and we test if 
- * we are able to identify the trends.
- *
- */
-
-/*
- * Copyright (C) 2019 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-var gM M;
-var gA A k c y;
-var P;             // follows M(-1)/A
-var W l d;         // follows M(-1)
-var R n;
-
-varexo e_a e_m;
-
-parameters alp bet gam mst rho psi del;
-
-alp = 0.33;
-bet = 0.99;
-gam = 0.003;
-mst = 1.011;
-rho = 0.7;
-psi = 0.787;
-del = 0.02;
-
-model;
-A = gA*A(-1);
-M = gM*M(-1);
-gA = exp(gam+e_a);
-log(gM) = (1-rho)*log(mst) + rho*log(gM(-1))+e_m;
-c+k = k(-1)^alp*(A*n)^(1-alp)+(1-del)*k(-1);
-P*c = M;
-P/(c(+1)*P(+1))=bet*P(+1)*(alp*k^(alp-1)*(A(+1)*n(+1))^(1-alp)+(1-del))/(c(+2)*P(+2));
-(psi/(1-psi))*(c*P/(1-n))=W;
-R = P*(1-alp)*k(-1)^alp*A^(1-alp)*n^(-alp)/W;
-W = l/n;
-M-M(-1)+d = l;
-1/(c*P)=bet*R/(c(+1)*P(+1));
-y = k(-1)^alp*(A*n)^(1-alp);
-end;
-
-verbatim;
-
-    bgp.write(M_);
-    options = optimoptions('fsolve','Display','iter','Algorithm','levenberg-marquardt','MaxFunctionEvaluations',1000000,'MaxIterations',100000,'SpecifyObjectiveGradient',true,'FunctionTolerance',1e-6,'StepTolerance',1e-6);
-    y = 1+(rand(M_.endo_nbr,1)-.5)*.5;
-    g = ones(M_.endo_nbr,1);% 1+(rand(M_.endo_nbr,1)-.5)*.1;
-    [y, fval, exitflag] = fsolve(@fs2000.bgpfun, [y;g], options);
-    y(1:M_.orig_endo_nbr)
-    y(M_.endo_nbr+(1:M_.orig_endo_nbr))
-
-end;
\ No newline at end of file
diff --git a/tests/bgp/nk-1/nk.mod b/tests/bgp/nk-1/nk.mod
deleted file mode 100644
index aa0c3721a9ab2bb250dfb14b31fe15d0813be09a..0000000000000000000000000000000000000000
--- a/tests/bgp/nk-1/nk.mod
+++ /dev/null
@@ -1,41 +0,0 @@
-var y i pi ;
-
-parameters a1 a2 a3 a4 a5;
-
-a1 = -.5;
-a2 =  .1;
-a3 =  .9;
-a4 = 1.5;
-a5 = 0.5;
-
-model;
-
-y = y(1)*(i/pi(1))^a1;
-
-pi = (y^a2)*(pi(1)^a3);
-
-i = (pi^a4)*(y^a5);
-
-end;
-
-verbatim;
-
-    bgp.write(M_);
-    if isoctave || matlab_ver_less_than('8.1')
-        options = optimset('Display', 'iter', 'MaxFunEvals', 1000000,'MaxIter',100000,'Jacobian','on','TolFun',1e-8,'TolX',1e-8);
-    else
-        options = optimoptions('fsolve','Display','iter','Algorithm','levenberg-marquardt','MaxFunctionEvaluations',1000000,'MaxIterations',100000,'SpecifyObjectiveGradient',true,'FunctionTolerance',1e-8,'StepTolerance',1e-8);
-    end
-    if isoctave
-        % Octave can't take a function handle of a function in a package
-        % See https://savannah.gnu.org/bugs/index.php?46659
-        fun = str2func('nk.bgpfun');
-    else
-        fun = @nk.bgpfun;
-    end
-    y = 1+(rand(3,1)-.5)*.5;
-    g = 1+(rand(3,1)-.5)*.1;
-    [y, fval, exitflag] = fsolve(fun, [y;g], options);
-    assert(max(abs(y-1))<1e-9);
-
-end;
diff --git a/tests/bgp/ramsey-1/ramsey.mod b/tests/bgp/ramsey-1/ramsey.mod
deleted file mode 100644
index 54e71e58378ff072fcfe890e6e451d2acd8c8d08..0000000000000000000000000000000000000000
--- a/tests/bgp/ramsey-1/ramsey.mod
+++ /dev/null
@@ -1,36 +0,0 @@
-var c k x;
-
-parameters alpha gam delta beta ;
-
-alpha=0.5;
-gam=0.5;
-delta=0.02;
-beta=0.05;
-
-model;
-x = x(-1)*1.02;
-c + k - x^(1-alpha)*k(-1)^alpha - (1-delta)*k(-1);
-c^(-gam) - (1+beta)^(-1)*(alpha*x(+1)^(1-alpha)*k^(alpha-1) + 1 - delta)*c(+1)^(-gam);
-end;
-
-verbatim;
-
-    bgp.write(M_);
-    if isoctave || matlab_ver_less_than('8.1')
-        options = optimset('Display', 'iter', 'MaxFunEvals', 1000000,'MaxIter',100000,'Jacobian','on','TolFun',1e-6,'TolX',1e-6);
-    else
-        options = optimoptions('fsolve','Display','iter','Algorithm','levenberg-marquardt','MaxFunctionEvaluations',1000000,'MaxIterations',100000,'SpecifyObjectiveGradient',true,'FunctionTolerance',1e-6,'StepTolerance',1e-6);
-    end
-    if isoctave
-        % Octave can't take a function handle of a function in a package
-        % See https://savannah.gnu.org/bugs/index.php?46659
-        fun = str2func('ramsey.bgpfun');
-    else
-        fun = @ramsey.bgpfun;
-    end
-    y = 1+(rand(M_.endo_nbr,1)-.5)*.5;
-    g = ones(M_.endo_nbr,1);% 1+(rand(M_.endo_nbr,1)-.5)*.1;
-    [y, fval, exitflag] = fsolve(fun, [y;g], options);
-    assert(max(abs(y(M_.endo_nbr+(1:M_.orig_endo_nbr))-1.02))<1e-6)
-
-end;
diff --git a/tests/bgp/solow-1/solow.mod b/tests/bgp/solow-1/solow.mod
deleted file mode 100644
index 7ef45993fc9699ce8f012359d3bc6b1ce4a77cd1..0000000000000000000000000000000000000000
--- a/tests/bgp/solow-1/solow.mod
+++ /dev/null
@@ -1,91 +0,0 @@
-// This mod file test the bgp.write function by characterizing the Balanced Growth Path of the Solow model. This is done 10000 times in a Monte Carlo loop
-// randomizing the initial guess of the nonlinear solver.
-
-var Efficiency
-    EfficiencyGrowth
-    Population
-    PopulationGrowth
-    Output
-    PhysicalCapitalStock ;
-
-varexo e_x
-       e_n ;
-
-parameters alpha
-	       delta
-	       s
-           rho_x
-           rho_n
-           EfficiencyGrowth_ss
-           PopulationGrowth_ss ;
-
-alpha = .33;
-delta = .02;
-s     = .20;
-rho_x = .90;
-rho_n = .95;
-EfficiencyGrowth_ss = 1.01;
-PopulationGrowth_ss = 1.01;
-
-model;
-    Efficiency = EfficiencyGrowth*Efficiency(-1);
-    EfficiencyGrowth/EfficiencyGrowth_ss = (EfficiencyGrowth(-1)/EfficiencyGrowth_ss)^(rho_x)*exp(e_x);
-    Population = PopulationGrowth*Population(-1);
-    PopulationGrowth/PopulationGrowth_ss = (PopulationGrowth(-1)/PopulationGrowth_ss)^(rho_n)*exp(e_n);
-    Output = PhysicalCapitalStock(-1)^alpha*(Efficiency*Population)^(1-alpha);
-    PhysicalCapitalStock = (1-delta)*PhysicalCapitalStock(-1) + s*Output;
-end;
-
-
-verbatim;
-
-    bgp.write(M_);
-    MC = 10000;
-    KY = NaN(MC,1);
-    GY = NaN(MC,1);
-    GK = NaN(MC,1);
-    EG = NaN(MC,1);
-    if isoctave || matlab_ver_less_than('8.1')
-        options = optimset('Display', 'off', 'MaxFunEvals', 1000000,'MaxIter',100000,'Jacobian','on','TolFun',1e-8,'TolX',1e-8);
-    else
-        options = optimoptions('fsolve','Display','off','Algorithm','levenberg-marquardt','MaxFunctionEvaluations',1000000,'MaxIterations',100000,'SpecifyObjectiveGradient',true,'FunctionTolerance',1e-8,'StepTolerance',1e-8);
-    end
-    if isoctave
-        % Octave can't take a function handle of a function in a package
-        % See https://savannah.gnu.org/bugs/index.php?46659
-        fun = str2func('solow.bgpfun');
-    else
-        fun = @solow.bgpfun;
-    end
-    for i=1:MC
-        y = 1+(rand(6,1)-.5)*.2;
-        g = ones(6,1);
-        [y, fval, exitflag] = fsolve(fun, [y;g], options);
-        if exitflag>0
-            KY(i) = y(6)/y(5);
-            GY(i) = y(11);
-            GK(i) = y(12);
-            EG(i) = y(2);
-        end
-    end
-    % Compute the physical capital stock over output ratio along the BGP as
-    % a function of the deep parameters...
-    theoretical_long_run_ratio = s*EfficiencyGrowth_ss*PopulationGrowth_ss/(EfficiencyGrowth_ss*PopulationGrowth_ss-1+delta);
-    % ... and compare to the average ratio in the Monte Carlo
-    mu = mean(KY(~isnan(KY)));
-    s2 = var(KY(~isnan(KY)));
-    nn = sum(isnan(KY))/MC;
-    disp(sprintf('Theoretical K/Y BGP ratio = %s.', num2str(theoretical_long_run_ratio)));
-    disp(sprintf('mean(K/Y) = %s.', num2str(mu)));
-    disp(sprintf('var(K/Y) = %s.', num2str(s2)));
-    disp(sprintf('Number of failures: %u (over %u problems).', sum(isnan(KY)), MC));
-    assert(abs(mu-theoretical_long_run_ratio)<1e-8);
-    assert(s2<1e-16);
-    assert(abs(mean(GY(~isnan(GY)))-EfficiencyGrowth_ss*PopulationGrowth_ss)<1e-8)
-    assert(var(GY(~isnan(GY)))<1e-16);
-    assert(abs(mean(GK(~isnan(GK)))-EfficiencyGrowth_ss*PopulationGrowth_ss)<1e-8)
-    assert(var(GK(~isnan(GK)))<1e-16);
-    assert(abs(mean(EG(~isnan(EG)))-EfficiencyGrowth_ss)<1e-8)
-    assert(var(EG(~isnan(EG)))<1e-16);
-
-end;