diff --git a/src/DataTree.cc b/src/DataTree.cc
index 363a1647..f2d17173 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::Float64) = (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::Float64, p::Float64, 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 9e444e3a..52d9af3b 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -1364,14 +1364,14 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
          << "  dynamicG2TT!    : Computes the dynamic model temporary terms for the Hessian" << endl
          << "  dynamicG3TT!    : Computes the dynamic model temporary terms for the third derivatives" << endl << endl
          << "## Function Arguments ##" << endl
-         << "  T            : Vector{Float64}(num_temp_terms), temporary terms" << endl
-         << "  y            : Vector{Float64}(num_dynamic_vars), endogenous variables in the order stored model_.lead_lag_incidence; see the manual" << endl
-         << "  x            : Matrix{Float64}(nperiods,model_.exo_nbr), exogenous variables (in declaration order) for all simulation periods" << endl
-         << "  params       : Vector{Float64}(model_.param_nbr), parameter values in declaration order" << endl
-         << "  steady_state : Vector{Float64}(model_endo_nbr)" << endl
+         << "  T            : Vector{<: Real}(num_temp_terms), temporary terms" << endl
+         << "  y            : Vector{<: Real}(num_dynamic_vars), endogenous variables in the order stored model_.lead_lag_incidence; see the manual" << endl
+         << "  x            : Matrix{<: Real}(nperiods,model_.exo_nbr), exogenous variables (in declaration order) for all simulation periods" << endl
+         << "  params       : Vector{<: Real}(model_.param_nbr), parameter values in declaration order" << endl
+         << "  steady_state : Vector{<: Real}(model_endo_nbr)" << endl
          << "  it_          : Int, time period for exogenous variables for which to evaluate the model" << endl
-         << "  residual     : Vector{Float64}(model_.eq_nbr), residuals of the dynamic model equations in order of declaration of the equations." << endl
-         << "  g1           : Matrix{Float64}(model_.eq_nbr, num_dynamic_vars), Jacobian matrix of the dynamic model equations" << endl
+         << "  residual     : Vector{<: Real}(model_.eq_nbr), residuals of the dynamic model equations in order of declaration of the equations." << endl
+         << "  g1           : Matrix{<: Real}(model_.eq_nbr, num_dynamic_vars), Jacobian matrix of the dynamic model equations" << endl
          << "                 The rows and columns respectively correspond to equations in order of declaration and variables in order" << endl
          << "                 stored in model_.lead_lag_incidence" << endl
          << "  g2           : spzeros(model_.eq_nbr, (num_dynamic_vars)^2) Hessian matrix of the dynamic model equations" << endl
@@ -1397,17 +1397,19 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
          << "tmp_nbr[4] = " << temporary_terms_derivatives[3].size() << "# Number of temporary terms for g3 (third order derivates)" << endl << endl;
 
   // dynamicResidTT!
-  output << "function dynamicResidTT!(T::Vector{Float64}," << endl
-         << "                         y::Vector{Float64}, x::Matrix{Float64}, "
-         << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+  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;
 
   // dynamic!
-  output << "function dynamicResid!(T::Vector{Float64}, residual::AbstractVector{Float64}," << endl
-         << "                       y::Vector{Float64}, x::Matrix{Float64}, "
-         << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl
+  output << "function dynamicResid!(T::Vector{<: Real}, residual::AbstractVector{<: Real}," << endl
+         << "                       y::Vector{<: Real}, x::Matrix{<: Real}, "
+         << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int, T_flag::Bool)" << endl
          << "    @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() << endl
          << "    @assert length(residual) == " << equations.size() << endl
          << "    @assert length(y)+size(x, 2) == " << getJacobianColsNbr() << endl
@@ -1415,23 +1417,27 @@ 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;
 
   // dynamicG1TT!
