diff --git a/.gitignore b/.gitignore
index 330c3f41f927d66c0275e9aa7d63aa828969accb..dfe1d234e7d229276e66f1c6e0cf79af605d53bc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -119,6 +119,9 @@ mex/build/matlab/run_m2html.m
 /matlab/preprocessor*
 /matlab/dynare_version.m
 
+# JULIA dir
+/julia/preprocessor*
+
 # DLL rules
 *.mex
 *.dll
@@ -200,3 +203,6 @@ mex/build/matlab/run_m2html.m
 # Reporting
 *synctex.gz
 tests/reporting/tmpRepDir
+
+# Julia Tests
+/tests/**/*.jl
diff --git a/julia/Dynare.jl b/julia/Dynare.jl
new file mode 100644
index 0000000000000000000000000000000000000000..240396231b9313c831167f4ab193db0ee0f3866f
--- /dev/null
+++ b/julia/Dynare.jl
@@ -0,0 +1,44 @@
+module Dynare
+##
+ # Copyright (C) 2015 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/>.
+##
+
+export dynare, @dynare
+
+function dynare(modfile)
+    # Add cd to path if not already there
+    if isempty(findin([pwd()], LOAD_PATH))
+        unshift!(LOAD_PATH, pwd())
+    end
+
+    # Process modfile
+    println(string("Using ", WORD_SIZE, "-bit preprocessor"))
+    preprocessor = string(dirname(@__FILE__()), "/preprocessor", WORD_SIZE, "/dynare_m")
+    run(`$preprocessor $modfile language=julia output=dynamic`)
+
+    # Load module created by preprocessor
+    basename = split(modfile, ".mod"; keep=false)
+    require(basename[1])
+end
+
+
+macro dynare(modelname)
+    :(dynare($modelname))
+end
+
+end
diff --git a/julia/DynareModel.jl b/julia/DynareModel.jl
new file mode 100644
index 0000000000000000000000000000000000000000..6b37b8d71927be2a6570a12af6d4adf4991e53c1
--- /dev/null
+++ b/julia/DynareModel.jl
@@ -0,0 +1,176 @@
+module DynareModel
+##
+ # Copyright (C) 2015 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/>.
+##
+
+
+export Endo, Exo, ExoDet, Param, dynare_model
+
+abstract Atom
+
+immutable Endo <: Atom
+    name::UTF8String
+    tex_name::UTF8String
+    long_name::UTF8String
+end
+
+immutable Exo <: Atom
+    name::UTF8String
+    tex_name::UTF8String
+    long_name::UTF8String
+end
+
+immutable ExoDet <: Atom
+    name::UTF8String
+    tex_name::UTF8String
+    long_name::UTF8String
+end
+
+immutable Param <: Atom
+    name::UTF8String
+    tex_name::UTF8String
+    long_name::UTF8String
+end
+
+immutable AuxVars
+    endo_index::Int
+    var_type::Int
+    orig_index::Int
+    orig_lead_lag::Int
+    eq_nbr::Int
+    orig_expr::UTF8String
+end
+
+immutable PredVars
+    index::Int
+end
+
+immutable ObsVars
+    index::Int
+end
+
+immutable DetShocks
+    exo_det::Int
+    exo_id::Int
+    multiplicative::Bool
+    periods::Vector{Int}
+    value::Float64
+end
+
+immutable EquationTag
+    eq_nbr::Int
+    name::UTF8String
+    value::UTF8String
+end
+
+type Model
+    fname::ASCIIString
+    dname::ASCIIString
+    dynare_version::ASCIIString
+    endo::Vector{Endo}
+    exo::Vector{Exo}
+    exo_det::Vector{ExoDet}
+    param::Vector{Param}
+    aux_vars::Vector{AuxVars}
+    pred_vars::Vector{Int}
+    obs_vars::Vector{Int}
+    orig_endo_nbr::Int
+    orig_eq_nbr::Int
+    eq_nbr::Int
+    ramsey_eq_nbr::Int
+    det_shocks::Vector{DetShocks}
+    nstatic::Int
+    nfwrd::Int
+    npred::Int
+    nboth::Int
+    nsfwrd::Int
+    nspred::Int
+    ndynamic::Int
+    maximum_lag::Int
+    maximum_lead::Int
+    maximum_endo_lag::Int
+    maximum_endo_lead::Int
+    maximum_exo_lag::Int
+    maximum_exo_lead::Int
+    lead_lag_incidence::Matrix{Int}
+    nnzderivatives::Vector{Int}
+    static_and_dynamic_models_differ::Bool
+    equation_tags::Vector{UTF8String}
+    exo_names_orig_ord::Vector{Int}
+    sigma_e::Matrix{Float64}
+    correlation_matrix::Matrix{Float64}
+    h::Matrix{Float64}
+    correlation_matrix_me::Matrix{Float64}
+    sigma_e_is_diagonal::Bool
+    params::Matrix{Float64}
+    static::Function
+    static_params_derivs::Function
+    dynamic::Function
+    dynamic_params_derivs::Function
+    steady_state::Function
+end
+
+function dynare_model()
+    return Model("",                    # fname
+                 "",                    # dname
+                 "",                    # dynare_version
+                 Array(Endo,0),         # endo
+                 Array(Exo,0),          # exo
+                 Array(ExoDet,0),       # exo_det
+                 Array(Param,0),        # param
+                 Array(AuxVars,0),      # aux_vars
+                 Array(Int,0),          # pred_vars
+                 Array(Int,0),          # obs_vars
+                 0,                     # orig_endo_nbr
+                 0,                     # orig_eq_nbr
+                 0,                     # eq_nbr
+                 0,                     # ramsey_eq_nbr
+                 Array(DetShocks,0),    # det_shocks
+                 0,                     # nstatic
+                 0,                     # nfwrd
+                 0,                     # npred
+                 0,                     # nboth
+                 0,                     # nsfwrd
+                 0,                     # nspred
+                 0,                     # ndynamic
+                 0,                     # maximum_lag
+                 0,                     # maximum_lead
+                 0,                     # maximum_endo_lag
+                 0,                     # maximum_endo_lead
+                 0,                     # maximum_exo_lag
+                 0,                     # maximum_exo_lead
+                 Array(Int, 3, 0),      # lead_lag_incidence
+                 zeros(Int, 3),         # nnzderivatives
+                 false,                 # static_and_dynamic_models_differ
+                 Array(ASCIIString,0),  # equation_tags
+                 Array(Int64,1),        # exo_names_orig_ord
+                 Array(Float64, 0, 0),  # sigma_e (Cov matrix of the structural innovations)
+                 Array(Float64, 0, 0),  # correlation_matrix (Corr matrix of the structural innovations)
+                 Array(Float64, 0, 0),  # h (Cov matrix of the measurement errors)
+                 Array(Float64, 0, 0),  # correlation_matrix_me (Cov matrix of the measurement errors)
+                 true,                  # sigma_e_is_diagonal
+                 Array(Float64, 0, 0),  # params
+                 function()end,         # static
+                 function()end,         # static_params_derivs
+                 function()end,         # dynamic
+                 function()end,         # dynamic_params_derivs
+                 function()end          # steady_state
+                )
+end
+
+end
diff --git a/julia/DynareOptions.jl b/julia/DynareOptions.jl
new file mode 100644
index 0000000000000000000000000000000000000000..bdcca5c1621b3334526537040265705f88100d25
--- /dev/null
+++ b/julia/DynareOptions.jl
@@ -0,0 +1,35 @@
+module DynareOptions
+##
+ # Copyright (C) 2015 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/>.
+##
+
+
+export dynare_options
+
+type Options
+    dynare_version::ASCIIString
+    linear::Bool
+end
+
+function dynare_options()
+    return Options("",          # dynare_version
+                   false        # linear
+                  )
+end
+
+end
diff --git a/julia/DynareOutput.jl b/julia/DynareOutput.jl
new file mode 100644
index 0000000000000000000000000000000000000000..f8aea529d75139ab8147a88c7838c618e99dbcde
--- /dev/null
+++ b/julia/DynareOutput.jl
@@ -0,0 +1,37 @@
+module DynareOutput
+##
+ # Copyright (C) 2015 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/>.
+##
+
+
+export dynare_output
+
+type Output
+    dynare_version::ASCIIString
+    steady_state::Matrix{Float64}
+    exo_steady_state::Matrix{Float64}
+end
+
+function dynare_output()
+    return Output("",                   # dynare_version
+                  Array(Float64, 0, 0), # steady_state
+                  Array(Float64, 0, 0)  # exo_steady_state
+                 )
+end
+
+end
diff --git a/julia/Utils.jl b/julia/Utils.jl
new file mode 100644
index 0000000000000000000000000000000000000000..6fd7eff5316a3f97410331793983eed1a5404f6b
--- /dev/null
+++ b/julia/Utils.jl
@@ -0,0 +1,36 @@
+module Utils
+##
+ # Copyright (C) 2015 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/>.
+##
+
+export get_power_deriv
+
+function get_power_deriv(x::Float64, p::Real, k::Int)
+    if abs(x)<1e-12 && p>0 && k>p && typeof(p)==Int
+        dxp = .0
+    else
+        dxp = x^(p-k)
+        for i = 0:k-1
+            dxp *= p
+            p -= 1
+        end
+    end
+    return dxp
+end
+
+end
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index eaeebf2bc27bde15f83b0a3b6af732edc0dd86ce..72a26f3506d0eda22b790d50ad527b2d2f26245a 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -1206,7 +1206,7 @@ PlannerObjectiveStatement::computingPass()
 void
 PlannerObjectiveStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
 {
-  model_tree->writeStaticFile(basename + "_objective", false, false, false);
+  model_tree->writeStaticFile(basename + "_objective", false, false, false, false);
 }
 
 BVARDensityStatement::BVARDensityStatement(int maxnlags_arg, const OptionsList &options_list_arg) :
diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc
index 97616691ae29afe443a1aa5e3304065a6f36cc97..0481ec73ab54a18f70f0f6d298ceffcb7a0651cb 100644
--- a/preprocessor/DynamicModel.cc
+++ b/preprocessor/DynamicModel.cc
@@ -1557,11 +1557,36 @@ DynamicModel::writeDynamicMFile(const string &dynamic_basename) const
                     << "% Warning : this file is generated automatically by Dynare" << endl
                     << "%           from model file (.mod)" << endl << endl;
 
-  writeDynamicModel(mDynamicModelFile, false);
+  writeDynamicModel(mDynamicModelFile, false, false);
   mDynamicModelFile << "end" << endl; // Close *_dynamic function
   mDynamicModelFile.close();
 }
 
