diff --git a/src/DataTree.cc b/src/DataTree.cc
index d96c965b5f54a37280b2916c9cc876e3ab7a083e..f2d1717367983f98d87e8bde8a244740eca06d85 100644
--- a/src/DataTree.cc
+++ b/src/DataTree.cc
@@ -916,9 +916,9 @@ void
 DataTree::writePowerDerivJulia(ostream &output) const
 {
   if (isBinaryOpUsed(BinaryOpcode::powerDeriv))
-    output << "nearbyint(x::<: Real) = (abs((x)-floor(x)) < abs((x)-ceil(x)) ? floor(x) : ceil(x))" << endl
+    output << "nearbyint(x::T) where T <: Real  = (abs((x)-floor(x)) < abs((x)-ceil(x)) ? floor(x) : ceil(x))" << endl
 	   << endl
-	   << "function get_power_deriv(x::<: Real, p::<: Real, k::Int64)" << endl
+	   << "function get_power_deriv(x::T, p::T, k::Int64) where T <: Real" << endl
 	   << "    if (abs(x) < 1e-12 && p > 0 && k > p && abs(p-nearbyint(p)) < 1e-12 )" << endl
 	   << "        return 0.0" << endl
 	   << "    else" << endl
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 98636796086d9d7584153790c53693e0ab842e67..52d9af3ba3be8052e421dad59c7754b0976d48bd 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -1400,7 +1400,9 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
   output << "function dynamicResidTT!(T::Vector{<: Real}," << endl
          << "                         y::Vector{<: Real}, x::Matrix{<: Real}, "
          << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int)" << endl
+         << "@inbounds begin" << endl
          << tt_output[0].str()
+	 << "end" << endl
          << "    return nothing" << endl
          << "end" << endl << endl;
 
@@ -1415,7 +1417,9 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
          << "    if T_flag" << endl
          << "        dynamicResidTT!(T, y, x, params, steady_state, it_)" << endl
          << "    end" << endl
+         << "@inbounds begin" << endl
          << d_output[0].str()
+	 << "end" << endl
          << "    return nothing" << endl
          << "end" << endl << endl;
 
@@ -1424,7 +1428,9 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
          << "                      y::Vector{<: Real}, x::Matrix{<: Real}, "
          << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int)" << endl
          << "    dynamicResidTT!(T, y, x, params, steady_state, it_)" << endl
+         << "@inbounds begin" << endl
          << tt_output[1].str()
+	 << "end" << endl
          << "    return nothing" << endl
          << "end" << endl << endl;
 
@@ -1441,7 +1447,9 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
          << "        dynamicG1TT!(T, y, x, params, steady_state, it_)" << endl
          << "    end" << endl
          << "    fill!(g1, 0.0)" << endl
+         << "@inbounds begin" << endl
          << d_output[1].str()
+	 << "end" << endl
          << "    return nothing" << endl
          << "end" << endl << endl;
 
@@ -1450,7 +1458,9 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
          << "                      y::Vector{<: Real}, x::Matrix{<: Real}, "
          << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int)" << endl
          << "    dynamicG1TT!(T, y, x, params, steady_state, it_)" << endl
+         << "@inbounds begin" << endl
          << tt_output[2].str()
+	 << "end" << endl
          << "    return nothing" << endl
          << "end" << endl << endl;
 
@@ -1467,7 +1477,9 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
          << "        dynamicG2TT!(T, y, x, params, steady_state, it_)" << endl
          << "    end" << endl
          << "    fill!(g2, 0.0)" << endl
+         << "@inbounds begin" << endl
          << d_output[2].str()
+	 << "end" << endl
          << "    return nothing" << endl
          << "end" << endl << endl;
 
@@ -1476,7 +1488,9 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
          << "                      y::Vector{<: Real}, x::Matrix{<: Real}, "
          << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int)" << endl
          << "    dynamicG2TT!(T, y, x, params, steady_state, it_)" << endl
+         << "@inbounds begin" << endl
          << tt_output[3].str()
+	 << "end" << endl
          << "    return nothing" << endl
          << "end" << endl << endl;
 