-  output << "function dynamicG1TT!(T::Vector{Float64}," << endl
-         << "                      y::Vector{Float64}, x::Matrix{Float64}, "
-         << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+  output << "function dynamicG1TT!(T::Vector{<: Real}," << endl
+         << "                      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;
 
   // dynamicG1!
-  output << "function dynamicG1!(T::Vector{Float64}, g1::Matrix{Float64}," << endl
-         << "                    y::Vector{Float64}, x::Matrix{Float64}, "
-         << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl
+  output << "function dynamicG1!(T::Vector{<: Real}, g1::Matrix{<: Real}," << endl
+         << "                    y::Vector{<: Real}, x::Matrix{<: Real}, "
+         << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int, T_flag::Bool)" << endl
          << "    @assert length(T) >= "
          << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() << endl
          << "    @assert size(g1) == (" << equations.size() << ", " << getJacobianColsNbr() << ")" << endl
@@ -1441,24 +1447,28 @@ 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;
 
   // dynamicG2TT!
-  output << "function dynamicG2TT!(T::Vector{Float64}," << endl
-         << "                      y::Vector{Float64}, x::Matrix{Float64}, "
-         << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+  output << "function dynamicG2TT!(T::Vector{<: Real}," << endl
+         << "                      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;
 
   // dynamicG2!
   int hessianColsNbr{getJacobianColsNbr() * getJacobianColsNbr()};
-  output << "function dynamicG2!(T::Vector{Float64}, g2::Matrix{Float64}," << endl
-         << "                    y::Vector{Float64}, x::Matrix{Float64}, "
-         << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl
+  output << "function dynamicG2!(T::Vector{<: Real}, g2::Matrix{<: Real}," << endl
+         << "                    y::Vector{<: Real}, x::Matrix{<: Real}, "
+         << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int, T_flag::Bool)" << endl
          << "    @assert length(T) >= " << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() << endl
          << "    @assert size(g2) == (" << equations.size() << ", " << hessianColsNbr << ")" << endl
          << "    @assert length(y)+size(x, 2) == " << getJacobianColsNbr() << endl
@@ -1467,24 +1477,28 @@ 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;
 
   // dynamicG3TT!
-  output << "function dynamicG3TT!(T::Vector{Float64}," << endl
-         << "                      y::Vector{Float64}, x::Matrix{Float64}, "
-         << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+  output << "function dynamicG3TT!(T::Vector{<: Real}," << endl
+         << "                      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;
 
   // dynamicG3!
   int ncols{hessianColsNbr * getJacobianColsNbr()};
-  output << "function dynamicG3!(T::Vector{Float64}, g3::Matrix{Float64}," << endl
-         << "                    y::Vector{Float64}, x::Matrix{Float64}, "
-         << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int, T_flag::Bool)" << endl
+  output << "function dynamicG3!(T::Vector{<: Real}, g3::Matrix{<: Real}," << endl
+         << "                    y::Vector{<: Real}, x::Matrix{<: Real}, "
+         << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int, T_flag::Bool)" << endl
          << "    @assert length(T) >= "
          << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size() << endl
          << "    @assert size(g3) == (" << equations.size() << ", " << ncols << ")" << endl
@@ -1494,38 +1508,40 @@ 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;
 
   // dynamic!
-  output << "function dynamic!(T::Vector{Float64}, residual::AbstractVector{Float64}," << endl
-         << "                  y::Vector{Float64}, x::Matrix{Float64}, "
-         << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+  output << "function dynamic!(T::Vector{<: Real}, residual::AbstractVector{<: Real}," << endl
+         << "                  y::Vector{<: Real}, x::Matrix{<: Real}, "
+         << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int)" << endl
          << "    dynamicResid!(T, residual, y, x, params, steady_state, it_, true)" << endl
          << "    return nothing" << endl
          << "end" << endl
          << endl