+void
+DynamicModel::writeDynamicJuliaFile(const string &basename) const
+{
+  string filename = basename + "Dynamic.jl";
+
+  ofstream output;
+  output.open(filename.c_str(), ios::out | ios::binary);
+  if (!output.is_open())
+    {
+      cerr << "Error: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  output << "module " << basename << "Dynamic" << endl
+         << "#" << endl
+         << "# NB: this file was automatically generated by Dynare" << endl
+         << "#     from " << basename << ".mod" << endl
+         << "#" << endl
+         << "using Utils" << endl << endl
+         << "export dynamic!" << endl << endl;
+  writeDynamicModel(output, false, true);
+  output << "end" << endl;
+  output.close();
+}
+
 void
 DynamicModel::writeDynamicCFile(const string &dynamic_basename, const int order) const
 {
@@ -1596,7 +1621,7 @@ DynamicModel::writeDynamicCFile(const string &dynamic_basename, const int order)
   writePowerDerivCHeader(mDynamicModelFile);
 
   // Writing the function body
-  writeDynamicModel(mDynamicModelFile, true);
+  writeDynamicModel(mDynamicModelFile, true, false);
 
   writePowerDeriv(mDynamicModelFile, true);
   mDynamicModelFile.close();
@@ -2085,21 +2110,23 @@ DynamicModel::writeSparseDynamicMFile(const string &dynamic_basename, const stri
 }
 
 void
-DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
+DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia) const
 {
-  ostringstream model_output;    // Used for storing model equations
+  ostringstream model_output;    // Used for storing model
+  ostringstream model_eq_output; // Used for storing model equations
   ostringstream jacobian_output; // Used for storing jacobian equations
   ostringstream hessian_output;  // Used for storing Hessian equations
   ostringstream third_derivatives_output;
 
-  ExprNodeOutputType output_type = (use_dll ? oCDynamicModel : oMatlabDynamicModel);
+  ExprNodeOutputType output_type = (use_dll ? oCDynamicModel :
+                                    julia ? oJuliaDynamicModel : oMatlabDynamicModel);
 
   deriv_node_temp_terms_t tef_terms;
   writeModelLocalVariables(model_output, output_type, tef_terms);
 
   writeTemporaryTerms(temporary_terms, model_output, output_type, tef_terms);
 
-  writeModelEquations(model_output, output_type);
+  writeModelEquations(model_eq_output, output_type);
 
   int nrows = equations.size();
   int hessianColsNbr = dynJacobianColsNbr * dynJacobianColsNbr;
@@ -2134,35 +2161,50 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
       int col_nb = id1 * dynJacobianColsNbr + id2;
       int col_nb_sym = id2 * dynJacobianColsNbr + id1;
 
-      sparseHelper(2, hessian_output, k, 0, output_type);
-      hessian_output << "=" << eq + 1 << ";" << endl;
-
-      sparseHelper(2, hessian_output, k, 1, output_type);
-      hessian_output << "=" << col_nb + 1 << ";" << endl;
-
-      sparseHelper(2, hessian_output, k, 2, output_type);
-      hessian_output << "=";
-      d2->writeOutput(hessian_output, output_type, temporary_terms, tef_terms);
-      hessian_output << ";" << endl;
-
-      k++;
-
-      // Treating symetric elements
-      if (id1 != id2)
+      ostringstream for_sym;
+      if (output_type == oJuliaDynamicModel)
+        {
+          for_sym << "g2[" << eq + 1 << "," << col_nb + 1 << "]";
+          hessian_output << "  @inbounds " << for_sym.str() << " = ";
+          d2->writeOutput(hessian_output, output_type, temporary_terms, tef_terms);
+          hessian_output << endl;
+        }
+      else
         {
           sparseHelper(2, hessian_output, k, 0, output_type);
           hessian_output << "=" << eq + 1 << ";" << endl;
 
           sparseHelper(2, hessian_output, k, 1, output_type);
-          hessian_output << "=" << col_nb_sym + 1 << ";" << endl;
+          hessian_output << "=" << col_nb + 1 << ";" << endl;
 
           sparseHelper(2, hessian_output, k, 2, output_type);
           hessian_output << "=";
-          sparseHelper(2, hessian_output, k-1, 2, output_type);
+          d2->writeOutput(hessian_output, output_type, temporary_terms, tef_terms);
           hessian_output << ";" << endl;
 
           k++;
         }
+
+      // Treating symetric elements
+      if (id1 != id2)
+        if (output_type == oJuliaDynamicModel)
+          hessian_output << "  @inbounds g2[" << eq + 1 << "," << col_nb_sym + 1 << "] = "
+                         << for_sym.str() << endl;
+        else
+          {
+            sparseHelper(2, hessian_output, k, 0, output_type);
+            hessian_output << "=" << eq + 1 << ";" << endl;
+
+            sparseHelper(2, hessian_output, k, 1, output_type);
+            hessian_output << "=" << col_nb_sym + 1 << ";" << endl;
+
+            sparseHelper(2, hessian_output, k, 2, output_type);
+            hessian_output << "=";
+            sparseHelper(2, hessian_output, k-1, 2, output_type);
+            hessian_output << ";" << endl;
+
+            k++;
+          }
     }
 
   // Writing third derivatives
@@ -2183,18 +2225,30 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
       // Reference column number for the g3 matrix
       int ref_col = id1 * hessianColsNbr + id2 * dynJacobianColsNbr + id3;
 
-      sparseHelper(3, third_derivatives_output, k, 0, output_type);
-      third_derivatives_output << "=" << eq + 1 << ";" << endl;
+      ostringstream for_sym;
+      if (output_type == oJuliaDynamicModel)
+        {
+          for_sym << "g3[" << eq + 1 << "," << ref_col + 1 << "]";
+          third_derivatives_output << "  @inbounds " << for_sym.str() << " = ";
+          d3->writeOutput(third_derivatives_output, output_type, temporary_terms, tef_terms);
+          third_derivatives_output << endl;
+        }
+      else
+        {
+          sparseHelper(3, third_derivatives_output, k, 0, output_type);
+          third_derivatives_output << "=" << eq + 1 << ";" << endl;
 
-      sparseHelper(3, third_derivatives_output, k, 1, output_type);
-      third_derivatives_output << "=" << ref_col + 1 << ";" << endl;
+          sparseHelper(3, third_derivatives_output, k, 1, output_type);
+          third_derivatives_output << "=" << ref_col + 1 << ";" << endl;
 
-      sparseHelper(3, third_derivatives_output, k, 2, output_type);
-      third_derivatives_output << "=";
-      d3->writeOutput(third_derivatives_output, output_type, temporary_terms, tef_terms);
-      third_derivatives_output << ";" << endl;
+          sparseHelper(3, third_derivatives_output, k, 2, output_type);
+          third_derivatives_output << "=";
+          d3->writeOutput(third_derivatives_output, output_type, temporary_terms, tef_terms);
+          third_derivatives_output << ";" << endl;
+        }
 
-      // Compute the column numbers for the 5 other permutations of (id1,id2,id3) and store them in a set (to avoid duplicates if two indexes are equal)
+      // Compute the column numbers for the 5 other permutations of (id1,id2,id3)
+      // and store them in a set (to avoid duplicates if two indexes are equal)
       set<int> cols;
       cols.insert(id1 * hessianColsNbr + id3 * dynJacobianColsNbr + id2);
       cols.insert(id2 * hessianColsNbr + id1 * dynJacobianColsNbr + id3);
@@ -2205,24 +2259,28 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
       int k2 = 1; // Keeps the offset of the permutation relative to k
       for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
         if (*it2 != ref_col)
-          {
-            sparseHelper(3, third_derivatives_output, k+k2, 0, output_type);
-            third_derivatives_output << "=" << eq + 1 << ";" << endl;
+          if (output_type == oJuliaDynamicModel)
+            third_derivatives_output << "  @inbounds g3[" << eq + 1 << "," << *it2 + 1 << "] = "
+                                     << for_sym.str() << endl;
+          else
+            {
+              sparseHelper(3, third_derivatives_output, k+k2, 0, output_type);
+              third_derivatives_output << "=" << eq + 1 << ";" << endl;
 
-            sparseHelper(3, third_derivatives_output, k+k2, 1, output_type);
-            third_derivatives_output << "=" << *it2 + 1 << ";" << endl;
+              sparseHelper(3, third_derivatives_output, k+k2, 1, output_type);
+              third_derivatives_output << "=" << *it2 + 1 << ";" << endl;
 
-            sparseHelper(3, third_derivatives_output, k+k2, 2, output_type);
-            third_derivatives_output << "=";
-            sparseHelper(3, third_derivatives_output, k, 2, output_type);
-            third_derivatives_output << ";" << endl;
+              sparseHelper(3, third_derivatives_output, k+k2, 2, output_type);
+              third_derivatives_output << "=";
+              sparseHelper(3, third_derivatives_output, k, 2, output_type);
+              third_derivatives_output << ";" << endl;
 
-            k2++;
-          }
+              k2++;
+            }
       k += k2;
     }
 
-  if (!use_dll)
+  if (output_type == oMatlabDynamicModel)
     {
       DynamicOutput << "%" << endl
                     << "% Model equations" << endl
@@ -2230,6 +2288,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
                     << endl
                     << "residual = zeros(" << nrows << ", 1);" << endl
                     << model_output.str()
+                    << model_eq_output.str()
         // Writing initialization instruction for matrix g1
                     << "if nargout >= 2," << endl
                     << "  g1 = zeros(" << nrows << ", " << dynJacobianColsNbr << ");" << endl
@@ -2271,7 +2330,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
 
       DynamicOutput << "end" << endl;
     }
-  else
+  else if (output_type == oCDynamicModel)
     {
       DynamicOutput << "void Dynamic(double *y, double *x, int nb_row_x, double *params, double *steady_state, int it_, double *residual, double *g1, double *v2, double *v3)" << endl
                     << "{" << endl
@@ -2279,6 +2338,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
                     << endl
                     << "  /* Residual equations */" << endl
                     << model_output.str()
+                    << model_eq_output.str()
                     << "  /* Jacobian  */" << endl
                     << "  if (g1 == NULL)" << endl
                     << "    return;" << endl
@@ -2307,10 +2367,106 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
 
       DynamicOutput << "}" << endl << endl;
     }
+  else
+    {
+      ostringstream comments;
+      comments << "## Function Arguments" << endl
+               << endl
+               << "## Input" << endl
+               << " 1 y:            Array{Float64, num_dynamic_vars, 1}             Vector of endogenous variables in the order stored" << endl
+               << "                                                                 in model.lead_lag_incidence; see the manual" << endl
+               << " 2 x:            Array{Float64, nperiods, length(model.exo)}     Matrix of exogenous variables (in declaration order)" << endl
+               << "                                                                 for all simulation periods" << endl
+               << " 3 params:       Array{Float64, length(model.param), 1}          Vector of parameter values in declaration order" << endl
+               << " 4 steady_state:" << endl
+               << " 5 it_:          Int                                             Time period for exogenous variables for which to evaluate the model" << endl
+               << endl
+               << "## Output" << endl
+               << " 6 residual:     Array(Float64, model.eq_nbr, 1)                 Vector of residuals of the dynamic model equations in" << endl
+               << "                                                                 order of declaration of the equations." << endl;
+
+      DynamicOutput << "function dynamic!(y::Vector{Float64}, x::Matrix{Float64}, "
+                    << "params::Vector{Float64}," << endl
+                    << "                  steady_state::Vector{Float64}, it_::Int, "
+                    << "residual::Vector{Float64})" << endl
+                    << "#=" << endl << comments.str() << "=#" << endl
+                    << "  @assert size(y) == " << dynJacobianColsNbr << endl
+                    << "  @assert size(params) == " << symbol_table.param_nbr() << endl
+                    << "  @assert size(residual) == " << nrows << endl
+                    << "  #" << endl
+                    << "  # Model equations" << endl
+                    << "  #" << endl
+                    << model_output.str()
+                    << model_eq_output.str()
+                    << "end" << endl << endl
+                    << "function dynamic!(y::Vector{Float64}, x::Matrix{Float64}, "
+                    << "params::Vector{Float64}," << endl
+                    << "                  steady_state::Vector{Float64}, it_::Int, "
+                    << "residual::Vector{Float64}," << endl
+                    << "                  g1::Matrix{Float64})" << endl;
+
+      comments << " 7 g1:           Array(Float64, model.eq_nbr, num_dynamic_vars)  Jacobian matrix of the dynamic model equations;" << endl
+               << "                                                                 rows: equations in order of declaration" << endl
+               << "                                                                 columns: variables in order stored in M_.lead_lag_incidence" << endl;
+
+      DynamicOutput << "#=" << endl << comments.str() << "=#" << endl
+                    << "  @assert size(g1) == (" << nrows << ", " << dynJacobianColsNbr << ")" << endl
+                    << "  fill!(g1, 0.0)" << endl
+                    << "  dynamic!(y, x, params, steady_state, it_, residual)" << endl
+                    << model_output.str()
+                    << "  #" << endl
+                    << "  # Jacobian matrix" << endl
+                    << "  #" << endl
+                    << jacobian_output.str()
+                    << "end" << endl << endl
+                    << "function dynamic!(y::Vector{Float64}, x::Matrix{Float64}, "
+                    << "params::Vector{Float64}," << endl
+                    << "                  steady_state::Vector{Float64}, it_::Int, "
+                    << "residual::Vector{Float64}," << endl
+                    << "                  g1::Matrix{Float64}, g2::Matrix{Float64})" << endl;
+
+      comments << " 8 g2:           spzeros(model.eq_nbr, (num_dynamic_vars)^2)     Hessian matrix of the dynamic model equations;" << endl
+               << "                                                                 rows: equations in order of declaration" << endl
+               << "                                                                 columns: variables in order stored in M_.lead_lag_incidence" << endl;
+
+      DynamicOutput << "#=" << endl << comments.str() << "=#" << endl
+                    << "  @assert size(g2) == (" << nrows << ", " << hessianColsNbr << ")" << endl
+                    << "  dynamic!(y, x, params, steady_state, it_, residual, g1)" << endl;
+      if (second_derivatives.size())
+        DynamicOutput << model_output.str()
+                      << "  #" << endl
+                      << "  # Hessian matrix" << endl
+                      << "  #" << endl
+                      << hessian_output.str();
+
+      // Initialize g3 matrix
+      int ncols = hessianColsNbr * dynJacobianColsNbr;
+      DynamicOutput << "end" << endl << endl
+                    << "function dynamic!(y::Vector{Float64}, x::Matrix{Float64}, "
+                    << "params::Vector{Float64}," << endl
+                    << "                  steady_state::Vector{Float64}, it_::Int, "
+                    << "residual::Vector{Float64}," << endl
+                    << "                  g1::Matrix{Float64}, g2::Matrix{Float64}, g3::Matrix{Float64})" << endl;
+
+      comments << " 9 g3:           spzeros(model.eq_nbr, (num_dynamic_vars)^3)     Third order derivative matrix of the dynamic model equations;" << endl
+               << "                                                                 rows: equations in order of declaration" << endl
+               << "                                                                 columns: variables in order stored in M_.lead_lag_incidence" << endl;
+
+      DynamicOutput << "#=" << endl << comments.str() << "=#" << endl
+                    << "  @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl
+                    << "  dynamic!(y, x, params, steady_state, it_, residual, g1, g2)" << endl;
+      if (third_derivatives.size())
+        DynamicOutput << model_output.str()
+                      << "  #" << endl
+                      << "  # Third order derivatives" << endl
+                      << "  #" << endl
+                      << third_derivatives_output.str();
+      DynamicOutput << "end" << endl;
+    }
 }
 
 void
-DynamicModel::writeOutput(ostream &output, const string &basename, bool block_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present) const
+DynamicModel::writeOutput(ostream &output, const string &basename, bool block_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present, bool julia) const
 {
   /* Writing initialisation for M_.lead_lag_incidence matrix
      M_.lead_lag_incidence is a matrix with as many columns as there are
@@ -2321,7 +2477,20 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
      model at a given period.
   */
 
-  output << "M_.lead_lag_incidence = [";
+  string modstruct;
+  string outstruct;
+  if (julia)
+    {
+      modstruct = "model.";
+      outstruct = "output.";
+    }
+  else
+    {
+      modstruct = "M_.";
+      outstruct = "oo_.";
+    }
+
+  output << modstruct << "lead_lag_incidence = [";
   // Loop on endogenous variables
   int nstatic = 0, 
       nfwrd   = 0,
@@ -2373,26 +2542,41 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
       output << ";";
     }
   output << "]';" << endl;
-  output << "M_.nstatic = " << nstatic << ";" << endl
-         << "M_.nfwrd   = " << nfwrd   << ";" << endl
-         << "M_.npred   = " << npred   << ";" << endl
-         << "M_.nboth   = " << nboth   << ";" << endl
-         << "M_.nsfwrd   = " << nfwrd+nboth   << ";" << endl
-         << "M_.nspred   = " << npred+nboth   << ";" << endl
-         << "M_.ndynamic   = " << npred+nboth+nfwrd << ";" << endl;
+  output << modstruct << "nstatic = " << nstatic << ";" << endl
+         << modstruct << "nfwrd   = " << nfwrd   << ";" << endl
+         << modstruct << "npred   = " << npred   << ";" << endl
+         << modstruct << "nboth   = " << nboth   << ";" << endl
+         << modstruct << "nsfwrd   = " << nfwrd+nboth   << ";" << endl
+         << modstruct << "nspred   = " << npred+nboth   << ";" << endl
+         << modstruct << "ndynamic   = " << npred+nboth+nfwrd << ";" << endl;
 
   // Write equation tags
-  output << "M_.equations_tags = {" << endl;
-  for (size_t i = 0; i < equation_tags.size(); i++)
-    output << "  " << equation_tags[i].first + 1 << " , '"
-           << equation_tags[i].second.first << "' , '"
-           << equation_tags[i].second.second << "' ;" << endl;
-  output << "};" << endl;
+  if (julia)
+    {
+      output << modstruct << "equation_tags = [" << endl;
+      for (size_t i = 0; i < equation_tags.size(); i++)
+        output << "                       EquationTag("
+               << equation_tags[i].first + 1 << " , \""
+               << equation_tags[i].second.first << "\" , \""
+               << equation_tags[i].second.second << "\")" << endl;
+      output << "                      ]" << endl;
+    }
+  else
+    {
+      output << modstruct << "equations_tags = {" << endl;
+      for (size_t i = 0; i < equation_tags.size(); i++)
+        output << "  " << equation_tags[i].first + 1 << " , '"
+               << equation_tags[i].second.first << "' , '"
+               << equation_tags[i].second.second << "' ;" << endl;
+      output << "};" << endl;
+    }
 
   /* Say if static and dynamic models differ (because of [static] and [dynamic]
      equation tags) */
-  output << "M_.static_and_dynamic_models_differ = "
-         << (static_only_equations.size() > 0 ? "1" : "0")
+  output << modstruct << "static_and_dynamic_models_differ = "
+         << (static_only_equations.size() > 0 ?
+             (julia ? "true" : "1") :
+             (julia ? "false" : "0"))
          << ";" << endl;
 
   //In case of sparse model, writes the block_decomposition structure of the model