@@ -1494,7 +1508,9 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
          << "      dynamicG3TT!(T, y, x, params, steady_state, it_)" << endl
          << "    end" << endl
          << "    fill!(g3, 0.0)" << endl
+         << "@inbounds begin" << endl
          << d_output[3].str()
+	 << "end" << endl
          << "    return nothing" << endl
          << "end" << endl << endl;
 
@@ -4331,8 +4347,10 @@ DynamicModel::writeSetAuxiliaryVariables(const string &basename, bool julia) con
          << comment << endl
          << comment << " Warning : this file is generated automatically by Dynare" << endl
          << comment << "           from model file (.mod)" << endl << endl
+         << "@inbounds begin" << endl    
          << output_func_body.str()
-         << "end" << endl;
+    	 << "end" << endl
+     << "end" << endl;
   if (julia)
     output << "end" << endl;
 
@@ -5049,20 +5067,34 @@ DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) con
                      << "export params_derivs" << endl << endl
                      << "function params_derivs(y, x, paramssteady_state, it_, "
                      << "ss_param_deriv, ss_param_2nd_deriv)" << endl
+		     << "@inbounds begin" << endl
                      << tt_output.str()
+		     << "end" << endl
                      << "rp = zeros(" << equations.size() << ", "
                      << symbol_table.param_nbr() << ");" << endl
+		     << "@inbounds begin" << endl
                      << rp_output.str()
+		     << "end" << endl
                      << "gp = zeros(" << equations.size() << ", " << getJacobianColsNbr() << ", " << symbol_table.param_nbr() << ");" << endl
+		     << "@inbounds begin" << endl
                      << gp_output.str()
+		     << "end" << endl
                      << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
+		     << "@inbounds begin" << endl
                      << rpp_output.str()
+		     << "end" << endl
                      << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
+		     << "@inbounds begin" << endl
                      << gpp_output.str()
+		     << "end" << endl
                      << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
+		     << "@inbounds begin" << endl
                      << hp_output.str()
+		     << "end" << endl
                      << "g3p = zeros(" << params_derivatives.find({ 3, 1 })->second.size() << ",6);" << endl
+		     << "@inbounds begin" << endl
                      << g3p_output.str()
+		     << "end" << endl
                      << "(rp, gp, rpp, gpp, hp, g3p)" << endl
                      << "end" << endl
                      << "end" << endl;
diff --git a/src/ModelEquationBlock.cc b/src/ModelEquationBlock.cc
index cd47e8725d94c72eab33e9ba14b846b0a446d74b..9b0551f46b7677308ca8f0789f6fbc3da28230c7 100644
--- a/src/ModelEquationBlock.cc
+++ b/src/ModelEquationBlock.cc
@@ -197,7 +197,8 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool julia) const
            << "#" << endl
            << "export steady_state!" << endl << endl
            << "function steady_state!(ys_::Vector{<: Real}, exo_::Vector{<: Real}, "
-           << "params::Vector{<: Real})" << endl;
+           << "params::Vector{<: Real})" << endl
+	   << "@inbounds begin" << endl;
 
   for (const auto & [symb_ids, value] : def_table)
     {
@@ -225,7 +226,7 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool julia) const
 
   output << "end" << endl;
   if (julia)