-         << "function dynamic!(T::Vector{Float64}, residual::AbstractVector{Float64}, g1::Matrix{Float64}," << endl
-         << "                  y::Vector{Float64}, x::Matrix{Float64}, "
-         << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+         << "function dynamic!(T::Vector{<: Real}, residual::AbstractVector{<: Real}, g1::Matrix{<: Real}," << endl
+         << "                  y::Vector{<: Real}, x::Matrix{<: Real}, "
+         << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int)" << endl
          << "    dynamicG1!(T, g1, y, x, params, steady_state, it_, true)" << endl
          << "    dynamicResid!(T, residual, y, x, params, steady_state, it_, false)" << endl
          << "    return nothing" << endl
          << "end" << endl
          << endl
-         << "function dynamic!(T::Vector{Float64}, residual::AbstractVector{Float64}, g1::Matrix{Float64}, g2::Matrix{Float64}," << endl
-         << "                  y::Vector{Float64}, x::Matrix{Float64}, "
-         << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+         << "function dynamic!(T::Vector{<: Real}, residual::AbstractVector{<: Real}, g1::Matrix{<: Real}, g2::Matrix{<: Real}," << endl
+         << "                  y::Vector{<: Real}, x::Matrix{<: Real}, "
+         << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int)" << endl
          << "    dynamicG2!(T, g2, y, x, params, steady_state, it_, true)" << endl
          << "    dynamicG1!(T, g1, y, x, params, steady_state, it_, false)" << endl
          << "    dynamicResid!(T, residual, y, x, params, steady_state, it_, false)" << endl
          << "    return nothing" << endl
          << "end" << endl
          << endl
-         << "function dynamic!(T::Vector{Float64}, residual::AbstractVector{Float64}, g1::Matrix{Float64}, g2::Matrix{Float64}, g3::Matrix{Float64}," << endl
-         << "                  y::Vector{Float64}, x::Matrix{Float64}, "
-         << "params::Vector{Float64}, steady_state::Vector{Float64}, it_::Int)" << endl
+         << "function dynamic!(T::Vector{<: Real}, residual::AbstractVector{<: Real}, g1::Matrix{<: Real}, g2::Matrix{<: Real}, g3::Matrix{<: Real}," << endl
+         << "                  y::Vector{<: Real}, x::Matrix{<: Real}, "
+         << "params::Vector{<: Real}, steady_state::Vector{<: Real}, it_::Int)" << endl
          << "    dynamicG3!(T, g3, y, x, params, steady_state, it_, true)" << endl
          << "    dynamicG2!(T, g2, y, x, params, steady_state, it_, false)" << endl
          << "    dynamicG1!(T, g1, y, x, params, steady_state, it_, false)" << 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 02290549..9b0551f4 100644
--- a/src/ModelEquationBlock.cc
+++ b/src/ModelEquationBlock.cc
@@ -196,8 +196,9 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool julia) const
            << "#     from " << basename << ".mod" << endl
            << "#" << endl
            << "export steady_state!" << endl << endl
-           << "function steady_state!(ys_::Vector{Float64}, exo_::Vector{Float64}, "
-           << "params::Vector{Float64})" << endl;
+           << "function steady_state!(ys_::Vector{<: Real}, exo_::Vector{<: Real}, "
+           << "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 a28f9df0..2fc143fb 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 c10592a3..e0456fbd 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 5777967f..d769bcc4 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -1308,14 +1308,14 @@ StaticModel::writeStaticJuliaFile(const string &basename) const
          << "  staticG2TT!    : Computes the static model temporary terms for the Hessian" << endl
          << "  staticG3TT!    : Computes the static model temporary terms for the third derivatives" << endl << endl
          << "## Function Arguments ##" << endl