@@ -2636,14 +2820,14 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
           output << "block_structure.block(" << block+1 << ").n_backward = " << n_backward << ";\n";
           output << "block_structure.block(" << block+1 << ").n_mixed = " << n_mixed << ";\n";
         }
-      output << "M_.block_structure.block = block_structure.block;\n";
+      output << modstruct << "block_structure.block = block_structure.block;\n";
       string cst_s;
       int nb_endo = symbol_table.endo_nbr();
-      output << "M_.block_structure.variable_reordered = [";
+      output << modstruct << "block_structure.variable_reordered = [";
       for (int i = 0; i < nb_endo; i++)
         output << " " << variable_reordered[i]+1;
       output << "];\n";
-      output << "M_.block_structure.equation_reordered = [";
+      output << modstruct << "block_structure.equation_reordered = [";
       for (int i = 0; i < nb_endo; i++)
         output << " " << equation_reordered[i]+1;
       output << "];\n";
@@ -2677,8 +2861,8 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
               if (prev_lag != -1000000)
                 output << "];\n";
               prev_lag = it->first.first;
-              output << "M_.block_structure.incidence(" << max_endo_lag+it->first.first+1 << ").lead_lag = " << prev_lag << ";\n";
-              output << "M_.block_structure.incidence(" << max_endo_lag+it->first.first+1 << ").sparse_IM = [";
+              output << modstruct << "block_structure.incidence(" << max_endo_lag+it->first.first+1 << ").lead_lag = " << prev_lag << ";\n";
+              output << modstruct << "block_structure.incidence(" << max_endo_lag+it->first.first+1 << ").sparse_IM = [";
             }
           output << it->first.second.first+1 << " " << it->first.second.second+1 << ";\n";
         }
@@ -2696,7 +2880,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
               n_obs--;
 
           int n = n_obs + n_state;
-          output << "M_.nobs_non_statevar = " << n_obs << ";" << endl;
+          output << modstruct << "nobs_non_statevar = " << n_obs << ";" << endl;
           int nb_diag = 0;
           //map<pair<int,int>, int>::const_iterator  row_state_var_incidence_it = row_state_var_incidence.begin();
 
@@ -2783,11 +2967,11 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
                 i_nz_state_var[lp + i] = lp + nze; 
               lp += nze; 
             }
-          output << "M_.nz_state_var = [";
+          output << modstruct << "nz_state_var = [";
           for (unsigned int i = 0; i < lp; i++)
             output << i_nz_state_var[i] << " ";
           output << "];" << endl;
-          output << "M_.n_diag = " << nb_diag << ";" << endl;
+          output << modstruct << "n_diag = " << nb_diag << ";" << endl;
           KF_index_file.write(reinterpret_cast<char *>(&nb_diag), sizeof(nb_diag));
           
           
@@ -2831,7 +3015,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
             KF_index_file.write(reinterpret_cast<char *>(&(*it)), sizeof(index_KF));      
           KF_index_file.close();
         }
-        output << "M_.state_var = [";
+        output << modstruct << "state_var = [";
 
         for (vector<int>::const_iterator it=state_var.begin(); it != state_var.end(); it++)
           output << *it << " ";
@@ -2839,30 +3023,36 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
     }
 
   // Writing initialization for some other variables
-  output << "M_.exo_names_orig_ord = [1:" << symbol_table.exo_nbr() << "];" << endl
-         << "M_.maximum_lag = " << max_lag << ";" << endl
-         << "M_.maximum_lead = " << max_lead << ";" << endl;
+  if (!julia)
+    output << modstruct << "exo_names_orig_ord = [1:" << symbol_table.exo_nbr() << "];" << endl;
+  else
+    output << modstruct << "exo_names_orig_ord = collect(1:" << symbol_table.exo_nbr() << ");" << endl;
 
-  output << "M_.maximum_endo_lag = " << max_endo_lag << ";" << endl
-         << "M_.maximum_endo_lead = " << max_endo_lead << ";" << endl
-         << "oo_.steady_state = zeros(" << symbol_table.endo_nbr() << ", 1);" << endl;
+  output << modstruct << "maximum_lag = " << max_lag << ";" << endl
+         << modstruct << "maximum_lead = " << max_lead << ";" << endl;
 
-  output << "M_.maximum_exo_lag = " << max_exo_lag << ";" << endl
-         << "M_.maximum_exo_lead = " << max_exo_lead << ";" << endl
-         << "oo_.exo_steady_state = zeros(" << symbol_table.exo_nbr() << ", 1);" << endl;
+  output << modstruct << "maximum_endo_lag = " << max_endo_lag << ";" << endl
+         << modstruct << "maximum_endo_lead = " << max_endo_lead << ";" << endl
+         << outstruct << "steady_state = zeros(" << symbol_table.endo_nbr() << ", 1);" << endl;
+
+  output << modstruct << "maximum_exo_lag = " << max_exo_lag << ";" << endl
+         << modstruct << "maximum_exo_lead = " << max_exo_lead << ";" << endl
+         << outstruct << "exo_steady_state = zeros(" << symbol_table.exo_nbr() << ", 1);" << endl;
 
   if (symbol_table.exo_det_nbr())
     {
-      output << "M_.maximum_exo_det_lag = " << max_exo_det_lag << ";" << endl
-             << "M_.maximum_exo_det_lead = " << max_exo_det_lead << ";" << endl
-             << "oo_.exo_det_steady_state = zeros(" << symbol_table.exo_det_nbr() << ", 1);" << endl;
+      output << modstruct << "maximum_exo_det_lag = " << max_exo_det_lag << ";" << endl
+             << modstruct << "maximum_exo_det_lead = " << max_exo_det_lead << ";" << endl
+             << outstruct << "exo_det_steady_state = zeros(" << symbol_table.exo_det_nbr() << ", 1);" << endl;
     }
 
-  output << "M_.params = NaN(" << symbol_table.param_nbr() << ", 1);" << endl;
+  output << modstruct << "params = " << (julia ? "fill(NaN, " : "NaN(")
+         << symbol_table.param_nbr() << ", 1);" << endl;
 
   // Write number of non-zero derivatives
   // Use -1 if the derivatives have not been computed
-  output << "M_.NNZDerivatives = [" << NNZDerivatives[0] << "; ";
+  output << modstruct << (julia ? "nnzderivatives" : "NNZDerivatives")
+                          << " = [" << NNZDerivatives[0] << "; ";
   if (order > 1)
     output << NNZDerivatives[1] << "; ";
   else
@@ -3306,7 +3496,7 @@ DynamicModel::collectBlockVariables()
 }
 
 void
-DynamicModel::writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, int order) const
+DynamicModel::writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, int order, bool julia) const
 {
   int r;
   string t_basename = basename + "_dynamic";
@@ -3330,6 +3520,8 @@ DynamicModel::writeDynamicFile(const string &basename, bool block, bool bytecode
     }
   else if (use_dll)
     writeDynamicCFile(t_basename, order);
+  else if (julia)
+    writeDynamicJuliaFile(basename);
   else
     writeDynamicMFile(t_basename);
 }
@@ -3740,7 +3932,7 @@ DynamicModel::testTrendDerivativesEqualToZero(const eval_context_t &eval_context
 }
 
 void
-DynamicModel::writeParamsDerivativesFile(const string &basename) const
+DynamicModel::writeParamsDerivativesFile(const string &basename, bool julia) const
 {
   if (!residuals_params_derivatives.size()
       && !residuals_params_second_derivatives.size()
@@ -3749,8 +3941,7 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
       && !hessian_params_derivatives.size())
     return;
 
-  string filename = basename + "_params_derivs.m";
-
+  string filename = julia ? basename + "DynamicParamsDerivs.jl" : basename + "_params_derivs.m";
   ofstream paramsDerivsFile;
   paramsDerivsFile.open(filename.c_str(), ios::out | ios::binary);
   if (!paramsDerivsFile.is_open())
@@ -3758,15 +3949,28 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
       cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
       exit(EXIT_FAILURE);
     }
-  paramsDerivsFile << "function [rp, gp, rpp, gpp, hp] = " << basename << "_params_derivs(y, x, params, steady_state, it_, ss_param_deriv, ss_param_2nd_deriv)" << endl
-                   << "%" << endl
-                   << "% Warning : this file is generated automatically by Dynare" << endl
-                   << "%           from model file (.mod)" << endl << endl;
+
+  ExprNodeOutputType output_type = (julia ? oJuliaDynamicModel : oMatlabDynamicModel);
+
+  if (!julia)
+    paramsDerivsFile << "function [rp, gp, rpp, gpp, hp] = " << basename << "_params_derivs(y, x, params, steady_state, it_, ss_param_deriv, ss_param_2nd_deriv)" << endl
+                     << "%" << endl
+                     << "% Warning : this file is generated automatically by Dynare" << endl
+                     << "%           from model file (.mod)" << endl << endl;
+  else
+    paramsDerivsFile << "module " << basename << "DynamicParamsDerivs" << endl
+                     << "#" << endl
+                     << "# NB: this file was automatically generated by Dynare" << endl
+                     << "#     from " << basename << ".mod" << endl
+                     << "#" << endl
+                     << "export params_derivs" << endl << endl
+                     << "function params_derivs(y, x, paramssteady_state, it_, "
+                     << "ss_param_deriv, ss_param_2nd_deriv)" << endl;
 
   deriv_node_temp_terms_t tef_terms;
-  writeModelLocalVariables(paramsDerivsFile, oMatlabDynamicModel, tef_terms);
+  writeModelLocalVariables(paramsDerivsFile, output_type, tef_terms);
 
-  writeTemporaryTerms(params_derivs_temporary_terms, paramsDerivsFile, oMatlabDynamicModel, tef_terms);
+  writeTemporaryTerms(params_derivs_temporary_terms, paramsDerivsFile, output_type, tef_terms);
 
   // Write parameter derivative
   paramsDerivsFile << "rp = zeros(" << equation_number() << ", "
@@ -3781,8 +3985,9 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
 
       int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
 
-      paramsDerivsFile << "rp(" << eq+1 << ", " << param_col << ") = ";
-      d1->writeOutput(paramsDerivsFile, oMatlabDynamicModel, params_derivs_temporary_terms, tef_terms);
+      paramsDerivsFile << "rp" << LEFT_ARRAY_SUBSCRIPT(output_type) << eq+1 << ", " << param_col
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << " = ";
+      d1->writeOutput(paramsDerivsFile, output_type, params_derivs_temporary_terms, tef_terms);
       paramsDerivsFile << ";" << endl;
     }
 
@@ -3801,13 +4006,15 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
       int var_col = getDynJacobianCol(var) + 1;
       int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
 
-      paramsDerivsFile << "gp(" << eq+1 << ", " << var_col << ", " << param_col << ") = ";
-      d2->writeOutput(paramsDerivsFile, oMatlabDynamicModel, params_derivs_temporary_terms, tef_terms);
+      paramsDerivsFile << "gp" << LEFT_ARRAY_SUBSCRIPT(output_type) << eq+1 << ", " << var_col
+                       << ", " << param_col << RIGHT_ARRAY_SUBSCRIPT(output_type) << " = ";
+      d2->writeOutput(paramsDerivsFile, output_type, params_derivs_temporary_terms, tef_terms);
       paramsDerivsFile << ";" << endl;
     }
 
-  // If nargout >= 3...
-  paramsDerivsFile << "if nargout >= 3" << endl;
+  if (!julia)
+    // If nargout >= 3...
+    paramsDerivsFile << "if nargout >= 3" << endl;
 
   // Write parameter second derivatives (only if nargout >= 3)
   paramsDerivsFile << "rpp = zeros(" << residuals_params_second_derivatives.size()
@@ -3825,11 +4032,15 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
       int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
       int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1;
 
-      paramsDerivsFile << "rpp(" << i << ",1)=" << eq+1 << ";" << endl
-                       << "rpp(" << i << ",2)=" << param1_col << ";" << endl
-                       << "rpp(" << i << ",3)=" << param2_col << ";" << endl
-                       << "rpp(" << i << ",4)=";
-      d2->writeOutput(paramsDerivsFile, oMatlabDynamicModel, params_derivs_temporary_terms, tef_terms);
+      paramsDerivsFile << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
+                       << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
+                       << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
+                       << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+      d2->writeOutput(paramsDerivsFile, output_type, params_derivs_temporary_terms, tef_terms);
       paramsDerivsFile << ";" << endl;
     }
 
@@ -3851,18 +4062,24 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
       int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
       int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1;
 
-      paramsDerivsFile << "gpp(" << i << ",1)=" << eq+1 << ";" << endl
-                       << "gpp(" << i << ",2)=" << var_col << ";" << endl
-                       << "gpp(" << i << ",3)=" << param1_col << ";" << endl
-                       << "gpp(" << i << ",4)=" << param2_col << ";" << endl
-                       << "gpp(" << i << ",5)=";
-      d2->writeOutput(paramsDerivsFile, oMatlabDynamicModel, params_derivs_temporary_terms, tef_terms);
+      paramsDerivsFile << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
+                       << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var_col << ";" << endl
+                       << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
+                       << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
+                       << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+      d2->writeOutput(paramsDerivsFile, output_type, params_derivs_temporary_terms, tef_terms);
       paramsDerivsFile << ";" << endl;
     }
 
-  // If nargout >= 5...
-  paramsDerivsFile << "end" << endl
-                   << "if nargout >= 5" << endl;
+  if (!julia)
+    // If nargout >= 5...
+    paramsDerivsFile << "end" << endl
+                     << "if nargout >= 5" << endl;
 
   // Write hessian derivatives (only if nargout >= 5)
   paramsDerivsFile << "hp = zeros(" << hessian_params_derivatives.size() << ",5);" << endl;
@@ -3881,15 +4098,22 @@ DynamicModel::writeParamsDerivativesFile(const string &basename) const
       int var2_col = getDynJacobianCol(var2) + 1;
       int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
 
-      paramsDerivsFile << "hp(" << i << ",1)=" << eq+1 << ";" << endl
-                       << "hp(" << i << ",2)=" << var1_col << ";" << endl
-                       << "hp(" << i << ",3)=" << var2_col << ";" << endl
-                       << "hp(" << i << ",4)=" << param_col << ";" << endl
-                       << "hp(" << i << ",5)=";
-      d2->writeOutput(paramsDerivsFile, oMatlabDynamicModel, params_derivs_temporary_terms, tef_terms);
+      paramsDerivsFile << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
+                       << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
+                       << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
+                       << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
+                       << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+      d2->writeOutput(paramsDerivsFile, output_type, params_derivs_temporary_terms, tef_terms);
       paramsDerivsFile << ";" << endl;
     }
 
+  if (julia)
+    paramsDerivsFile << "(rp, gp, rpp, gpp, hp)" << endl;
   paramsDerivsFile << "end" << endl
                    << "end" << endl;
   paramsDerivsFile.close();