-    output << "end" << endl;
+    output << "end" << endl << "end" << endl;
 
   writeToFileIfModified(output, julia ? basename + "SteadyState2.jl" : packageDir(basename) + "/steadystate.m");
 }
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index a28f9df0cfa74c3ab9a7074fd2fc138768911877..2fc143fbe55d46d9041c7ea2bf2e96094de4d967 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -1046,7 +1046,7 @@ ModelTree::writeModelLocalVariableTemporaryTerms(temporary_terms_t &temp_term_un
       value->writeExternalFunctionOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
 
       if (isJuliaOutput(output_type))
-        output << "    @inbounds const ";
+        output << "    const ";
 
       mlv->writeOutput(output, output_type, tto, tt_idxs, tef_terms);
       output << " = ";
@@ -1074,9 +1074,6 @@ ModelTree::writeTemporaryTerms(const temporary_terms_t &tt,
       if (dynamic_cast<AbstractExternalFunctionNode *>(it))
         it->writeExternalFunctionOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
 
-      if (isJuliaOutput(output_type))
-        output << "    @inbounds ";
-
       it->writeOutput(output, output_type, tt, tt_idxs, tef_terms);
       output << " = ";
       it->writeOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
@@ -1348,7 +1345,7 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type,
       if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
         if (isJuliaOutput(output_type))
           {
-            output << "    @inbounds residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
+            output << "    residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
                    << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
                    << RIGHT_ARRAY_SUBSCRIPT(output_type)
                    << " = (";
@@ -1372,8 +1369,6 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type,
           }
       else // The right hand side of the equation is empty ==> residual=lhs;
         {
-          if (isJuliaOutput(output_type))
-            output << "    @inbounds ";
           output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
                  << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
                  << RIGHT_ARRAY_SUBSCRIPT(output_type)
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index c10592a3ed020f927538ce66faa1af4082124cc4..e0456fbdd310399e0552c49b24951e856b54736f 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -545,8 +545,6 @@ ModelTree::writeModelFileHelper() const
         {
           auto [eq, var] = vectorToTuple<2>(indices);
 
-          if constexpr(isJuliaOutput(output_type))
-            d_output[1] << "    @inbounds ";
           d_output[1] << "g1" << LEFT_ARRAY_SUBSCRIPT(output_type);
           if constexpr(isMatlabOutput(output_type) || isJuliaOutput(output_type))
             d_output[1] << eq + 1 << "," << getJacobianCol(var) + 1;
@@ -589,7 +587,7 @@ ModelTree::writeModelFileHelper() const
 
             if constexpr(isJuliaOutput(output_type))
               {
-                d_output[i] << "    @inbounds " << "g" << i << "[" << eq + 1 << "," << col_idx + 1 << "] = ";
+                d_output[i] << "    g" << i << "[" << eq + 1 << "," << col_idx + 1 << "] = ";
                 d->writeOutput(d_output[i], output_type, temp_term_union, temporary_terms_idxs, tef_terms);
                 d_output[i] << endl;
               }
@@ -618,7 +616,7 @@ ModelTree::writeModelFileHelper() const
                 int col_idx_sym{getJacobianCol(vidx[2]) * getJacobianColsNbr() + getJacobianCol(vidx[1])};
 
                 if constexpr(isJuliaOutput(output_type))
-                  d_output[2] << "    @inbounds g2[" << eq + 1 << "," << col_idx_sym + 1 << "] = "
+                  d_output[2] << "    g2[" << eq + 1 << "," << col_idx_sym + 1 << "] = "
                               << "g2[" << eq + 1 << "," << col_idx + 1 << "]" << endl;
                 else
                   {
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index b33eae8e96f66846af7815d0ed3ed7ea7d82eb33..d769bcc4f35a95f44ecf2153e62a9b7700e09cf4 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -1343,7 +1343,9 @@ StaticModel::writeStaticJuliaFile(const string &basename) const
   output << "function staticResidTT!(T::Vector{<: Real}," << endl
          << "                        y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real})" << endl
          << "    @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size()  << endl
+         << "    @inbounds begin" << endl
          << tt_output[0].str()
+	 << "    end" << endl
          << "    return nothing" << endl
          << "end" << endl << endl;
 
@@ -1357,7 +1359,9 @@ StaticModel::writeStaticJuliaFile(const string &basename) const
          << "    if T0_flag" << endl
          << "        staticResidTT!(T, y, x, params)" << endl
          << "    end" << endl
+         << "    @inbounds begin" << endl
          << d_output[0].str()
+	 << "    end" << endl
          << "    if ~isreal(residual)" << endl
          << "        residual = real(residual)+imag(residual).^2;" << endl
          << "    end" << endl
@@ -1370,7 +1374,9 @@ StaticModel::writeStaticJuliaFile(const string &basename) const
          << "    if T0_flag" << endl
          << "        staticResidTT!(T, y, x, params)" << endl
          << "    end" << endl
+         << "    @inbounds begin" << endl
          << tt_output[1].str()
+	 << "    end" << endl
          << "    return nothing" << endl
          << "end" << endl << endl;
 
@@ -1387,7 +1393,9 @@ StaticModel::writeStaticJuliaFile(const string &basename) const
          << "        staticG1TT!(T, y, x, params, T0_flag)" << endl
          << "    end" << endl
          << "    fill!(g1, 0.0)" << endl
+         << "     @inbounds begin" << endl
          << d_output[1].str()
+	 << "    end" << endl
          << "    if ~isreal(g1)" << endl
          << "        g1 = real(g1)+2*imag(g1);" << endl
          << "    end" << endl
@@ -1400,7 +1408,9 @@ StaticModel::writeStaticJuliaFile(const string &basename) const
          << "    if T1_flag" << endl
          << "        staticG1TT!(T, y, x, params, TO_flag)" << endl
          << "    end" << endl
+         << "    @inbounds begin" << endl
          << tt_output[2].str()
+	 << "    end" << endl
          << "    return nothing" << endl
          << "end" << endl << endl;
 
@@ -1418,7 +1428,9 @@ StaticModel::writeStaticJuliaFile(const string &basename) const
          << "        staticG2TT!(T, y, x, params, T1_flag, T0_flag)" << endl
          << "    end" << endl
          << "    fill!(g2, 0.0)" << endl
+         << "    @inbounds begin" << endl
          << d_output[2].str()
+	 << "    end" << endl
          << "    return nothing" << endl
          << "end" << endl << endl;
 
@@ -1428,7 +1440,9 @@ StaticModel::writeStaticJuliaFile(const string &basename) const
          << "    if T2_flag" << endl
          << "        staticG2TT!(T, y, x, params, T1_flag, T0_flag)" << endl
          << "    end" << endl
+         << "    @inbounds begin" << endl
          << tt_output[3].str()
+	 << "    end" << endl
          << "    return nothing" << endl
          << "end" << endl << endl;
 
@@ -1446,7 +1460,9 @@ StaticModel::writeStaticJuliaFile(const string &basename) const
          << "        staticG3TT!(T, y, x, params, T2_flag, T1_flag, T0_flag)" << endl
          << "    end" << endl
          << "    fill!(g3, 0.0)" << endl
+         << "    @inbounds begin" << endl
          << d_output[3].str()
+	 << "    end" << endl
          << "    return nothing" << endl
          << "end" << endl << endl;
 
@@ -1845,7 +1861,9 @@ StaticModel::writeSetAuxiliaryVariables(const string &basename, bool julia) cons
          << comment << endl
          << comment << " Warning : this file is generated automatically by Dynare" << endl
          << comment << "           from model file (.mod)" << endl << endl
+	 << "@inbounds begin" << endl
          << output_func_body.str()
+         << "end" << endl
          << "end" << endl;
   if (julia)
     output << "end" << endl;
@@ -2171,19 +2189,31 @@ StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) cons
                      << "#" << endl
                      << "export params_derivs" << endl << endl
                      << "function params_derivs(y, x, params)" << endl
+		     << "@inbounds begin" << endl
                      << tt_output.str()
-                     << "rp = zeros(" << equations.size() << ", "
+		     << "end" << endl
+		     << "rp = zeros(" << equations.size() << ", "
                      << symbol_table.param_nbr() << ");" << endl
+		     << "@inbounds begin" << endl
                      << jacobian_output.str()
+		     << "end" << endl
                      << "gp = zeros(" << equations.size() << ", " << symbol_table.endo_nbr() << ", "
                      << symbol_table.param_nbr() << ");" << endl
+		     << "@inbounds begin" << endl
                      << hessian_output.str()
+		     << "end" << endl
                      << "rpp = zeros(" << params_derivatives.find({ 0, 2 })->second.size() << ",4);" << endl
+		     << "@inbounds begin" << endl
                      << hessian1_output.str()
+		     << "end" << endl
                      << "gpp = zeros(" << params_derivatives.find({ 1, 2 })->second.size() << ",5);" << endl
+		     << "@inbounds begin" << endl
                      << third_derivs_output.str()
+		     << "end" << endl
                      << "hp = zeros(" << params_derivatives.find({ 2, 1 })->second.size() << ",5);" << endl
+		     << "@inbounds begin" << endl
                      << third_derivs1_output.str()
+		     << "end" << endl
                      << "(rp, gp, rpp, gpp, hp)" << endl
                      << "end" << endl
                      << "end" << endl;