diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index ef0f9f064efa313613c6d1323aedfcc8873cce83..478e9df8d818dd24443624276df0f5c89ac18553 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -4595,7 +4595,8 @@ void
 DynamicModel::writeSetAuxiliaryVariables(const string &basename, bool julia) const
 {
   ostringstream output_func_body;
-  writeAuxVarRecursiveDefinitions(output_func_body, ExprNodeOutputType::matlabDseries);
+  writeAuxVarRecursiveDefinitions(output_func_body, julia ? ExprNodeOutputType::juliaTimeDataFrame
+                                  : ExprNodeOutputType::matlabDseries);
 
   if (output_func_body.str().empty())
     return;
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index c6f3997edc1ed3c6a5bce210319b85d502745526..ae80658fc99a9fbbe763db1b29460d725d11de5e 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -948,6 +948,19 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       return;
     }
 
+  auto juliaTimeDataFrameHelper = [&]()
+  {
+    if (lag != 0)
+      output << "lag(";
+    output << "ds." << datatree.symbol_table.getName(symb_id);
+    if (lag != 0)
+      {
+        if (lag != -1)
+          output << "," << -lag;
+        output << ")";
+      }
+  };
+
   int i;
   switch (type)
     {
@@ -1012,6 +1025,9 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
           if (lag != 0)
             output << LEFT_ARRAY_SUBSCRIPT(output_type) << lag << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
+        case ExprNodeOutputType::juliaTimeDataFrame:
+          juliaTimeDataFrameHelper();
+          break;
         case ExprNodeOutputType::epilogueFile:
           output << "ds." << datatree.symbol_table.getName(symb_id);
           output << LEFT_ARRAY_SUBSCRIPT(output_type) << "t";
@@ -1073,6 +1089,9 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
           if (lag != 0)
             output << LEFT_ARRAY_SUBSCRIPT(output_type) << lag << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
+        case ExprNodeOutputType::juliaTimeDataFrame:
+          juliaTimeDataFrameHelper();
+          break;
         case ExprNodeOutputType::epilogueFile:
           output << "ds." << datatree.symbol_table.getName(symb_id);
           output << LEFT_ARRAY_SUBSCRIPT(output_type) << "t";
@@ -1131,6 +1150,9 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
           if (lag != 0)
             output << LEFT_ARRAY_SUBSCRIPT(output_type) << lag << RIGHT_ARRAY_SUBSCRIPT(output_type);
           break;
+        case ExprNodeOutputType::juliaTimeDataFrame:
+          juliaTimeDataFrameHelper();
+          break;
         case ExprNodeOutputType::epilogueFile:
           output << "ds." << datatree.symbol_table.getName(symb_id);
           output << LEFT_ARRAY_SUBSCRIPT(output_type) << "t";
@@ -1152,7 +1174,8 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
             output << lag;
           output << RIGHT_ARRAY_SUBSCRIPT(output_type);
         }
-      else if (output_type == ExprNodeOutputType::matlabDseries)
+      else if (output_type == ExprNodeOutputType::matlabDseries
+               || output_type == ExprNodeOutputType::juliaTimeDataFrame)
         // Only writing dseries for epilogue_static, hence no need to check lag
         output << "ds." << datatree.symbol_table.getName(symb_id);
       else
@@ -2842,6 +2865,10 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       break;
     }
 