diff --git a/preprocessor/DynamicModel.hh b/preprocessor/DynamicModel.hh
index 61b758e229c149c158ccbc7bf926c961285cb718..fbb18b222022faffa8d420d2b5852b9e911f2284 100644
--- a/preprocessor/DynamicModel.hh
+++ b/preprocessor/DynamicModel.hh
@@ -76,6 +76,8 @@ private:
 
   //! Writes dynamic model file (Matlab version)
   void writeDynamicMFile(const string &dynamic_basename) const;
+  //! Writes dynamic model file (Julia version)
+  void writeDynamicJuliaFile(const string &dynamic_basename) const;
   //! Writes dynamic model file (C version)
   /*! \todo add third derivatives handling */
   void writeDynamicCFile(const string &dynamic_basename, const int order) const;
@@ -83,7 +85,7 @@ private:
   void writeSparseDynamicMFile(const string &dynamic_basename, const string &basename) const;
   //! Writes the dynamic model equations and its derivatives
   /*! \todo add third derivatives handling in C output */
-  void writeDynamicModel(ostream &DynamicOutput, bool use_dll) const;
+  void writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia) const;
   //! Writes the Block reordred structure of the model in M output
   void writeModelEquationsOrdered_M(const string &dynamic_basename) const;
   //! Writes the code of the Block reordred structure of the model in virtual machine bytecode
@@ -213,15 +215,15 @@ public:
   void computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, bool paramsDerivatives,
                      const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode);
   //! Writes model initialization and lead/lag incidence matrix to output
-  void writeOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present) const;
+  void writeOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present, bool julia) const;
 
   //! Adds informations for simulation in a binary file
   void Write_Inf_To_Bin_File_Block(const string &dynamic_basename, const string &bin_basename,
                                    const int &num, int &u_count_int, bool &file_open, bool is_two_boundaries) const;
   //! Writes dynamic model file
-  void writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, int order) const;
+  void writeDynamicFile(const string &basename, bool block, bool bytecode, bool use_dll, int order, bool julia) const;
   //! Writes file containing parameters derivatives
-  void writeParamsDerivativesFile(const string &basename) const;
+  void writeParamsDerivativesFile(const string &basename, bool julia) const;
 
   //! Converts to static model (only the equations)
   /*! It assumes that the static model given in argument has just been allocated */
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index b4f2a5668fef8f784a00ab85cf36fa87d6a300ef..a1b92b1e243c24b95544b5ebdddcbcf7ebe17a2f 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -783,7 +783,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
 <DYNARE_STATEMENT>max_dim_cova_group {return token::MAX_DIM_COVA_GROUP;}
 <DYNARE_STATEMENT>gsa_sample_file {return token::GSA_SAMPLE_FILE;}
 
-<DYNARE_STATEMENT,DYNARE_BLOCK>[A-Za-z_][A-Za-z0-9_]* {
+<DYNARE_STATEMENT,DYNARE_BLOCK>[A-Za-z_\x80-\xf3][A-Za-z0-9_\x80-\xf3]* {
   yylval->string_val = new string(yytext);
   return token::NAME;
 }
@@ -845,7 +845,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
     element in initval (in which case Dynare recognizes the matrix name as an external
     function symbol), and may want to modify the matrix later with Matlab statements.
  */
-<INITIAL>[A-Za-z_][A-Za-z0-9_]* {
+<INITIAL>[A-Za-z_\x80-\xf3][A-Za-z0-9_\x80-\xf3]* {
   if (driver.symbol_exists_and_is_not_modfile_local_or_external_function(yytext))
     {
       BEGIN DYNARE_STATEMENT;
@@ -863,7 +863,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
  /* For joint prior statement, match [symbol, symbol, ...]
    If no match, begin native and push everything back on stack
  */
-<INITIAL>\[([[:space:]]*[A-Za-z_][A-Za-z0-9_]*[[:space:]]*,{1}[[:space:]]*)*([[:space:]]*[A-Za-z_][A-Za-z0-9_]*[[:space:]]*){1}\] {
+<INITIAL>\[([[:space:]]*[A-Za-z_\x80-\xf3][A-Za-z0-9_\x80-\xf3]*[[:space:]]*,{1}[[:space:]]*)*([[:space:]]*[A-Za-z_\x80-\xf3][A-Za-z0-9_\x80-\xf3]*[[:space:]]*){1}\] {
   string yytextcpy = string(yytext);
   yytextcpy.erase(remove(yytextcpy.begin(), yytextcpy.end(), '['), yytextcpy.end());
   yytextcpy.erase(remove(yytextcpy.begin(), yytextcpy.end(), ']'), yytextcpy.end());
diff --git a/preprocessor/DynareMain.cc b/preprocessor/DynareMain.cc
index e78fe4f490c963299f2d862b629c21ba0e0f9358..8a4c8e39dbf73c7aa1b30b4f7eef9e23392c3bf0 100644
--- a/preprocessor/DynareMain.cc
+++ b/preprocessor/DynareMain.cc
@@ -55,7 +55,7 @@ usage()
 {
   cerr << "Dynare usage: dynare mod_file [debug] [noclearall] [onlyclearglobals] [savemacro[=macro_file]] [onlymacro] [nolinemacro] [notmpterms] [nolog] [warn_uninit]"
        << " [console] [nograph] [nointeractive] [parallel[=cluster_name]] [conffile=parallel_config_path_and_filename] [parallel_slave_open_mode] [parallel_test]"
-       << " [-D<variable>[=<value>]] [-I/path] [nostrict] [fast] [minimal_workspace] [output=dynamic|first|second|third] [language=C|C++]"
+       << " [-D<variable>[=<value>]] [-I/path] [nostrict] [fast] [minimal_workspace] [output=dynamic|first|second|third] [language=C|C++|julia]"
 #if defined(_WIN32) || defined(__CYGWIN32__)
        << " [cygwin] [msvc]"
 #endif
diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc
index d2b779a449caa423998f8c9efcc0d03d6eeb6f1c..7b52a0bdaf6b06591d160939e2b3ac6969279614 100644
--- a/preprocessor/ExprNode.cc
+++ b/preprocessor/ExprNode.cc
@@ -654,6 +654,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
     case eEndogenous:
       switch (output_type)
         {
+        case oJuliaDynamicModel:
         case oMatlabDynamicModel:
         case oCDynamicModel:
           i = datatree.getDynJacobianCol(datatree.getDerivID(symb_id, lag)) + ARRAY_SUBSCRIPT_OFFSET(output_type);
@@ -664,6 +665,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
           output <<  "y" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
         case oCStaticModel:
+        case oJuliaStaticModel:
         case oMatlabStaticModel:
         case oMatlabStaticModelSparse:
           i = tsid + ARRAY_SUBSCRIPT_OFFSET(output_type);
@@ -681,15 +683,17 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         case oMatlabOutsideModel:
           output << "oo_.steady_state(" << tsid + 1 << ")";
           break;
+        case oJuliaDynamicSteadyStateOperator:
         case oMatlabDynamicSteadyStateOperator:
         case oMatlabDynamicSparseSteadyStateOperator:
-          output << "steady_state(" << tsid + 1 << ")";
+          output << "steady_state" << LEFT_ARRAY_SUBSCRIPT(output_type) << tsid + 1 << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
         case oCDynamicSteadyStateOperator:
           output << "steady_state[" << tsid << "]";
           break;
+        case oJuliaSteadyStateFile:
         case oSteadyStateFile:
-          output << "ys_(" << tsid + 1 << ")";
+          output << "ys_" << LEFT_ARRAY_SUBSCRIPT(output_type) << tsid + 1 << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
         case oCSteadyStateFile:
           output << "ys_[" << tsid << "]";
@@ -704,14 +708,18 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       i = tsid + ARRAY_SUBSCRIPT_OFFSET(output_type);
       switch (output_type)
         {
+        case oJuliaDynamicModel:
         case oMatlabDynamicModel:
         case oMatlabDynamicModelSparse:
           if (lag > 0)
-            output <<  "x(it_+" << lag << ", " << i << ")";
+            output <<  "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_+" << lag << ", " << i
+                   << RIGHT_ARRAY_SUBSCRIPT(output_type);
           else if (lag < 0)
-            output <<  "x(it_" << lag << ", " << i << ")";
+            output <<  "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_" << lag << ", " << i
+                   << RIGHT_ARRAY_SUBSCRIPT(output_type);
           else
-            output <<  "x(it_, " << i << ")";
+            output <<  "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_, " << i
+                   << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
         case oCDynamicModel:
         case oCDynamic2Model:
@@ -723,6 +731,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
             output <<  "x[it_" << lag << "+" << i << "*nb_row_x]";
           break;
         case oCStaticModel:
+        case oJuliaStaticModel:
         case oMatlabStaticModel:
         case oMatlabStaticModelSparse:
           output << "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
@@ -734,8 +743,9 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         case oMatlabDynamicSteadyStateOperator:
           output <<  "oo_.exo_steady_state(" << i << ")";
           break;
+        case oJuliaSteadyStateFile:
         case oSteadyStateFile:
-          output << "exo_(" << i << ")";
+          output << "exo_" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
         case oCSteadyStateFile:
           output << "exo_[" << i - 1 << "]";
@@ -750,14 +760,18 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       i = tsid + datatree.symbol_table.exo_nbr() + ARRAY_SUBSCRIPT_OFFSET(output_type);
       switch (output_type)
         {
+        case oJuliaDynamicModel:
         case oMatlabDynamicModel:
         case oMatlabDynamicModelSparse:
           if (lag > 0)
-            output <<  "x(it_+" << lag << ", " << i << ")";
+            output <<  "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_+" << lag << ", " << i
+                   << RIGHT_ARRAY_SUBSCRIPT(output_type);
           else if (lag < 0)
-            output <<  "x(it_" << lag << ", " << i << ")";
+            output <<  "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_" << lag << ", " << i
+                   << RIGHT_ARRAY_SUBSCRIPT(output_type);
           else
-            output <<  "x(it_, " << i << ")";
+            output <<  "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << "it_, " << i
+                   << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
         case oCDynamicModel:
         case oCDynamic2Model:
@@ -769,6 +783,7 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
             output <<  "x[it_" << lag << "+" << i << "*nb_row_x]";
           break;
         case oCStaticModel:
+        case oJuliaStaticModel:
         case oMatlabStaticModel:
         case oMatlabStaticModelSparse:
           output << "x" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
@@ -780,8 +795,9 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         case oMatlabDynamicSteadyStateOperator:
           output <<  "oo_.exo_det_steady_state(" << tsid + 1 << ")";
           break;
+        case oJuliaSteadyStateFile:
         case oSteadyStateFile:
-          output << "exo_(" << i << ")";
+          output << "exo_" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
         case oCSteadyStateFile:
           output << "exo_[" << i - 1 << "]";
@@ -1836,6 +1852,9 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         case oCDynamicModel:
           new_output_type = oCDynamicSteadyStateOperator;
           break;
+        case oJuliaDynamicModel:
+          new_output_type = oJuliaDynamicSteadyStateOperator;
+          break;
         case oMatlabDynamicModelSparse:
           new_output_type = oMatlabDynamicSparseSteadyStateOperator;
           break;
@@ -2939,7 +2958,10 @@ BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         unpackPowerDeriv()->writeOutput(output, output_type, temporary_terms, tef_terms);
       else
         {
-          output << "getPowerDeriv(";
+          if (output_type == oJuliaStaticModel || output_type == oJuliaDynamicModel)
+            output << "get_power_deriv(";
+          else
+            output << "getPowerDeriv(";
           arg1->writeOutput(output, output_type, temporary_terms, tef_terms);
           output << ",";
           arg2->writeOutput(output, output_type, temporary_terms, tef_terms);
@@ -3047,7 +3069,7 @@ BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         output << "~=";
       else
         {
-          if (IS_C(output_type))
+          if (IS_C(output_type) || IS_JULIA(output_type))
             output << "!=";
           else
             output << "\\neq ";
@@ -4809,7 +4831,8 @@ ExternalFunctionNode::writeOutput(ostream &output, ExprNodeOutputType output_typ
                                   deriv_node_temp_terms_t &tef_terms) const
 {
   if (output_type == oMatlabOutsideModel || output_type == oSteadyStateFile
-      || output_type == oCSteadyStateFile || IS_LATEX(output_type))
+      || output_type == oCSteadyStateFile || output_type == oJuliaSteadyStateFile
+      || IS_LATEX(output_type))
     {
       string name = IS_LATEX(output_type) ? datatree.symbol_table.getTeXName(symb_id)
         : datatree.symbol_table.getName(symb_id);
diff --git a/preprocessor/ExprNode.hh b/preprocessor/ExprNode.hh
index 394bf586a0eca6e7afb1adab54fb41cc16510d0d..a8f48e306dce6b72e91e99d711f448e848360529 100644
--- a/preprocessor/ExprNode.hh
+++ b/preprocessor/ExprNode.hh
@@ -67,6 +67,8 @@ enum ExprNodeOutputType
     oCDynamicModel,                               //!< C code, dynamic model
     oCDynamic2Model,                              //!< C code, dynamic model, alternative numbering of endogenous variables
     oCStaticModel,                                //!< C code, static model
+    oJuliaStaticModel,                            //!< Julia code, static model
+    oJuliaDynamicModel,                           //!< Julia code, dynamic model
     oMatlabOutsideModel,                          //!< Matlab code, outside model block (for example in initval)
     oLatexStaticModel,                            //!< LaTeX code, static model
     oLatexDynamicModel,                           //!< LaTeX code, dynamic model
@@ -74,8 +76,10 @@ enum ExprNodeOutputType
     oMatlabDynamicSteadyStateOperator,            //!< Matlab code, dynamic model, inside a steady state operator
     oMatlabDynamicSparseSteadyStateOperator,      //!< Matlab code, dynamic block decomposed model, inside a steady state operator
     oCDynamicSteadyStateOperator,                 //!< C code, dynamic model, inside a steady state operator
-    oSteadyStateFile,                              //!< Matlab code, in the generated steady state file
-    oCSteadyStateFile                             //!< C code, in the generated steady state file
+    oJuliaDynamicSteadyStateOperator,             //!< Julia code, dynamic model, inside a steady state operator
+    oSteadyStateFile,                             //!< Matlab code, in the generated steady state file
+    oCSteadyStateFile,                            //!< C code, in the generated steady state file
+    oJuliaSteadyStateFile                         //!< Julia code, in the generated steady state file
   };
 
 #define IS_MATLAB(output_type) ((output_type) == oMatlabStaticModel     \
@@ -87,6 +91,10 @@ enum ExprNodeOutputType
                                 || (output_type) == oMatlabDynamicSparseSteadyStateOperator \
                                 || (output_type) == oSteadyStateFile)
 
+#define IS_JULIA(output_type) ((output_type) == oJuliaStaticModel     \
+                               || (output_type) == oJuliaDynamicModel  \
+                               || (output_type) == oJuliaDynamicSteadyStateOperator)
+
 #define IS_C(output_type) ((output_type) == oCDynamicModel \
 			   || (output_type) == oCDynamic2Model \
 			   || (output_type) == oCStaticModel \
@@ -97,9 +105,9 @@ enum ExprNodeOutputType
                                || (output_type) == oLatexDynamicModel   \
                                || (output_type) == oLatexDynamicSteadyStateOperator)
 
-/* Equal to 1 for Matlab langage, or to 0 for C language. Not defined for LaTeX.
-   In Matlab, array indexes begin at 1, while they begin at 0 in C */
-#define ARRAY_SUBSCRIPT_OFFSET(output_type) ((int) IS_MATLAB(output_type))
+/* Equal to 1 for Matlab langage or Julia, or to 0 for C language. Not defined for LaTeX.
+   In Matlab and Julia, array indexes begin at 1, while they begin at 0 in C */
+#define ARRAY_SUBSCRIPT_OFFSET(output_type) ((int) (IS_MATLAB(output_type) || IS_JULIA(output_type)))
 
 // Left and right array subscript delimiters: '(' and ')' for Matlab, '[' and ']' for C
 #define LEFT_ARRAY_SUBSCRIPT(output_type) (IS_MATLAB(output_type) ? '(' : '[')
diff --git a/preprocessor/ExtendedPreprocessorTypes.hh b/preprocessor/ExtendedPreprocessorTypes.hh
index d9a210e7706492cfc6bae18c6880b90bfd0453b4..e0a955f2c8222b611512d31a296902b9e4e65fbd 100644
--- a/preprocessor/ExtendedPreprocessorTypes.hh
+++ b/preprocessor/ExtendedPreprocessorTypes.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2014 Dynare Team
+ * Copyright (C) 2014-2015 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -35,7 +35,7 @@ enum LanguageOutputType
     c,                                // outputs files for C
     cpp,                              // outputs files for C++
     cuda,                             // outputs files for CUDA (not yet implemented)
-    julia,                            // outputs files for Julia (not yet implemented)
+    julia,                            // outputs files for Julia
     python,                           // outputs files for Python (not yet implemented) (not yet implemented)
   };
 #endif
diff --git a/preprocessor/Makefile.am b/preprocessor/Makefile.am
index 40083fdcccd6c2d128f423b51bcd6ff5bcba1e45..9fd67c4fad4520ead0aa7602db1e5b3ad8d1ca16 100644
--- a/preprocessor/Makefile.am
+++ b/preprocessor/Makefile.am
@@ -76,7 +76,9 @@ all-local: $(PROGRAMS)
 	  ARCH="64"; \
 	fi; \
 	mkdir -p ../matlab/preprocessor$$ARCH ; \
-	cd ../matlab/preprocessor$$ARCH && $(LN_S) -f $(abs_srcdir)/$(PROGRAMS) $(PROGRAMS)
+	cd ../matlab/preprocessor$$ARCH && $(LN_S) -f $(abs_srcdir)/$(PROGRAMS) $(PROGRAMS) ; \
+	mkdir -p ../../julia/preprocessor$$ARCH ; \
+	cd ../../julia/preprocessor$$ARCH && $(LN_S) -f $(abs_srcdir)/$(PROGRAMS) $(PROGRAMS)
 
 if HAVE_DOXYGEN
 html-local:
diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc
index 1e8fb25b2813566dcfc47d0ef200716148f683d1..f14c7b479faa1608ace8a987e0ad569e2ce8c276 100644
--- a/preprocessor/ModFile.cc
+++ b/preprocessor/ModFile.cc
@@ -742,7 +742,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
 
   if (dynamic_model.equation_number() > 0)
     {
-      dynamic_model.writeOutput(mOutputFile, basename, block, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present);
+      dynamic_model.writeOutput(mOutputFile, basename, block, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present, false);
       if (!no_static)
         static_model.writeOutput(mOutputFile, block);
     }
@@ -817,18 +817,18 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
 	{
 	  if (!no_static)
 	    {
-	      static_model.writeStaticFile(basename, block, byte_code, use_dll);
-	      static_model.writeParamsDerivativesFile(basename);
+	      static_model.writeStaticFile(basename, block, byte_code, use_dll, false);
+	      static_model.writeParamsDerivativesFile(basename, false);
 	    }
 
-	  dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option);
-	  dynamic_model.writeParamsDerivativesFile(basename);
+	  dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option, false);
+	  dynamic_model.writeParamsDerivativesFile(basename, false);
 	}
 
       // Create steady state file
-      steady_state_model.writeSteadyStateFile(basename, mod_file_struct.ramsey_model_present);
+      steady_state_model.writeSteadyStateFile(basename, mod_file_struct.ramsey_model_present, false);
     }
-  
+
   cout << "done" << endl;
 }
 
@@ -843,6 +843,9 @@ ModFile::writeExternalFiles(const string &basename, FileOutputType output, Langu
     case cpp:
       writeExternalFilesCC(basename, output);
       break;
+    case julia:
+      writeExternalFilesJulia(basename, output);
+      break;
     default:
       cerr << "This case shouldn't happen. Contact the authors of Dynare" << endl;
       exit(EXIT_FAILURE);
@@ -855,10 +858,10 @@ ModFile::writeExternalFilesC(const string &basename, FileOutputType output) cons
   writeModelC(basename);
   steady_state_model.writeSteadyStateFileC(basename, mod_file_struct.ramsey_model_present);
 
-  dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option);
+  dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option, false);
 
   if (!no_static)
-    static_model.writeStaticFile(basename, false, false, true);
+    static_model.writeStaticFile(basename, false, false, true, false);
 
 
   //  static_model.writeStaticCFile(basename, block, byte_code, use_dll);
@@ -960,10 +963,10 @@ ModFile::writeExternalFilesCC(const string &basename, FileOutputType output) con
   writeModelCC(basename);
   steady_state_model.writeSteadyStateFileC(basename, mod_file_struct.ramsey_model_present);
 
-  dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option);
+  dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option, false);
 
   if (!no_static)
-    static_model.writeStaticFile(basename, false, false, true);
+    static_model.writeStaticFile(basename, false, false, true, false);
 
   //  static_model.writeStaticCFile(basename, block, byte_code, use_dll);
   //  static_model.writeParamsDerivativesFileC(basename, cuda);
@@ -1059,3 +1062,113 @@ ModFile::writeModelCC(const string &basename) const
   mOutputFile.close();
 }
 