-         << "  T        : Vector{Float64}(num_temp_terms) temporary terms" << endl
-         << "  y        : Vector{Float64}(model_.endo_nbr) endogenous variables in declaration order" << endl
-         << "  x        : Vector{Float64}(model_.exo_nbr) exogenous variables in declaration order" << endl
-         << "  params   : Vector{Float64}(model_.param) parameter values in declaration order" << endl
-         << "  residual : Vector{Float64}(model_.eq_nbr) residuals of the static model equations" << endl
+         << "  T        : Vector{<: Real}(num_temp_terms) temporary terms" << endl
+         << "  y        : Vector{<: Real}(model_.endo_nbr) endogenous variables in declaration order" << endl
+         << "  x        : Vector{<: Real}(model_.exo_nbr) exogenous variables in declaration order" << endl
+         << "  params   : Vector{<: Real}(model_.param) parameter values in declaration order" << endl
+         << "  residual : Vector{<: Real}(model_.eq_nbr) residuals of the static model equations" << endl
          << "             in order of declaration of the equations. Dynare may prepend auxiliary equations," << endl
          << "             see model.aux_vars" << endl
-         << "  g1       : Matrix{Float64}(model.eq_nbr,model_.endo_nbr) Jacobian matrix of the static model equations" << endl
+         << "  g1       : Matrix{<: Real}(model.eq_nbr,model_.endo_nbr) Jacobian matrix of the static model equations" << endl
          << "             The columns and rows respectively correspond to the variables in declaration order and the" << endl
          << "             equations in order of declaration" << endl
          << "  g2       : spzeros(model.eq_nbr, model_.endo^2) Hessian matrix of the static model equations" << endl
@@ -1340,16 +1340,18 @@ StaticModel::writeStaticJuliaFile(const string &basename) const
          << "tmp_nbr[4] = " << temporary_terms_derivatives[3].size() << "# Number of temporary terms for g3 (third order derivates)" << endl << endl;
 
   // staticResidTT!
-  output << "function staticResidTT!(T::Vector{Float64}," << endl
-         << "                        y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
+  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;
 
   // static!
-  output << "function staticResid!(T::Vector{Float64}, residual::Vector{Float64}," << endl
-         << "                      y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T0_flag::Bool)" << endl
+  output << "function staticResid!(T::Vector{<: Real}, residual::Vector{<: Real}," << endl
+         << "                      y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}, T0_flag::Bool)" << endl
          << "    @assert length(y) == " << symbol_table.endo_nbr() << endl
          << "    @assert length(x) == " << symbol_table.exo_nbr() << endl
          << "    @assert length(params) == " << symbol_table.param_nbr() << 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
@@ -1365,18 +1369,20 @@ StaticModel::writeStaticJuliaFile(const string &basename) const
          << "end" << endl << endl;
 
   // staticG1TT!
-  output << "function staticG1TT!(T::Vector{Float64}," << endl
-         << "                     y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T0_flag::Bool)" << endl
+  output << "function staticG1TT!(T::Vector{<: Real}," << endl
+         << "                     y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}, T0_flag::Bool)" << endl
          << "    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;
 
   // staticG1!
-  output << "function staticG1!(T::Vector{Float64}, g1::Matrix{Float64}," << endl
-         << "                   y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T1_flag::Bool, T0_flag::Bool)" << endl
+  output << "function staticG1!(T::Vector{<: Real}, g1::Matrix{<: Real}," << endl
+         << "                   y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}, T1_flag::Bool, T0_flag::Bool)" << endl
          << "    @assert length(T) >= "
          << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() << endl
          << "    @assert size(g1) == (" << equations.size() << ", " << symbol_table.endo_nbr() << ")" << 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
@@ -1395,19 +1403,21 @@ StaticModel::writeStaticJuliaFile(const string &basename) const
          << "end" << endl << endl;
 
   // staticG2TT!
-  output << "function staticG2TT!(T::Vector{Float64}," << endl
-         << "                     y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T1_flag::Bool, T0_flag::Bool)" << endl
+  output << "function staticG2TT!(T::Vector{<: Real}," << endl
+         << "                     y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}, T1_flag::Bool, T0_flag::Bool)" << endl
          << "    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;
 
   // staticG2!
   int hessianColsNbr{symbol_table.endo_nbr() * symbol_table.endo_nbr()};