+  if (output_type == ExprNodeOutputType::juliaTimeDataFrame
+      && op_code != UnaryOpcode::uminus)
+    output << "."; // Use vectorized form of the function
+
   bool close_parenthesis = false;
 
   /* Enclose argument with parentheses if:
@@ -4569,60 +4596,81 @@ BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
     case BinaryOpcode::times:
       if (isLatexOutput(output_type))
         output << R"(\, )";
-      else if (output_type == ExprNodeOutputType::occbinDifferenceFile)
-        output << ".*"; // This file operates on vectors, see dynare#1826
+      else if (output_type == ExprNodeOutputType::occbinDifferenceFile // This file operates on vectors, see dynare#1826
+               || output_type == ExprNodeOutputType::juliaTimeDataFrame)
+        output << ".*";
       else
         output << "*";
       break;
     case BinaryOpcode::divide:
       if (!isLatexOutput(output_type))
         {
-          if (output_type == ExprNodeOutputType::occbinDifferenceFile)
+          if (output_type == ExprNodeOutputType::occbinDifferenceFile // This file operates on vectors, see dynare#1826
+               || output_type == ExprNodeOutputType::juliaTimeDataFrame)
             output << "./"; // This file operates on vectors, see dynare#1826
           else
             output << "/";
         }
       break;
     case BinaryOpcode::power:
-      if (output_type == ExprNodeOutputType::occbinDifferenceFile)
+      if (output_type == ExprNodeOutputType::occbinDifferenceFile // This file operates on vectors, see dynare#1826
+          || output_type == ExprNodeOutputType::juliaTimeDataFrame)
         output << ".^"; // This file operates on vectors, see dynare#1826
       else
         output << "^";
       break;
     case BinaryOpcode::less:
-      output << "<";
+      if (output_type == ExprNodeOutputType::juliaTimeDataFrame)
+        output << ".<";
+      else
+        output << "<";
       break;
     case BinaryOpcode::greater:
-      output << ">";
+      if (output_type == ExprNodeOutputType::juliaTimeDataFrame)
+        output << ".>";
+      else
+        output << ">";
       break;
     case BinaryOpcode::lessEqual:
       if (isLatexOutput(output_type))
         output << R"(\leq )";
+      else if (output_type == ExprNodeOutputType::juliaTimeDataFrame)
+        output << ".<=";
       else
         output << "<=";
       break;
     case BinaryOpcode::greaterEqual:
       if (isLatexOutput(output_type))
         output << R"(\geq )";
+      else if (output_type == ExprNodeOutputType::juliaTimeDataFrame)
+        output << ".>=";
       else
         output << ">=";
       break;
     case BinaryOpcode::equalEqual:
-      output << "==";
+      if (output_type == ExprNodeOutputType::juliaTimeDataFrame)
+        output << ".==";
+      else
+        output << "==";
       break;
     case BinaryOpcode::different:
       if (isMatlabOutput(output_type))
         output << "~=";
       else
         {
-          if (isCOutput(output_type) || isJuliaOutput(output_type))
+          if (output_type == ExprNodeOutputType::juliaTimeDataFrame)
+            output << ".!=";
+          else if (isCOutput(output_type) || isJuliaOutput(output_type))
             output << "!=";
           else
             output << R"(\neq )";
         }
       break;
     case BinaryOpcode::equal:
-      output << "=";
+      if (output_type == ExprNodeOutputType::juliaTimeDataFrame)
+        output << ".=";
+      else
+        output << "=";
       break;
     default:
       ;
@@ -6029,7 +6077,10 @@ TrinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       else if (isJuliaOutput(output_type))
 	{
 	  // Julia API is normcdf(mu, sigma, x) !
-          output << "normcdf(";
+          output << "normcdf";
+          if (output_type == ExprNodeOutputType::juliaTimeDataFrame)
+            output << ".";
+          output << "(";
           arg2->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ",";
           arg3->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
@@ -6075,7 +6126,10 @@ TrinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
         }
       else
         {
-          output << "normpdf(";
+          output << "normpdf";
+          if (output_type == ExprNodeOutputType::juliaTimeDataFrame)
+            output << ".";
+          output << "(";
           arg1->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
           output << ",";
           arg2->writeOutput(output, output_type, temporary_terms, temporary_terms_idxs, tef_terms);
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 35fdc45c10f923106e5c67f3feed3d64905f7843..9ea42414c961f72e01f687c381b3f13add621678 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -97,6 +97,7 @@ enum class ExprNodeOutputType
    steadyStateFile, //!< Matlab code, in the generated steady state file
    juliaSteadyStateFile, //!< Julia code, in the generated steady state file
    matlabDseries, //!< Matlab code for dseries
+   juliaTimeDataFrame, //!< Julia code for TimeDataFrame objects
    epilogueFile, //!< Matlab code, in the generated epilogue file
    occbinDifferenceFile //!< MATLAB, in the generated occbin_difference file
   };
@@ -120,7 +121,8 @@ isJuliaOutput(ExprNodeOutputType output_type)
   return output_type == ExprNodeOutputType::juliaStaticModel
     || output_type == ExprNodeOutputType::juliaDynamicModel
     || output_type == ExprNodeOutputType::juliaDynamicSteadyStateOperator
-    || output_type == ExprNodeOutputType::juliaSteadyStateFile;
+    || output_type == ExprNodeOutputType::juliaSteadyStateFile
+    || output_type == ExprNodeOutputType::juliaTimeDataFrame;
 }
 
 inline bool