+void
+ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output) const
+{
+  ofstream jlOutputFile;
+  if (basename.size())
+    {
+      string fname(basename);
+      fname += ".jl";
+      jlOutputFile.open(fname.c_str(), ios::out | ios::binary);
+      if (!jlOutputFile.is_open())
+        {
+          cerr << "ERROR: Can't open file " << fname
+               << " for writing" << endl;
+          exit(EXIT_FAILURE);
+        }
+    }
+  else
+    {
+      cerr << "ERROR: Missing file name" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  jlOutputFile << "module " << basename << endl
+               << "#" << endl
+               << "# NB: this file was automatically generated by Dynare" << endl
+               << "#     from " << basename << ".mod" << endl
+               << "#" << endl
+               << "using DynareModel" << endl
+               << "using DynareOptions" << endl
+               << "using DynareOutput" << endl
+               << "using Utils" << endl
+               << "using " << basename << "Static" << endl
+               << "using " << basename << "Dynamic" << endl
+               << "using " << basename << "SteadyState2" << endl << endl
+               << "export model" << endl;
+
+  // Write Output
+  jlOutputFile << endl
+               << "output = dynare_output()" << endl
+               << "output.dynare_version = \"" << PACKAGE_VERSION << "\"" << endl;
+
+  // Write Options
+  jlOutputFile << endl
+               << "options = dynare_options()" << endl
+               << "options.dynare_version = \"" << PACKAGE_VERSION << "\"" << endl;
+  if (linear == 1)
+    jlOutputFile << "options.linear = true" << endl;
+
+  // Write Model
+  jlOutputFile << endl
+               << "model = dynare_model()" << endl
+               << "model.fname = \"" << basename << "\"" << endl
+               << "model.dynare_version = \"" << PACKAGE_VERSION << "\"" << endl
+               << "model.sigma_e = zeros(Float64, " << symbol_table.exo_nbr() << ", "
+               << symbol_table.exo_nbr() << ")" << endl
+               << "model.correlation_matrix = ones(Float64, " << symbol_table.exo_nbr() << ", "
+               << symbol_table.exo_nbr() << ")" << endl
+               << "model.orig_eq_nbr = " << orig_eqn_nbr << endl
+               << "model.eq_nbr = " << dynamic_model.equation_number() << endl
+               << "model.ramsey_eq_nbr = " << ramsey_eqn_nbr << endl;
+
+  if (mod_file_struct.calibrated_measurement_errors)
+    jlOutputFile << "model.h = zeros(Float64,"
+                 << symbol_table.observedVariablesNbr() << ", "
+                 << symbol_table.observedVariablesNbr() << ");" << endl
+                 << "model.correlation_matrix_me = ones(Float64, "
+                 << symbol_table.observedVariablesNbr() << ", "
+                 << symbol_table.observedVariablesNbr() << ");" << endl;
+  else
+    jlOutputFile << "model.h = zeros(Float64, 1, 1)" << endl
+                 << "model.correlation_matrix_me = ones(Float64, 1, 1)" << endl;
+
+  cout << "Processing outputs ..." << endl;
+  symbol_table.writeJuliaOutput(jlOutputFile);
+
+  if (dynamic_model.equation_number() > 0)
+    {
+      dynamic_model.writeOutput(jlOutputFile, basename, false, false, false,
+                                mod_file_struct.order_option,
+                                mod_file_struct.estimation_present, true);
+      if (!no_static)
+        {
+          static_model.writeStaticFile(basename, false, false, false, true);
+          static_model.writeParamsDerivativesFile(basename, true);
+        }
+      dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll,
+                                     mod_file_struct.order_option, true);
+      dynamic_model.writeParamsDerivativesFile(basename, true);
+    }
+  steady_state_model.writeSteadyStateFile(basename, mod_file_struct.ramsey_model_present, true);
+
+  jlOutputFile << "model.static = " << basename << "Static.static!" << endl
+               << "model.dynamic = " << basename << "Dynamic.dynamic!" << endl
+               << "model.steady_state = " << basename << "SteadyState2.steady_state!" << endl
+               << "try" << endl
+               << "    using " << basename << "StaticParamsDerivs" << endl
+               << "    model.static_params_derivs = " << basename
+               << "StaticParamsDerivs.params_derivs" << endl
+               << "catch" << endl
+               << "end" << endl
+               << "try" << endl
+               << "    using " << basename << "DynamicParamsDerivs" << endl
+               << "    model.dynamic_params_derivs = " << basename
+               << "DynamicParamsDerivs.params_derivs" << endl
+               << "catch" << endl
+               << "end" << endl
+               << "end" << endl;
+  jlOutputFile.close();
+  cout << "done" << endl;
+}
diff --git a/preprocessor/ModFile.hh b/preprocessor/ModFile.hh
index 23ea3371e530a1fa61491a5241e23b71e4c79176..b53cca5ec27872a41ccca2ca69d3acd94c3023ba 100644
--- a/preprocessor/ModFile.hh
+++ b/preprocessor/ModFile.hh
@@ -158,6 +158,7 @@ public:
   void writeExternalFiles(const string &basename, FileOutputType output, LanguageOutputType language) const;
   void writeExternalFilesC(const string &basename, FileOutputType output) const;
   void writeExternalFilesCC(const string &basename, FileOutputType output) const;
+  void writeExternalFilesJulia(const string &basename, FileOutputType output) const;
   //! Writes C output files only => No further Matlab processing
   void writeCOutputFiles(const string &basename) const;
   void writeModelC(const string &basename) const;
diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc
index 14149f59a3a5ef6bebe1354bbf79c51f80200819..92b3996386a1021d5c54e3ae5b45da417d707cb0 100644
--- a/preprocessor/ModelTree.cc
+++ b/preprocessor/ModelTree.cc
@@ -1127,6 +1127,8 @@ ModelTree::writeTemporaryTerms(const temporary_terms_t &tt, ostream &output,
 
       if (IS_C(output_type))
         output << "double ";
+      else if (IS_JULIA(output_type))
+        output << "  @inbounds ";
 
       (*it)->writeOutput(output, output_type, tt, tef_terms);
       output << " = ";
@@ -1199,6 +1201,8 @@ ModelTree::writeModelLocalVariables(ostream &output, ExprNodeOutputType output_t
 
       if (IS_C(output_type))
         output << "double ";
+      else if (IS_JULIA(output_type))
+        output << "  @inbounds ";
 
       /* We append underscores to avoid name clashes with "g1" or "oo_" (see
          also VariableNode::writeOutput) */
@@ -1229,14 +1233,20 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type)
 
       if (vrhs != 0) // The right hand side of the equation is not empty ==> residual=lhs-rhs;
         {
+          if (IS_JULIA(output_type))
+            output << "  @inbounds ";
           output << "lhs =";
           lhs->writeOutput(output, output_type, temporary_terms);
           output << ";" << endl;
 
+          if (IS_JULIA(output_type))
+            output << "  @inbounds ";
           output << "rhs =";
           rhs->writeOutput(output, output_type, temporary_terms);
           output << ";" << endl;
 
+          if (IS_JULIA(output_type))
+            output << "  @inbounds ";
           output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
                  << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
                  << RIGHT_ARRAY_SUBSCRIPT(output_type)
@@ -1244,6 +1254,8 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type)
         }
       else // The right hand side of the equation is empty ==> residual=lhs;
         {
+          if (IS_JULIA(output_type))
+            output << "  @inbounds ";
           output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type)
                  << eq + ARRAY_SUBSCRIPT_OFFSET(output_type)
                  << RIGHT_ARRAY_SUBSCRIPT(output_type)