-  output << "function staticG2!(T::Vector{Float64}, g2::Matrix{Float64}," << endl
-         << "                   y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
+  output << "function staticG2!(T::Vector{<: Real}, g2::Matrix{<: Real}," << endl
+         << "                   y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
          << "    @assert length(T) >= "
          << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() << endl
          << "    @assert size(g2) == (" << equations.size() << ", " << hessianColsNbr << ")" << endl
@@ -1418,24 +1428,28 @@ 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;
 
   // staticG3TT!
-  output << "function staticG3TT!(T::Vector{Float64}," << endl
-         << "                     y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
+  output << "function staticG3TT!(T::Vector{<: Real}," << endl
+         << "                     y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
          << "    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;
 
   // staticG3!
   int ncols{hessianColsNbr * symbol_table.endo_nbr()};
-  output << "function staticG3!(T::Vector{Float64}, g3::Matrix{Float64}," << endl
-         << "                   y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64}, T3_flag::Bool, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
+  output << "function staticG3!(T::Vector{<: Real}, g3::Matrix{<: Real}," << endl
+         << "                   y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real}, T3_flag::Bool, T2_flag::Bool, T1_flag::Bool, T0_flag::Bool)" << endl
          << "    @assert length(T) >= "
          << temporary_terms_mlv.size() + temporary_terms_derivatives[0].size() + temporary_terms_derivatives[1].size() + temporary_terms_derivatives[2].size() + temporary_terms_derivatives[3].size() << endl
          << "    @assert size(g3) == (" << equations.size() << ", " << ncols << ")" << endl
@@ -1446,32 +1460,34 @@ 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;
 
   // static!
-  output << "function static!(T::Vector{Float64}, residual::Vector{Float64}," << endl
-         << "                  y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
+  output << "function static!(T::Vector{<: Real}, residual::Vector{<: Real}," << endl
+         << "                  y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real})" << endl
          << "    staticResid!(T, residual, y, x, params, true)" << endl
          << "    return nothing" << endl
          << "end" << endl
          << endl
-         << "function static!(T::Vector{Float64}, residual::Vector{Float64}, g1::Matrix{Float64}," << endl
-         << "                 y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
+         << "function static!(T::Vector{<: Real}, residual::Vector{<: Real}, g1::Matrix{<: Real}," << endl
+         << "                 y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real})" << endl
          << "    staticG1!(T, g1, y, x, params, true, true)" << endl
          << "    staticResid!(T, residual, y, x, params, false)" << endl
          << "    return nothing" << endl
          << "end" << endl
          << endl
-         << "function static!(T::Vector{Float64}, g1::Matrix{Float64}," << endl
-         << "                 y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
+         << "function static!(T::Vector{<: Real}, g1::Matrix{<: Real}," << endl
+         << "                 y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real})" << endl
          << "    staticG1!(T, g1, y, x, params, true, false)" << endl
          << "    return nothing" << endl
          << "end" << endl
          << endl
-         << "function static!(T::Vector{Float64}, residual::Vector{Float64}, g1::Matrix{Float64}, g2::Matrix{Float64}," << endl
-         << "                 y::Vector{Float64}, x::Vector{Float64}, params::Vector{Float64})" << endl
+         << "function static!(T::Vector{<: Real}, residual::Vector{<: Real}, g1::Matrix{<: Real}, g2::Matrix{<: Real}," << endl
+         << "                 y::Vector{<: Real}, x::Vector{<: Real}, params::Vector{<: Real})" << endl
          << "    staticG2!(T, g2, y, x, params, true, true, true)" << endl
          << "    staticG1!(T, g1, y, x, params, false, false)" << endl
          << "    staticResid!(T, residual, y, x, params, false)" << 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;