@@ -1470,8 +1482,11 @@ ModelTree::set_cutoff_to_zero()
 void
 ModelTree::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const
 {
-  output << "  g1" << LEFT_ARRAY_SUBSCRIPT(output_type);
-  if (IS_MATLAB(output_type))
+  output << "  ";
+  if (IS_JULIA(output_type))
+    output << "@inbounds ";
+  output << "g1" << LEFT_ARRAY_SUBSCRIPT(output_type);
+  if (IS_MATLAB(output_type) || IS_JULIA(output_type))
     output << eq_nb + 1 << "," << col_nb + 1;
   else
     output << eq_nb + col_nb *equations.size();
@@ -1482,7 +1497,7 @@ void
 ModelTree::sparseHelper(int order, ostream &output, int row_nb, int col_nb, ExprNodeOutputType output_type) const
 {
   output << "  v" << order << LEFT_ARRAY_SUBSCRIPT(output_type);
-  if (IS_MATLAB(output_type))
+  if (IS_MATLAB(output_type) || IS_JULIA(output_type))
     output << row_nb + 1 << "," << col_nb + 1;
   else
     output << row_nb + col_nb * NNZDerivatives[order-1];
diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc
index 523c0ab93dadeb32037ebad04c84b75bb6cf07b0..4fea0335eab4d76a2d3e6028743b1fad65ef2131 100644
--- a/preprocessor/StaticModel.cc
+++ b/preprocessor/StaticModel.cc
@@ -1175,26 +1175,29 @@ StaticModel::writeStaticMFile(const string &func_name) const
          << "% Warning : this file is generated automatically by Dynare" << endl
          << "%           from model file (.mod)" << endl << endl;
 
-  writeStaticModel(output, false);
+  writeStaticModel(output, false, false);
   output << "end" << endl;
   output.close();
 }
 
 void
-StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const
+StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) const
 {
-  ostringstream model_output;    // Used for storing model equations
+  ostringstream model_output;    // Used for storing model
+  ostringstream model_eq_output; // Used for storing model equations
   ostringstream jacobian_output; // Used for storing jacobian equations
   ostringstream hessian_output;  // Used for storing Hessian equations
   ostringstream third_derivatives_output;  // Used for storing third order derivatives equations
-  ExprNodeOutputType output_type = (use_dll ? oCStaticModel : oMatlabStaticModel);
+  ostringstream for_sym;
+  ExprNodeOutputType output_type = (use_dll ? oCStaticModel :
+                                    julia ? oJuliaStaticModel : oMatlabStaticModel);
 
   deriv_node_temp_terms_t tef_terms;
   writeModelLocalVariables(model_output, output_type, tef_terms);
 
   writeTemporaryTerms(temporary_terms, model_output, output_type, tef_terms);
 
-  writeModelEquations(model_output, output_type);
+  writeModelEquations(model_eq_output, output_type);
 
   int nrows = equations.size();
   int JacobianColsNbr = symbol_table.endo_nbr();
@@ -1231,35 +1234,49 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const
       int col_nb = tsid1*symbol_table.endo_nbr()+tsid2;
       int col_nb_sym = tsid2*symbol_table.endo_nbr()+tsid1;
 
-      sparseHelper(2, hessian_output, k, 0, output_type);
-      hessian_output << "=" << eq + 1 << ";" << endl;
-
-      sparseHelper(2, hessian_output, k, 1, output_type);
-      hessian_output << "=" << col_nb + 1 << ";" << endl;
-
-      sparseHelper(2, hessian_output, k, 2, output_type);
-      hessian_output << "=";
-      d2->writeOutput(hessian_output, output_type, temporary_terms, tef_terms);
-      hessian_output << ";" << endl;
-
-      k++;
-
-      // Treating symetric elements
-      if (symb_id1 != symb_id2)
+      if (output_type == oJuliaDynamicModel)
+        {
+          for_sym << "g2[" << eq + 1 << "," << col_nb + 1 << "]";
+          hessian_output << "  @inbounds " << for_sym.str() << " = ";
+          d2->writeOutput(hessian_output, output_type, temporary_terms, tef_terms);
+          hessian_output << endl;
+        }
+      else
         {
           sparseHelper(2, hessian_output, k, 0, output_type);
           hessian_output << "=" << eq + 1 << ";" << endl;
 
           sparseHelper(2, hessian_output, k, 1, output_type);
-          hessian_output << "=" << col_nb_sym + 1 << ";" << endl;
+          hessian_output << "=" << col_nb + 1 << ";" << endl;
 
           sparseHelper(2, hessian_output, k, 2, output_type);
           hessian_output << "=";
-          sparseHelper(2, hessian_output, k-1, 2, output_type);
+          d2->writeOutput(hessian_output, output_type, temporary_terms, tef_terms);
           hessian_output << ";" << endl;
 
           k++;
         }
+
+      // Treating symetric elements
+      if (symb_id1 != symb_id2)
+        if (output_type == oJuliaDynamicModel)
+          hessian_output << "  @inbounds g2[" << eq + 1 << "," << col_nb_sym + 1 << "] = "
+                         << for_sym.str() << endl;
+        else
+          {
+            sparseHelper(2, hessian_output, k, 0, output_type);
+            hessian_output << "=" << eq + 1 << ";" << endl;
+
+            sparseHelper(2, hessian_output, k, 1, output_type);
+            hessian_output << "=" << col_nb_sym + 1 << ";" << endl;
+
+            sparseHelper(2, hessian_output, k, 2, output_type);
+            hessian_output << "=";
+            sparseHelper(2, hessian_output, k-1, 2, output_type);
+            hessian_output << ";" << endl;
+
+            k++;
+          }
     }
 
   // Writing third derivatives
@@ -1281,18 +1298,29 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const
       // Reference column number for the g3 matrix
       int ref_col = id1 * hessianColsNbr + id2 * JacobianColsNbr + id3;
 
-      sparseHelper(3, third_derivatives_output, k, 0, output_type);
-      third_derivatives_output << "=" << eq + 1 << ";" << endl;
+      if (output_type == oJuliaDynamicModel)
+        {
+          for_sym << "g3[" << eq + 1 << "," << ref_col + 1 << "]";
+          third_derivatives_output << "  @inbounds " << for_sym.str() << " = ";
+          d3->writeOutput(third_derivatives_output, output_type, temporary_terms, tef_terms);
+          third_derivatives_output << endl;
+        }
+      else
+        {
+          sparseHelper(3, third_derivatives_output, k, 0, output_type);
+          third_derivatives_output << "=" << eq + 1 << ";" << endl;
 
-      sparseHelper(3, third_derivatives_output, k, 1, output_type);
-      third_derivatives_output << "=" << ref_col + 1 << ";" << endl;
+          sparseHelper(3, third_derivatives_output, k, 1, output_type);
+          third_derivatives_output << "=" << ref_col + 1 << ";" << endl;
 
-      sparseHelper(3, third_derivatives_output, k, 2, output_type);
-      third_derivatives_output << "=";
-      d3->writeOutput(third_derivatives_output, output_type, temporary_terms, tef_terms);
-      third_derivatives_output << ";" << endl;
+          sparseHelper(3, third_derivatives_output, k, 2, output_type);
+          third_derivatives_output << "=";
+          d3->writeOutput(third_derivatives_output, output_type, temporary_terms, tef_terms);
+          third_derivatives_output << ";" << endl;
+        }
 
-      // Compute the column numbers for the 5 other permutations of (id1,id2,id3) and store them in a set (to avoid duplicates if two indexes are equal)
+      // Compute the column numbers for the 5 other permutations of (id1,id2,id3)
+      // and store them in a set (to avoid duplicates if two indexes are equal)
       set<int> cols;
       cols.insert(id1 * hessianColsNbr + id3 * JacobianColsNbr + id2);
       cols.insert(id2 * hessianColsNbr + id1 * JacobianColsNbr + id3);
@@ -1303,30 +1331,35 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const
       int k2 = 1; // Keeps the offset of the permutation relative to k
       for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
         if (*it2 != ref_col)
-          {
-            sparseHelper(3, third_derivatives_output, k+k2, 0, output_type);
-            third_derivatives_output << "=" << eq + 1 << ";" << endl;
+          if (output_type == oJuliaDynamicModel)
+            third_derivatives_output << "  @inbounds g3[" << eq + 1 << "," << *it2 + 1 << "] = "
+                                     << for_sym.str() << endl;
+          else
+            {
+              sparseHelper(3, third_derivatives_output, k+k2, 0, output_type);
+              third_derivatives_output << "=" << eq + 1 << ";" << endl;
 
-            sparseHelper(3, third_derivatives_output, k+k2, 1, output_type);
-            third_derivatives_output << "=" << *it2 + 1 << ";" << endl;
+              sparseHelper(3, third_derivatives_output, k+k2, 1, output_type);
+              third_derivatives_output << "=" << *it2 + 1 << ";" << endl;
 
-            sparseHelper(3, third_derivatives_output, k+k2, 2, output_type);
-            third_derivatives_output << "=";
-            sparseHelper(3, third_derivatives_output, k, 2, output_type);
-            third_derivatives_output << ";" << endl;
+              sparseHelper(3, third_derivatives_output, k+k2, 2, output_type);
+              third_derivatives_output << "=";
+              sparseHelper(3, third_derivatives_output, k, 2, output_type);
+              third_derivatives_output << ";" << endl;
 
-            k2++;
-          }
+              k2++;
+            }
       k += k2;
     }
 
-  if (!use_dll)
+  if (output_type == oMatlabStaticModel)
     {
       StaticOutput << "residual = zeros( " << equations.size() << ", 1);" << endl << endl
                    << "%" << endl
                    << "% Model equations" << endl
                    << "%" << endl << endl
                    << model_output.str()
+                   << model_eq_output.str()
                    << "if ~isreal(residual)" << endl
                    << "  residual = real(residual)+imag(residual).^2;" << endl
                    << "end" << endl
@@ -1366,9 +1399,8 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const
                       << "  g3 = sparse(v3(:,1),v3(:,2),v3(:,3)," << nrows << "," << ncols << ");" << endl;
       else // Either 3rd derivatives is all zero, or we didn't compute it
         StaticOutput << "  g3 = sparse([],[],[]," << nrows << "," << ncols << ");" << endl;
-
     }
-  else
+  else if (output_type == oCStaticModel)
     {
       StaticOutput << "void Static(double *y, double *x, int nb_row_x, double *params, double *residual, double *g1, double *v2)" << endl
                    << "{" << endl
@@ -1376,6 +1408,7 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const
                    << endl
                    << "  /* Residual equations */" << endl
                    << model_output.str()
+                   << model_eq_output.str()
                    << "  /* Jacobian  */" << endl
                    << "  if (g1 == NULL)" << endl
                    << "    return;" << endl
@@ -1401,6 +1434,103 @@ StaticModel::writeStaticModel(ostream &StaticOutput, bool use_dll) const
                       << third_derivatives_output.str()
                       << "    }" << endl;
     }
+  else
+    {
+      ostringstream comments;
+      comments << "## Function Arguments" << endl
+               << endl
+               << "## Input" << endl
+               << " 1 y:         Array{Float64, length(model.endo), 1}             Vector of endogenous variables in declaration order" << endl
+               << " 2 x:         Array{Float64, length(model.exo), 1}              Vector of exogenous variables in declaration order" << endl
+               << " 3 params:    Array{Float64, length(model.param), 1}            Vector of parameter values in declaration order" << endl
+               << endl
+               << "## Output" << endl
+               << " 4 residual:  Array(Float64, model.eq_nbr, 1)                   Vector of residuals of the static model equations" << endl
+               << "                                                                in order of declaration of the equations." << endl
+               << "                                                                Dynare may prepend auxiliary equations, see model.aux_vars" << endl;
+
+      StaticOutput << "function static!(y::Vector{Float64}, x::Vector{Float64}, "
+                   << "params::Vector{Float64}," << endl
+                   << "                 residual::Vector{Float64})" << endl
+                   << "#=" << endl << comments.str() << "=#" << endl
+                   << "  @assert size(y) == " << symbol_table.endo_nbr() << endl
+                   << "  @assert size(x) == " << symbol_table.exo_nbr() << endl
+                   << "  @assert size(params) == " << symbol_table.param_nbr() << endl
+                   << "  @assert size(residual) == " << equations.size() << endl
+                   << "  #" << endl
+                   << "  # Model equations" << endl
+                   << "  #" << endl
+                   << model_output.str()
+                   << model_eq_output.str()
+                   << "if ~isreal(residual)" << endl
+                   << "  residual = real(residual)+imag(residual).^2;" << endl
+                   << "end" << endl
+                   << "end" << endl << endl
+                   << "function static!(y::Vector{Float64}, x::Vector{Float64}, "
+                   << "params::Vector{Float64}," << endl
+                   << "                 residual::Vector{Float64}, g1::Matrix{Float64})" << endl;
+
+      comments << " 5 g1:        Array(Float64, model.eq_nbr, length(model.endo))  Jacobian matrix of the static model equations;" << endl
+               << "                                                                columns: variables in declaration order" << endl
+               << "                                                                rows: equations in order of declaration" << endl;
+
+      StaticOutput << "#=" << endl << comments.str() << "=#" << endl
+                   << "  @assert size(g1) == (" << equations.size() << ", " << symbol_table.endo_nbr()
+                   << ")" << endl
+                   << "  fill!(g1, 0.0)" << endl
+                   << "  static!(y, x, params, residual)" << endl
+                   << model_output.str()
+                   << "  #" << endl
+                   << "  # Jacobian matrix" << endl
+                   << "  #" << endl
+                   << jacobian_output.str()
+                   << "  if ~isreal(g1)" << endl
+                   << "    g1 = real(g1)+2*imag(g1);" << endl
+                   << "  end" << endl
+                   << "end" << endl << endl
+                   << "function static!(y::Vector{Float64}, x::Vector{Float64}, "
+                   << "params::Vector{Float64}," << endl
+                   << "                 residual::Vector{Float64}, g1::Matrix{Float64}, "
+                   << "g2::Matrix{Float64})" << endl;
+
+      comments << " 6 g2:        spzeros(model.eq_nbr, length(model.endo)^2)       Hessian matrix of the static model equations;" << endl
+               << "                                                                columns: variables in declaration order" << endl
+               << "                                                                rows: equations in order of declaration" << endl;
+
+      StaticOutput << "#=" << endl << comments.str() << "=#" << endl
+                   << "  @assert size(g2) == (" << equations.size() << ", " << g2ncols << ")" << endl
+                   << "  static!(y, x, params, residual, g1)" << endl;
+      if (second_derivatives.size())
+        StaticOutput << model_output.str()
+                     << "  #" << endl
+                     << "  # Hessian matrix" << endl
+                     << "  #" << endl
+                     << hessian_output.str();
+
+      // Initialize g3 matrix
+      int ncols = hessianColsNbr * JacobianColsNbr;
+      StaticOutput << "end" << endl << endl
+                   << "function static!(y::Vector{Float64}, x::Vector{Float64}, "
+                   << "params::Vector{Float64}," << endl
+                   << "                 residual::Vector{Float64}, g1::Matrix{Float64}, "
+                   << "g2::Matrix{Float64}," << endl
+                   << "                 g3::Matrix{Float64})" << endl;
+
+      comments << " 7 g3:        spzeros(model.eq_nbr, length(model.endo)^3)       Third derivatives matrix of the static model equations;" << endl
+               << "                                                                columns: variables in declaration order" << endl
+               << "                                                                rows: equations in order of declaration" << endl;
+
+      StaticOutput << "#=" << endl << comments.str() << "=#" << endl
+                   << "  @assert size(g3) == (" << nrows << ", " << ncols << ")" << endl
+                   << "  static!(y, x, params, residual, g1, g2)" << endl;
+      if (third_derivatives.size())
+        StaticOutput << model_output.str()
+                     << "  #" << endl
+                     << "  # Third order derivatives" << endl
+                     << "  #" << endl
+                     << third_derivatives_output.str();
+      StaticOutput << "end" << endl;
+    }
 }
 
 void
@@ -1440,7 +1570,7 @@ StaticModel::writeStaticCFile(const string &func_name) const
   writePowerDerivCHeader(output);
 
   // Writing the function body
-  writeStaticModel(output, true);
+  writeStaticModel(output, true, false);
   output << "}" << endl << endl;
 
   writePowerDeriv(output, true);
@@ -1515,7 +1645,30 @@ StaticModel::writeStaticCFile(const string &func_name) const
 }
 
 void
-StaticModel::writeStaticFile(const string &basename, bool block, bool bytecode, bool use_dll) const
+StaticModel::writeStaticJuliaFile(const string &basename) const
+{
+  string filename = basename + "Static.jl";
+  ofstream output;
+  output.open(filename.c_str(), ios::out | ios::binary);
+  if (!output.is_open())
+    {
+      cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  output << "module " << basename << "Static" << endl
+         << "#" << endl
+         << "# NB: this file was automatically generated by Dynare" << endl
+         << "#     from " << basename << ".mod" << endl
+         << "#" << endl
+         << "using Utils" << endl << endl
+         << "export static!" << endl << endl;
+  writeStaticModel(output, false, true);
+  output << "end" << endl;
+}
+
+void
+StaticModel::writeStaticFile(const string &basename, bool block, bool bytecode, bool use_dll, bool julia) const
 {
   int r;
 
@@ -1544,9 +1697,11 @@ StaticModel::writeStaticFile(const string &basename, bool block, bool bytecode,
     }
   else if(use_dll)
     writeStaticCFile(basename);
+  else if (julia)
+    writeStaticJuliaFile(basename);
   else
     writeStaticMFile(basename);
-  writeAuxVarRecursiveDefinitions(basename);
+  writeAuxVarRecursiveDefinitions(basename, julia);
 }
 
 void
@@ -1906,10 +2061,11 @@ StaticModel::writeAuxVarInitval(ostream &output, ExprNodeOutputType output_type)
     }
 }
 
-void StaticModel::writeAuxVarRecursiveDefinitions(const string &basename) const
+void StaticModel::writeAuxVarRecursiveDefinitions(const string &basename, const bool julia) const
 {
   string func_name = basename + "_set_auxiliary_variables";
-  string filename = func_name + ".m";
+  string filename = julia ? func_name + ".jl" : func_name + ".m";
+  string comment = julia ? "#" : "%";
 
   ofstream output;
   output.open(filename.c_str(), ios::out | ios::binary);
@@ -1920,11 +2076,11 @@ void StaticModel::writeAuxVarRecursiveDefinitions(const string &basename) const
     }
 
   output << "function y = " << func_name + "(y, x, params)" << endl
-         << "%" << endl
-         << "% Status : Computes static model for Dynare" << endl
-         << "%" << endl
-         << "% Warning : this file is generated automatically by Dynare" << endl
-         << "%           from model file (.mod)" << endl
+         << comment << endl
+         << comment << " Status : Computes static model for Dynare" << endl
+         << comment << endl
+         << comment << " Warning : this file is generated automatically by Dynare" << endl
+         << comment << "           from model file (.mod)" << endl
          << endl;
 
   deriv_node_temp_terms_t tef_terms;
@@ -1942,7 +2098,7 @@ void StaticModel::writeAuxVarRecursiveDefinitions(const string &basename) const
 }
 
 void
-StaticModel::writeParamsDerivativesFile(const string &basename) const
+StaticModel::writeParamsDerivativesFile(const string &basename, bool julia) const
 {
   if (!residuals_params_derivatives.size()
       && !residuals_params_second_derivatives.size()
@@ -1951,8 +2107,7 @@ StaticModel::writeParamsDerivativesFile(const string &basename) const
       && !hessian_params_derivatives.size())
     return;
 
-  string filename = basename + "_static_params_derivs.m";
-
+  string filename = julia ? basename + "StaticParamsDerivs.jl" : basename + "_static_params_derivs.m";
   ofstream paramsDerivsFile;
   paramsDerivsFile.open(filename.c_str(), ios::out | ios::binary);
   if (!paramsDerivsFile.is_open())
@@ -1960,15 +2115,27 @@ StaticModel::writeParamsDerivativesFile(const string &basename) const
       cerr << "ERROR: Can't open file " << filename << " for writing" << endl;
       exit(EXIT_FAILURE);
     }
-  paramsDerivsFile << "function [rp, gp, rpp, gpp, hp] = " << basename << "_static_params_derivs(y, x, params)" << endl
-                   << "%" << endl
-                   << "% Warning : this file is generated automatically by Dynare" << endl
-                   << "%           from model file (.mod)" << endl << endl;
+
+  ExprNodeOutputType output_type = (julia ? oJuliaStaticModel : oMatlabStaticModel);
+
+  if (!julia)
+    paramsDerivsFile << "function [rp, gp, rpp, gpp, hp] = " << basename << "_static_params_derivs(y, x, params)" << endl
+                     << "%" << endl
+                     << "% Warning : this file is generated automatically by Dynare" << endl
+                     << "%           from model file (.mod)" << endl << endl;
+  else
+    paramsDerivsFile << "module " << basename << "StaticParamsDerivs" << endl
+                     << "#" << endl
+                     << "# NB: this file was automatically generated by Dynare" << endl
+                     << "#     from " << basename << ".mod" << endl
+                     << "#" << endl
+                     << "export params_derivs" << endl << endl
+                     << "function params_derivs(y, x, params)" << endl;
 
   deriv_node_temp_terms_t tef_terms;
-  writeModelLocalVariables(paramsDerivsFile, oMatlabStaticModel, tef_terms);
+  writeModelLocalVariables(paramsDerivsFile, output_type, tef_terms);
 
-  writeTemporaryTerms(params_derivs_temporary_terms, paramsDerivsFile, oMatlabStaticModel, tef_terms);
+  writeTemporaryTerms(params_derivs_temporary_terms, paramsDerivsFile, output_type, tef_terms);
 
   // Write parameter derivative
   paramsDerivsFile << "rp = zeros(" << equation_number() << ", "
@@ -1983,8 +2150,10 @@ StaticModel::writeParamsDerivativesFile(const string &basename) const
 
       int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
 
-      paramsDerivsFile << "rp(" << eq+1 << ", " << param_col << ") = ";
-      d1->writeOutput(paramsDerivsFile, oMatlabStaticModel, params_derivs_temporary_terms, tef_terms);
+      paramsDerivsFile << "rp" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                       <<  eq+1 << ", " << param_col
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << " = ";
+      d1->writeOutput(paramsDerivsFile, output_type, params_derivs_temporary_terms, tef_terms);
       paramsDerivsFile << ";" << endl;
     }
 
@@ -2003,13 +2172,16 @@ StaticModel::writeParamsDerivativesFile(const string &basename) const
       int var_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var)) + 1;
       int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
 
-      paramsDerivsFile << "gp(" << eq+1 << ", " << var_col << ", " << param_col << ") = ";
-      d2->writeOutput(paramsDerivsFile, oMatlabStaticModel, params_derivs_temporary_terms, tef_terms);
+      paramsDerivsFile << "gp" << LEFT_ARRAY_SUBSCRIPT(output_type)
+                       << eq+1 << ", " << var_col << ", " << param_col
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << " = ";
+      d2->writeOutput(paramsDerivsFile, output_type, params_derivs_temporary_terms, tef_terms);
       paramsDerivsFile << ";" << endl;
     }
 
-  // If nargout >= 3...
-  paramsDerivsFile << "if nargout >= 3" << endl;
+  if (!julia)
+    // If nargout >= 3...
+    paramsDerivsFile << "if nargout >= 3" << endl;
 
   // Write parameter second derivatives (only if nargout >= 3)
   paramsDerivsFile << "rpp = zeros(" << residuals_params_second_derivatives.size()
@@ -2027,11 +2199,15 @@ StaticModel::writeParamsDerivativesFile(const string &basename) const
       int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
       int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1;
 
-      paramsDerivsFile << "rpp(" << i << ",1)=" << eq+1 << ";" << endl
-                       << "rpp(" << i << ",2)=" << param1_col << ";" << endl
-                       << "rpp(" << i << ",3)=" << param2_col << ";" << endl
-                       << "rpp(" << i << ",4)=";
-      d2->writeOutput(paramsDerivsFile, oMatlabStaticModel, params_derivs_temporary_terms, tef_terms);
+      paramsDerivsFile << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
+                       << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
+                       << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
+                       << "rpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+      d2->writeOutput(paramsDerivsFile, output_type, params_derivs_temporary_terms, tef_terms);
       paramsDerivsFile << ";" << endl;
     }
 
@@ -2053,18 +2229,24 @@ StaticModel::writeParamsDerivativesFile(const string &basename) const
       int param1_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param1)) + 1;
       int param2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param2)) + 1;
 
-      paramsDerivsFile << "gpp(" << i << ",1)=" << eq+1 << ";" << endl
-                       << "gpp(" << i << ",2)=" << var_col << ";" << endl
-                       << "gpp(" << i << ",3)=" << param1_col << ";" << endl
-                       << "gpp(" << i << ",4)=" << param2_col << ";" << endl
-                       << "gpp(" << i << ",5)=";
-      d2->writeOutput(paramsDerivsFile, oMatlabStaticModel, params_derivs_temporary_terms, tef_terms);
+      paramsDerivsFile << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
+                       << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var_col << ";" << endl
+                       << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param1_col << ";" << endl
+                       << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param2_col << ";" << endl
+                       << "gpp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+      d2->writeOutput(paramsDerivsFile, output_type, params_derivs_temporary_terms, tef_terms);
       paramsDerivsFile << ";" << endl;
     }
 
-  // If nargout >= 5...
-  paramsDerivsFile << "end" << endl
-                   << "if nargout >= 5" << endl;
+  if (!julia)
+    // If nargout >= 5...
+    paramsDerivsFile << "end" << endl
+                     << "if nargout >= 5" << endl;
 
   // Write hessian derivatives (only if nargout >= 5)
   paramsDerivsFile << "hp = zeros(" << hessian_params_derivatives.size() << ",5);" << endl;
@@ -2083,15 +2265,22 @@ StaticModel::writeParamsDerivativesFile(const string &basename) const
       int var2_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(var2)) + 1;
       int param_col = symbol_table.getTypeSpecificID(getSymbIDByDerivID(param)) + 1;
 
-      paramsDerivsFile << "hp(" << i << ",1)=" << eq+1 << ";" << endl
-                       << "hp(" << i << ",2)=" << var1_col << ";" << endl
-                       << "hp(" << i << ",3)=" << var2_col << ";" << endl
-                       << "hp(" << i << ",4)=" << param_col << ";" << endl
-                       << "hp(" << i << ",5)=";
-      d2->writeOutput(paramsDerivsFile, oMatlabStaticModel, params_derivs_temporary_terms, tef_terms);
+      paramsDerivsFile << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",1"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << eq+1 << ";" << endl
+                       << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",2"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var1_col << ";" << endl
+                       << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",3"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << var2_col << ";" << endl
+                       << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",4"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=" << param_col << ";" << endl
+                       << "hp" << LEFT_ARRAY_SUBSCRIPT(output_type) << i << ",5"
+                       << RIGHT_ARRAY_SUBSCRIPT(output_type) << "=";
+      d2->writeOutput(paramsDerivsFile, output_type, params_derivs_temporary_terms, tef_terms);
       paramsDerivsFile << ";" << endl;
     }
 
+  if (julia)
+    paramsDerivsFile << "(rp, gp, rpp, gpp, hp)" << endl;
   paramsDerivsFile << "end" << endl
                    << "end" << endl;
   paramsDerivsFile.close();
diff --git a/preprocessor/StaticModel.hh b/preprocessor/StaticModel.hh
index e76a5f0fe3d1a0c1ac4a0b1c7ae9775c54c50f98..90a3aebb7962c27ec1b8952d8cddf8f17103e1f4 100644
--- a/preprocessor/StaticModel.hh
+++ b/preprocessor/StaticModel.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2012 Dynare Team
+ * Copyright (C) 2003-2015 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -50,8 +50,11 @@ private:
   //! Writes static model file (C version)
   void writeStaticCFile(const string &func_name) const;
 
+  //! Writes static model file (Julia version)
+  void writeStaticJuliaFile(const string &basename) const;
+
   //! Writes the static model equations and its derivatives
-  void writeStaticModel(ostream &StaticOutput, bool use_dll) const;
+  void writeStaticModel(ostream &StaticOutput, bool use_dll, bool julia) const;
 
   //! Writes the static function calling the block to solve (Matlab version)
   void writeStaticBlockMFSFile(const string &basename) const;
@@ -168,10 +171,10 @@ public:
                                    int &u_count_int, bool &file_open) const;
 
   //! Writes static model file
-  void writeStaticFile(const string &basename, bool block, bool bytecode, bool use_dll) const;
+  void writeStaticFile(const string &basename, bool block, bool bytecode, bool use_dll, bool julia) const;
 
   //! Writes file containing static parameters derivatives
-  void writeParamsDerivativesFile(const string &basename) const;
+  void writeParamsDerivativesFile(const string &basename, bool julia) const;
 
   //! Writes LaTeX file with the equations of the static model
   void writeLatexFile(const string &basename) const;
@@ -179,8 +182,8 @@ public:
   //! Writes initializations in oo_.steady_state or steady state file for the auxiliary variables
   void writeAuxVarInitval(ostream &output, ExprNodeOutputType output_type) const;
 
-  //! Writes definition of the auxiliary variables in a M file
-  void writeAuxVarRecursiveDefinitions(const string &basename) const;
+  //! Writes definition of the auxiliary variables in a .m or .jl file
+  void writeAuxVarRecursiveDefinitions(const string &basename, const bool julia) const;
 
   virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException);
   virtual void addAllParamDerivId(set<int> &deriv_id_set);
diff --git a/preprocessor/SteadyStateModel.cc b/preprocessor/SteadyStateModel.cc
index fe1a773b6f90af29cf4e25b6a92e0824085b10e8..c94590cf1be6a017efcb56648c8f37bf84f916cd 100644
--- a/preprocessor/SteadyStateModel.cc
+++ b/preprocessor/SteadyStateModel.cc
@@ -105,13 +105,12 @@ SteadyStateModel::checkPass(bool ramsey_model, WarningConsolidation &warnings) c
 }
 
 void
-SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_model) const
+SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_model, bool julia) const
 {
   if (def_table.size() == 0)
     return;
 
-  string filename = basename + "_steadystate2.m";
-
+  string filename = julia ? basename + "SteadyState2.jl" : basename + "_steadystate2.m";
   ofstream output;
   output.open(filename.c_str(), ios::out | ios::binary);
   if (!output.is_open())
@@ -120,10 +119,22 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_model
       exit(EXIT_FAILURE);
     }
 
-  output << "function [ys_, params, info] = " << basename << "_steadystate2("
-	 << "ys_, exo_, params)" << endl
-         << "% Steady state generated by Dynare preprocessor" << endl
-         << "    info = 0;" << endl;
+  ExprNodeOutputType output_type = (julia ? oJuliaSteadyStateFile : oSteadyStateFile);
+
+  if (!julia)
+    output << "function [ys_, params, info] = " << basename << "_steadystate2("
+           << "ys_, exo_, params)" << endl
+           << "% Steady state generated by Dynare preprocessor" << endl
+           << "    info = 0;" << endl;
+  else
+    output << "module " << basename << "SteadyState2" << endl
+           << "#" << endl
+           << "# NB: this file was automatically generated by Dynare" << endl
+           << "#     from " << basename << ".mod" << endl
+           << "#" << endl
+           << "export steady_state!" << endl << endl
+           << "function steady_state!(ys_::Vector{Float64}, exo_::Vector{Float64}, "
+           << "params::Vector{Float64})" << endl;
 
   for (size_t i = 0; i < def_table.size(); i++)
     {
@@ -135,7 +146,7 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_model
         {
           variable_node_map_t::const_iterator it = variable_node_map.find(make_pair(symb_ids[j], 0));
           assert(it != variable_node_map.end());
-          dynamic_cast<ExprNode *>(it->second)->writeOutput(output, oSteadyStateFile);
+          dynamic_cast<ExprNode *>(it->second)->writeOutput(output, output_type);
           if (j < symb_ids.size()-1)
             output << ",";
         }
@@ -143,13 +154,21 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_model
         output << "]";
 
       output << "=";
-      def_table[i].second->writeOutput(output, oSteadyStateFile);
+      def_table[i].second->writeOutput(output, output_type);
       output << ";" << endl;
     }
-  output << "    % Auxiliary equations" << endl;
-  static_model.writeAuxVarInitval(output, oSteadyStateFile);
-  output << "    check_=0;" << endl
-         << "end" << endl;
+  if (!julia)
+    output << "    % Auxiliary equations" << endl;
+  else
+    output << "    # Auxiliary equations" << endl;
+  static_model.writeAuxVarInitval(output, output_type);
+
+  if (!julia)
+    output << "    check_=0;" << endl;
+
+  output << "end" << endl;
+  if (julia)
+    output << "end" << endl;
 }
 
 void
diff --git a/preprocessor/SteadyStateModel.hh b/preprocessor/SteadyStateModel.hh
index 2dd752ba94be75820324daf6d4f0d3efdbce3543..f79b358f5c6bebd9c5c29259df6ff87eafd22bfe 100644
--- a/preprocessor/SteadyStateModel.hh
+++ b/preprocessor/SteadyStateModel.hh
@@ -48,7 +48,7 @@ public:
   /*!
     \param[in] ramsey_model Is there a Ramsey model in the MOD file? If yes, then use the "ys" in argument of the steady state file as initial values
   */
-  void writeSteadyStateFile(const string &basename, bool ramsey_model) const;
+  void writeSteadyStateFile(const string &basename, bool ramsey_model, bool julia) const;
   void writeSteadyStateFileC(const string &basename, bool ramsey_model) const;
 };
 
diff --git a/preprocessor/SymbolTable.cc b/preprocessor/SymbolTable.cc
index 332c94a46595f7e9ea91820faaa89107a00cc833..5599482ebf248d7d21dbe21461d4f7452abcf7fc 100644
--- a/preprocessor/SymbolTable.cc
+++ b/preprocessor/SymbolTable.cc
@@ -709,3 +709,116 @@ SymbolTable::getOrigEndogenous() const
       origendogs.insert(it->second);
   return origendogs;
 }
+
+void
+SymbolTable::writeJuliaOutput(ostream &output) const throw (NotYetFrozenException)
+{
+  if (!frozen)
+    throw NotYetFrozenException();
+
+  output << "# Endogenous Variables" << endl
+         << "model.endo = [" << endl;
+  if (endo_nbr() > 0)
+    for (int id = 0; id < endo_nbr(); id++)
+      output << "              DynareModel.Endo(\""
+             << getName(endo_ids[id]) << "\", \""
+             << getTeXName(endo_ids[id]) << "\", \""
+             << getLongName(endo_ids[id]) << "\")" << endl;
+  output << "             ]" << endl;
+
+  output << "# Exogenous Variables" << endl
+         << "model.exo = [" << endl;
+  if (exo_nbr() > 0)
+    for (int id = 0; id < exo_nbr(); id++)
+      output << "             DynareModel.Exo(\""
+             << getName(exo_ids[id]) << "\", \""
+             << getTeXName(exo_ids[id]) << "\", \""
+             << getLongName(exo_ids[id]) << "\")" << endl;
+  output << "            ]" << endl;
+
+  if (exo_det_nbr() > 0)
+    {
+      output << "# Exogenous Deterministic Variables" << endl
+             << "model.exo_det = [" << endl;
+      if (exo_det_nbr() > 0)
+        for (int id = 0; id < exo_det_nbr(); id++)
+          output << "                 DynareModel.ExoDet(\""
+                 << getName(exo_det_ids[id]) << "\", \""
+                 << getTeXName(exo_det_ids[id]) << "\", \""
+                 << getLongName(exo_det_ids[id]) << "\")" << endl;
+      output << "                ]" << endl;
+    }
+
+  output << "# Parameters" << endl
+         << "model.param = [" << endl;
+  if (param_nbr() > 0)
+    for (int id = 0; id < param_nbr(); id++)
+      output << "               DynareModel.Param(\""
+             << getName(param_ids[id]) << "\", \""
+             << getTeXName(param_ids[id]) << "\", \""
+             << getLongName(param_ids[id]) << "\")" << endl;
+  output << "              ]" << endl;
+
+  output << "model.orig_endo_nbr = " << orig_endo_nbr() << endl;
+
+  if (aux_vars.size() > 0)
+    {
+      output << "# Auxiliary Variables" << endl
+             << "model.aux_vars = [" << endl;
+      for (int i = 0; i < (int) aux_vars.size(); i++)
+        {
+          output << "                   DynareModel.AuxVars("
+                 << getTypeSpecificID(aux_vars[i].get_symb_id()) + 1 << ", "
+                 << aux_vars[i].get_type() << ", ";
+          switch (aux_vars[i].get_type())
+            {
+            case avEndoLead:
+            case avExoLead:
+              break;
+            case avEndoLag:
+            case avExoLag:
+              output << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) + 1 << ", "
+                     << aux_vars[i].get_orig_lead_lag() << ", NaN, NaN";
+              break;
+            case avMultiplier:
+              output << "NaN, NaN, " << aux_vars[i].get_equation_number_for_multiplier() + 1
+                     << ", NaN";
+              break;
+            case avDiffForward:
+              output << getTypeSpecificID(aux_vars[i].get_orig_symb_id())+1 << ", NaN, ";
+              break;
+            case avExpectation:
+              output << "NaN, NaN, NaN, \"\\mathbb{E}_{t"
+                     << (aux_vars[i].get_information_set() < 0 ? "" : "+")
+                     << aux_vars[i].get_information_set() << "}(";
+              aux_vars[i].get_expectation_expr_node()->writeOutput(output, oLatexDynamicModel);
+              output << ")\"";
+              break;
+            }
+          output << ")" << endl;
+        }
+      output << "]" << endl;
+    }
+
+    if (predeterminedNbr() > 0)
+      {
+        output << "# Predetermined Variables" << endl
+               << "model.pred_vars = [ " << endl;
+        for (set<int>::const_iterator it = predetermined_variables.begin();
+             it != predetermined_variables.end(); it++)
+          output << "                   DynareModel.PredVars("
+                 << getTypeSpecificID(*it)+1 << ")" << endl;
+        output << "                  ]" << endl;
+      }
+
+    if (observedVariablesNbr() > 0)
+      {
+        output << "# Observed Variables" << endl
+               << "options.obs_vars = [" << endl;
+        for (vector<int>::const_iterator it = varobs.begin();
+             it != varobs.end(); it++)
+          output << "                    DynareModel.ObsVars("
+                 << getTypeSpecificID(*it)+1 << ")" << endl;
+        output << "                   ]" << endl;
+      }
+}
diff --git a/preprocessor/SymbolTable.hh b/preprocessor/SymbolTable.hh
index 4be960a1c87bb6265a977ee2af25b0ff166801a3..b93e1b3fe312fab52cf337c4aabcfea69c31d16e 100644
--- a/preprocessor/SymbolTable.hh
+++ b/preprocessor/SymbolTable.hh
@@ -283,6 +283,8 @@ public:
   inline int orig_endo_nbr() const throw (NotYetFrozenException);
   //! Write output of this class
   void writeOutput(ostream &output) const throw (NotYetFrozenException);
+  //! Write Julia output of this class
+  void writeJuliaOutput(ostream &output) const throw (NotYetFrozenException);
   //! Write C output of this class
   void writeCOutput(ostream &output) const throw (NotYetFrozenException);
   //! Write CC output of this class
diff --git a/preprocessor/macro/MacroFlex.ll b/preprocessor/macro/MacroFlex.ll
index 9e694eb485b30ba0ae1ad2b1b1d7a66f7e1b54a7..9229b049dca2d5b8dff8e5d69302e7b05bee34f6 100644
--- a/preprocessor/macro/MacroFlex.ll
+++ b/preprocessor/macro/MacroFlex.ll
@@ -245,7 +245,7 @@ CONT \\\\
 <STMT>echo                  { return token::ECHO_DIR; }
 <STMT>error                 { return token::ERROR; }
 
-<STMT,EXPR>[A-Za-z_][A-Za-z0-9_]* {
+<STMT,EXPR>[A-Za-z_\x80-\xf3][A-Za-z0-9_\x80-\xf3]* {
                               yylval->string_val = new string(yytext);
                               return token::NAME;
                             }
diff --git a/tests/julia/rbc/README.md b/tests/julia/rbc/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..b556485bb8cb59aee45fe75a709131963647c2bb
--- /dev/null
+++ b/tests/julia/rbc/README.md
@@ -0,0 +1,8 @@
+To compile the Dynare file ```rbc.mod``` and produce a julia module, just do
+
+```
+include("test.jl")
+```
+
+in a script or in julia's shell.
+	
diff --git a/tests/julia/rbc/clean b/tests/julia/rbc/clean
new file mode 100755
index 0000000000000000000000000000000000000000..a48414dc7261d9d9258788923df024f5e26fa22d
--- /dev/null
+++ b/tests/julia/rbc/clean
@@ -0,0 +1,8 @@
+#/bin/sh
+rm -f *~
+rm -f rbc.jl
+rm -f rbcDynamic.jl
+rm -f rbcSteadyState2.jl
+rm -rf rbc
+rm -f rbcStatic.jl
+rm -f rbc_set_auxiliary_variables.jl
diff --git a/tests/julia/rbc/rbc.mod b/tests/julia/rbc/rbc.mod
new file mode 100644
index 0000000000000000000000000000000000000000..fcf6632cfc51531b9285c0c7d51e1ef9106db27c
--- /dev/null
+++ b/tests/julia/rbc/rbc.mod
@@ -0,0 +1,76 @@
+var Capital , Output, Labour, Consumption, Efficiency, efficiency ;
+
+varexo EfficiencyInnovation;
+
+parameters beta, theta, tau, alpha, Epsilon, delta, rho, effstar, sigma;
+
+beta    =  0.990;
+theta   =  0.357;
+tau     = 30.000;
+alpha   =  0.450;
+delta   =  0.020;
+rho     =  0.950;
+effstar =  1.500;
+sigma   =  0.010;
+Epsilon =  0.500;
+
+model;
+
+#Psi = (Epsilon-1)/Epsilon;
+
+  // Eq. n°1:
+  efficiency = rho*efficiency(-1) + sigma*EfficiencyInnovation;
+
+  // Eq. n°2:
+  Efficiency = effstar*exp(efficiency);
+
+  // Eq. n°3:
+  Output = Efficiency*(alpha*Capital(-1)^Psi+(1-alpha)*Labour^Psi)^(1/Psi);
+
+  // Eq. n°4:
+  Consumption + Capital - Output - (1-delta)*Capital(-1);
+
+  // Eq. n°5:
+  ((1-theta)/theta)*(Consumption/(1-Labour)) - (1-alpha)*Efficiency^((1-Psi))*(alpha*(Capital(-1)/Labour)^Psi+1-alpha)^((1-Psi)/Psi);
+
+  // Eq. n°6:
+  (((Consumption^theta)*((1-Labour)^(1-theta)))^(1-tau))/Consumption
+  - beta*(Consumption(1)^theta*(1-Labour(1))^(1-theta))^(1-tau)/Consumption(1)*(alpha*Efficiency(1)^Psi*(Output(1)/Capital)^(1-Psi)+1-delta);
+
+end;
+
+
+
+steady_state_model;
+
+  efficiency = 0;
+  Efficiency = effstar;
+
+  psi = (Epsilon-1)/Epsilon;
+
+  Output_per_unit_of_Capital = ( (1/beta-1+delta) / (alpha*effstar^psi) )^(1/(1-psi));
+
+  Consumption_per_unit_of_Capital = Output_per_unit_of_Capital-delta;
+
+  Labour_per_unit_of_Capital =  ((Output_per_unit_of_Capital/Efficiency)^psi-alpha)^(1/psi)/(1-alpha)^(1/psi);
+
+  gamma_1 = theta*(1-alpha)/(1-theta)*(Output_per_unit_of_Capital/Labour_per_unit_of_Capital)^(1-psi);
+  gamma_2 = (Output_per_unit_of_Capital-delta)/Labour_per_unit_of_Capital;
+
+  Labour = 1/(1+gamma_2/gamma_1);
+
+  Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital;
+
+  Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital;
+
+  ShareOfCapital= alpha^(1/(1-psi))*effstar^psi/(1/beta-1+delta)^(psi/(1-psi));
+  
+  Consumption = Consumption_per_unit_of_Labour*Labour;
+  Capital = Labour/Labour_per_unit_of_Capital;
+  Output = Output_per_unit_of_Capital*Capital;
+
+end;
+
+shocks;
+var EfficiencyInnovation = 1;
+end;
diff --git a/tests/julia/rbc/test.jl b/tests/julia/rbc/test.jl
new file mode 100644
index 0000000000000000000000000000000000000000..339e6f44d49bfda040f98b7893ebfa71f51b1668
--- /dev/null
+++ b/tests/julia/rbc/test.jl
@@ -0,0 +1,13 @@
+#clear workspace
+workspace()
+
+# Modification of the path (for packages). Should be done in ~/.juliarc.jl with a fixed path instead.
+if isempty(findin([abspath("../../../julia")], LOAD_PATH))
+    unshift!(LOAD_PATH, abspath("../../../julia"))
+end
+
+# Load Dynare package
+importall Dynare
+
+# Compile the rbc.mod file -> produce a module with the model definition.
+@dynare "rbc.mod"