diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc
index c4f9ca5ae5c69667c617b6178755e34ad02b909d..5980658be5a39055008bb2961ca992e896397abc 100644
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -3880,8 +3880,8 @@ BasicPriorStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsoli
     }
 
   auto it_stdev = options_list.num_options.find("stdev");
-  if ((it_stdev == options_list.num_options.end() && variance == NULL)
-      || (it_stdev != options_list.num_options.end() && variance != NULL))
+  if ((it_stdev == options_list.num_options.end() && variance == nullptr)
+      || (it_stdev != options_list.num_options.end() && variance != nullptr))
     {
       cerr << "ERROR: You must pass exactly one of stdev and variance to the prior statement." << endl;
       exit(EXIT_FAILURE);
@@ -3973,7 +3973,7 @@ BasicPriorStatement::writeJsonPriorOutput(ostream &output) const
          << ", \"subsample\": \"" << subsample_name << "\""
          << ", ";
   writeJsonShape(output);
-  if (variance != NULL)
+  if (variance != nullptr)
     {
       output << ", \"variance\": \"";
       variance->writeJsonOutput(output, {}, {});
diff --git a/src/ConfigFile.cc b/src/ConfigFile.cc
index cc2c8182fc6448b10e5bc6856e04a5823eb4802b..a70e8945855810babc8dec552a1d2009a37ee395 100644
--- a/src/ConfigFile.cc
+++ b/src/ConfigFile.cc
@@ -122,7 +122,7 @@ ConfigFile::getConfigFileInfo(const string &config_file)
           defaultConfigFile += "\\dynare.ini";
         }
 #else
-      if (getenv("HOME") == NULL)
+      if (getenv("HOME") == nullptr)
         {
           if (parallel || parallel_test)
             cerr << "ERROR: ";
diff --git a/src/DataTree.cc b/src/DataTree.cc
index a0e429f9959a2d0a40936fc0352933cf7038f485..7793e42005d9831a13ad0677fed9db8d02b30a01 100644
--- a/src/DataTree.cc
+++ b/src/DataTree.cc
@@ -95,7 +95,7 @@ DataTree::AddPlus(expr_t iArg1, expr_t iArg2)
     {
       // Simplify x+(-y) in x-y
       auto *uarg2 = dynamic_cast<UnaryOpNode *>(iArg2);
-      if (uarg2 != NULL && uarg2->get_op_code() == oUminus)
+      if (uarg2 != nullptr && uarg2->get_op_code() == oUminus)
         return AddMinus(iArg1, uarg2->get_arg());
 
       // To treat commutativity of "+"
@@ -138,7 +138,7 @@ DataTree::AddUMinus(expr_t iArg1)
     {
       // Simplify -(-x) in x
       auto *uarg = dynamic_cast<UnaryOpNode *>(iArg1);
-      if (uarg != NULL && uarg->get_op_code() == oUminus)
+      if (uarg != nullptr && uarg->get_op_code() == oUminus)
         return uarg->get_arg();
 
       return AddUnaryOp(oUminus, iArg1);
diff --git a/src/DataTree.hh b/src/DataTree.hh
index b723f31c32562e567204ddaf33c50ef58311f93a..117daa0707b686fa14397b6bac8cf2c09f89b34c 100644
--- a/src/DataTree.hh
+++ b/src/DataTree.hh
@@ -337,7 +337,7 @@ DataTree::AddUnaryOp(UnaryOpcode op_code, expr_t arg, int arg_exp_info_set, int
   // Try to reduce to a constant
   // Case where arg is a constant and op_code == oUminus (i.e. we're adding a negative constant) is skipped
   auto *carg = dynamic_cast<NumConstNode *>(arg);
-  if (op_code != oUminus || carg == NULL)
+  if (op_code != oUminus || carg == nullptr)
     {
       try
         {
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index ab5ed672bccfe67f0b2c42e9827b08de65e997df..b4ea1895282bf66093c732e2d5be1deedcc7d0c7 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -215,7 +215,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
 {
   string tmp_s, sps;
   ostringstream tmp_output, tmp1_output, global_output;
-  expr_t lhs = NULL, rhs = NULL;
+  expr_t lhs = nullptr, rhs = nullptr;
   BinaryOpNode *eq_node;
   ostringstream Ufoss;
   vector<string> Uf(symbol_table.endo_nbr(), "");
@@ -457,7 +457,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &dynamic_basename) const
               output << "  " << "% //Temporary variables" << endl;
               for (auto it : v_temporary_terms[block][i])
                 {
-                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != NULL)
+                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
                     it->writeExternalFunctionOutput(output, local_output_type, tt2, temporary_terms_idxs, tef_terms);
 
                   output << "  " <<  sps;
@@ -1063,7 +1063,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
   ostringstream tmp_output;
   ofstream code_file;
   unsigned int instruction_number = 0;
-  expr_t lhs = NULL, rhs = NULL;
+  expr_t lhs = nullptr, rhs = nullptr;
   BinaryOpNode *eq_node;
   Uff Uf[symbol_table.endo_nbr()];
   map<expr_t, int> reference_count;
@@ -1198,7 +1198,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
             {
               for (auto it : v_temporary_terms[block][i])
                 {
-                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != NULL)
+                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
                     it->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, true, false, tef_terms);
 
                   FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find(it->idx)->second));
@@ -1266,7 +1266,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
               variable_ID = getBlockVariableID(block, i);
               equation_ID = getBlockEquationID(block, i);
               feedback_variables.push_back(variable_ID);
-              Uf[equation_ID].Ufl = NULL;
+              Uf[equation_ID].Ufl = nullptr;
               goto end;
             default:
             end:
@@ -1337,7 +1337,7 @@ DynamicModel::writeModelEquationsCode_Block(string &file_name, const string &bin
                           Uf[eqr].Ufl->pNext = (Uff_l *) malloc(sizeof(Uff_l));
                           Uf[eqr].Ufl = Uf[eqr].Ufl->pNext;
                         }
-                      Uf[eqr].Ufl->pNext = NULL;
+                      Uf[eqr].Ufl->pNext = nullptr;
                       Uf[eqr].Ufl->u = count_u;
                       Uf[eqr].Ufl->var = varr;
                       Uf[eqr].Ufl->lag = lag;
@@ -3742,7 +3742,7 @@ DynamicModel::getUndiffLHSForPac(vector<int> &lhs, vector<expr_t> &lhs_expr_t, v
 
       bool printerr = false;
       ExprNode::subst_table_t::const_iterator it1;
-      expr_t node = NULL;
+      expr_t node = nullptr;
       expr_t aux_var = lhs_expr_t.at(i);
       for (it1 = diff_subst_table.begin(); it1 != diff_subst_table.end(); it1++)
         if (it1->second == aux_var)
@@ -3751,7 +3751,7 @@ DynamicModel::getUndiffLHSForPac(vector<int> &lhs, vector<expr_t> &lhs_expr_t, v
             break;
           }
 
-      if (node == NULL)
+      if (node == nullptr)
         {
           cerr << "Unexpected error encountered." << endl;
           exit(EXIT_FAILURE);
@@ -3837,7 +3837,7 @@ DynamicModel::substitutePacExpectation()
   for (auto & equation : equations)
     {
       auto *substeq = dynamic_cast<BinaryOpNode *>(equation->substitutePacExpectation(subst_table));
-      assert(substeq != NULL);
+      assert(substeq != nullptr);
       equation = substeq;
     }
 
@@ -4507,7 +4507,7 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model, const boo
   for (i = 0; i < (int) equations.size(); i++)
     {
       auto *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->addMultipliersToConstraints(i));
-      assert(substeq != NULL);
+      assert(substeq != nullptr);
       equations[i] = substeq;
     }
   if (!nopreprocessoroutput)
@@ -4544,7 +4544,7 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model, const boo
   for (i = 0; i < (int) equations.size(); i++)
     for (int lag = -max_eq_lag; lag <= max_eq_lead; lag++)
       {
-        expr_t dfpower = NULL;
+        expr_t dfpower = nullptr;
         std::stringstream lagstream;
         lagstream << abs(lag);
         if (lag < 0)
@@ -5264,7 +5264,7 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
           exit(EXIT_FAILURE);
         }
       auto *substeq = dynamic_cast<BinaryOpNode *>(subst);
-      assert(substeq != NULL);
+      assert(substeq != nullptr);
       equation = substeq;
     }
 
@@ -5297,7 +5297,7 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
           exit(EXIT_FAILURE);
         }
       auto *substeq = dynamic_cast<BinaryOpNode *>(subst);
-      assert(substeq != NULL);
+      assert(substeq != nullptr);
       aux_equation = substeq;
     }
 
@@ -5330,7 +5330,7 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
           exit(EXIT_FAILURE);
         }
       auto *substeq = dynamic_cast<BinaryOpNode *>(subst);
-      assert(substeq != NULL);
+      assert(substeq != nullptr);
       diff_aux_equation = substeq;
     }
 
@@ -5404,7 +5404,7 @@ DynamicModel::substituteUnaryOps(StaticModel &static_model)
     {
       auto *substeq = dynamic_cast<BinaryOpNode *>(equation->
                                                            substituteUnaryOpNodes(static_model, nodes, subst_table, neweqs));
-      assert(substeq != NULL);
+      assert(substeq != nullptr);
       equation = substeq;
     }
 
@@ -5439,7 +5439,7 @@ DynamicModel::substituteDiff(StaticModel &static_model, ExprNode::subst_table_t
     {
       auto *substeq = dynamic_cast<BinaryOpNode *>(equation->
                                                            substituteDiff(static_model, diff_table, diff_subst_table, neweqs));
-      assert(substeq != NULL);
+      assert(substeq != nullptr);
       equation = substeq;
     }
 
@@ -5473,7 +5473,7 @@ DynamicModel::substituteExpectation(bool partial_information_model)
   for (auto & equation : equations)
     {
       auto *substeq = dynamic_cast<BinaryOpNode *>(equation->substituteExpectation(subst_table, neweqs, partial_information_model));
-      assert(substeq != NULL);
+      assert(substeq != nullptr);
       equation = substeq;
     }
 
@@ -5502,7 +5502,7 @@ DynamicModel::transformPredeterminedVariables()
   for (auto & equation : equations)
     {
       auto *substeq = dynamic_cast<BinaryOpNode *>(equation->decreaseLeadsLagsPredeterminedVariables());
-      assert(substeq != NULL);
+      assert(substeq != nullptr);
       equation = substeq;
     }
 }
@@ -5516,14 +5516,14 @@ DynamicModel::detrendEquations()
     for (auto & equation : equations)
       {
         auto *substeq = dynamic_cast<BinaryOpNode *>(equation->detrend(it->first, it->second.first, it->second.second));
-        assert(substeq != NULL);
+        assert(substeq != nullptr);
         equation = dynamic_cast<BinaryOpNode *>(substeq);
       }
 
   for (auto & equation : equations)
     {
       BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equation->removeTrendLeadLag(trend_symbols_map));
-      assert(substeq != NULL);
+      assert(substeq != nullptr);
       equation = dynamic_cast<BinaryOpNode *>(substeq);
     }
 }
@@ -5534,7 +5534,7 @@ DynamicModel::removeTrendVariableFromEquations()
   for (auto & equation : equations)
     {
       auto *substeq = dynamic_cast<BinaryOpNode *>(equation->replaceTrendVar());
-      assert(substeq != NULL);
+      assert(substeq != nullptr);
       equation = dynamic_cast<BinaryOpNode *>(substeq);
     }
 }
@@ -5553,7 +5553,7 @@ DynamicModel::fillEvalContext(eval_context_t &eval_context) const
     {
       assert(aux_equation->get_op_code() == oEqual);
       auto *auxvar = dynamic_cast<VariableNode *>(aux_equation->get_arg1());
-      assert(auxvar != NULL);
+      assert(auxvar != nullptr);
       try
         {
           double val = aux_equation->get_arg2()->eval(eval_context);
@@ -5604,7 +5604,7 @@ void
 DynamicModel::addStaticOnlyEquation(expr_t eq, int lineno, const vector<pair<string, string> > &eq_tags)
 {
   auto *beq = dynamic_cast<BinaryOpNode *>(eq);
-  assert(beq != NULL && beq->get_op_code() == oEqual);
+  assert(beq != nullptr && beq->get_op_code() == oEqual);
 
   vector<pair<string, string> > soe_eq_tags;
   for (const auto & eq_tag : eq_tags)
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 05e4dff790887426294f01b3ae1be5f0b51ef77a..e4a85a7d037891937e0fddc8e44c0d87ac185f67 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -177,7 +177,7 @@ pair<int, expr_t >
 ExprNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
 {
   /* nothing to do */
-  return (make_pair(0, (expr_t) NULL));
+  return (make_pair(0, (expr_t) nullptr));
 }
 
 void
@@ -257,7 +257,7 @@ ExprNode::createEndoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector
           int symb_id = datatree.symbol_table.addEndoLeadAuxiliaryVar(orig_expr->idx, substexpr);
           neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(datatree.AddVariable(symb_id, 0), substexpr)));
           substexpr = datatree.AddVariable(symb_id, +1);
-          assert(dynamic_cast<VariableNode *>(substexpr) != NULL);
+          assert(dynamic_cast<VariableNode *>(substexpr) != nullptr);
           subst_table[orig_expr] = dynamic_cast<VariableNode *>(substexpr);
         }
       else
@@ -293,7 +293,7 @@ ExprNode::createExoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector<
           int symb_id = datatree.symbol_table.addExoLeadAuxiliaryVar(orig_expr->idx, substexpr);
           neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(datatree.AddVariable(symb_id, 0), substexpr)));
           substexpr = datatree.AddVariable(symb_id, +1);
-          assert(dynamic_cast<VariableNode *>(substexpr) != NULL);
+          assert(dynamic_cast<VariableNode *>(substexpr) != nullptr);
           subst_table[orig_expr] = dynamic_cast<VariableNode *>(substexpr);
         }
       else
@@ -1184,7 +1184,7 @@ VariableNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr
     {
       if (datatree.symbol_table.getTypeSpecificID(symb_id) == var_endo && lag == 0)
         /* the endogenous variable */
-        return (make_pair(1, (expr_t) NULL));
+        return (make_pair(1, (expr_t) nullptr));
       else
         return (make_pair(0, datatree.AddVariableInternal(symb_id, lag)));
     }
@@ -1939,7 +1939,7 @@ UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id)
           if (datatree.getTypeByDerivID(deriv_id) == eParameter)
             {
               auto *varg = dynamic_cast<VariableNode *>(arg);
-              if (varg == NULL)
+              if (varg == nullptr)
                 {
                   cerr << "UnaryOpNode::composeDerivatives: STEADY_STATE() should only be used on "
                        << "standalone variables (like STEADY_STATE(y)) to be derivable w.r.t. parameters" << endl;
@@ -1960,7 +1960,7 @@ UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id)
       if (datatree.getTypeByDerivID(deriv_id) == eParameter)
         {
           auto *varg = dynamic_cast<VariableNode *>(arg);
-          assert(varg != NULL);
+          assert(varg != nullptr);
           assert(datatree.symbol_table.getType(varg->symb_id) == eEndogenous);
           return datatree.AddSteadyStateParam2ndDeriv(arg, param1_symb_id, datatree.getSymbIDByDerivID(deriv_id));
         }
@@ -2298,7 +2298,7 @@ UnaryOpNode::writeJsonOutput(ostream &output,
     case oSteadyStateParamDeriv:
       {
         auto *varg = dynamic_cast<VariableNode *>(arg);
-        assert(varg != NULL);
+        assert(varg != nullptr);
         assert(datatree.symbol_table.getType(varg->symb_id) == eEndogenous);
         assert(datatree.symbol_table.getType(param1_symb_id) == eParameter);
         int tsid_endo = datatree.symbol_table.getTypeSpecificID(varg->symb_id);
@@ -2309,7 +2309,7 @@ UnaryOpNode::writeJsonOutput(ostream &output,
     case oSteadyStateParam2ndDeriv:
       {
         auto *varg = dynamic_cast<VariableNode *>(arg);
-        assert(varg != NULL);
+        assert(varg != nullptr);
         assert(datatree.symbol_table.getType(varg->symb_id) == eEndogenous);
         assert(datatree.symbol_table.getType(param1_symb_id) == eParameter);
         assert(datatree.symbol_table.getType(param2_symb_id) == eParameter);
@@ -2464,7 +2464,7 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
     case oSteadyStateParamDeriv:
       {
         auto *varg = dynamic_cast<VariableNode *>(arg);
-        assert(varg != NULL);
+        assert(varg != nullptr);
         assert(datatree.symbol_table.getType(varg->symb_id) == eEndogenous);
         assert(datatree.symbol_table.getType(param1_symb_id) == eParameter);
         int tsid_endo = datatree.symbol_table.getTypeSpecificID(varg->symb_id);
@@ -2476,7 +2476,7 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
     case oSteadyStateParam2ndDeriv:
       {
         auto *varg = dynamic_cast<VariableNode *>(arg);
-        assert(varg != NULL);
+        assert(varg != nullptr);
         assert(datatree.symbol_table.getType(varg->symb_id) == eEndogenous);
         assert(datatree.symbol_table.getType(param1_symb_id) == eParameter);
         assert(datatree.symbol_table.getType(param2_symb_id) == eParameter);
@@ -2694,7 +2694,7 @@ UnaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_
   expr_t New_expr_t = res.second;
 
   if (is_endogenous_present == 2) /* The equation could not be normalized and the process is given-up*/
-    return (make_pair(2, (expr_t) NULL));
+    return (make_pair(2, (expr_t) nullptr));
   else if (is_endogenous_present) /* The argument of the function contains the current values of
                                      the endogenous variable associated to the equation.
                                      In order to normalized, we have to apply the invert function to the RHS.*/
@@ -2702,67 +2702,67 @@ UnaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_
       switch (op_code)
         {
         case oUminus:
-          List_of_Op_RHS.emplace_back(oUminus, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oUminus, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oExp:
-          List_of_Op_RHS.emplace_back(oLog, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oLog, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oLog:
-          List_of_Op_RHS.emplace_back(oExp, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oExp, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oLog10:
-          List_of_Op_RHS.emplace_back(oPower, make_pair((expr_t) NULL, datatree.AddNonNegativeConstant("10")));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oPower, make_pair((expr_t) nullptr, datatree.AddNonNegativeConstant("10")));
+          return (make_pair(1, (expr_t) nullptr));
         case oCos:
-          List_of_Op_RHS.emplace_back(oAcos, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oAcos, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oSin:
-          List_of_Op_RHS.emplace_back(oAsin, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oAsin, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oTan:
-          List_of_Op_RHS.emplace_back(oAtan, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oAtan, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oAcos:
-          List_of_Op_RHS.emplace_back(oCos, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oCos, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oAsin:
-          List_of_Op_RHS.emplace_back(oSin, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oSin, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oAtan:
-          List_of_Op_RHS.emplace_back(oTan, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oTan, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oCosh:
-          List_of_Op_RHS.emplace_back(oAcosh, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oAcosh, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oSinh:
-          List_of_Op_RHS.emplace_back(oAsinh, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oAsinh, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oTanh:
-          List_of_Op_RHS.emplace_back(oAtanh, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oAtanh, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oAcosh:
-          List_of_Op_RHS.emplace_back(oCosh, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oCosh, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oAsinh:
-          List_of_Op_RHS.emplace_back(oSinh, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oSinh, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oAtanh:
-          List_of_Op_RHS.emplace_back(oTanh, make_pair((expr_t) NULL, (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oTanh, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         case oSqrt:
-          List_of_Op_RHS.emplace_back(oPower, make_pair((expr_t) NULL, datatree.Two));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oPower, make_pair((expr_t) nullptr, datatree.Two));
+          return (make_pair(1, (expr_t) nullptr));
         case oAbs:
-          return (make_pair(2, (expr_t) NULL));
+          return (make_pair(2, (expr_t) nullptr));
         case oSign:
-          return (make_pair(2, (expr_t) NULL));
+          return (make_pair(2, (expr_t) nullptr));
         case oSteadyState:
-          return (make_pair(2, (expr_t) NULL));
+          return (make_pair(2, (expr_t) nullptr));
         case oErf:
-          return (make_pair(2, (expr_t) NULL));
+          return (make_pair(2, (expr_t) nullptr));
         default:
           cerr << "Unary operator not handled during the normalization process" << endl;
-          return (make_pair(2, (expr_t) NULL)); // Could not be normalized
+          return (make_pair(2, (expr_t) nullptr)); // Could not be normalized
         }
     }
   else
@@ -2814,7 +2814,7 @@ UnaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_
           return (make_pair(0, datatree.AddErf(New_expr_t)));
         default:
           cerr << "Unary operator not handled during the normalization process" << endl;
-          return (make_pair(2, (expr_t) NULL)); // Could not be normalized
+          return (make_pair(2, (expr_t) nullptr)); // Could not be normalized
         }
     }
   cerr << "UnaryOpNode::normalizeEquation: impossible case" << endl;
@@ -3005,7 +3005,7 @@ UnaryOpNode::substituteAdl() const
     }
 
   expr_t arg1subst = arg->substituteAdl();
-  expr_t retval = NULL;
+  expr_t retval = nullptr;
   ostringstream inttostr;
 
   for (auto it = adl_lags.begin(); it != adl_lags.end(); it++)
@@ -3136,7 +3136,7 @@ UnaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table,
     }
 
   int last_arg_max_lag = 0;
-  VariableNode *last_aux_var = NULL;
+  VariableNode *last_aux_var = nullptr;
   for (auto rit = it->second.rbegin();
        rit != it->second.rend(); rit++)
     {
@@ -3146,7 +3146,7 @@ UnaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table,
       auto *vn = dynamic_cast<VariableNode *>(argsubst);
       if (rit == it->second.rbegin())
         {
-          if (vn != NULL)
+          if (vn != nullptr)
             symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst, vn->get_symb_id(), vn->get_lag());
           else
             {
@@ -3168,7 +3168,7 @@ UnaryOpNode::substituteDiff(DataTree &static_datatree, diff_table_t &diff_table,
       else
         {
           // just add equation of form: AUX_DIFF = LAST_AUX_VAR(-1)
-          VariableNode *new_aux_var = NULL;
+          VariableNode *new_aux_var = nullptr;
           for (int i = last_arg_max_lag; i > rit->first; i--)
             {
               if (i == last_arg_max_lag)
@@ -3205,7 +3205,7 @@ UnaryOpNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nod
       return buildSimilarUnaryOpNode(argsubst, datatree);
     }
 
-  VariableNode *aux_var = NULL;
+  VariableNode *aux_var = nullptr;
   for (auto rit = it->second.rbegin();
        rit != it->second.rend(); rit++)
     {
@@ -3314,7 +3314,7 @@ UnaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNo
       expr_t newAuxE = datatree.AddVariable(symb_id, 0);
 
       if (partial_information_model && expectation_information_set == 0)
-        if (dynamic_cast<VariableNode *>(arg) == NULL)
+        if (dynamic_cast<VariableNode *>(arg) == nullptr)
           {
             cerr << "ERROR: In Partial Information models, EXPECTATION(0)(X) "
                  << "can only be used when X is a single variable." << endl;
@@ -3324,11 +3324,11 @@ UnaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNo
       //take care of any nested expectation operators by calling arg->substituteExpectation(.), then decreaseLeadsLags for this oExpectation operator
       //arg(lag-period) (holds entire subtree of arg(lag-period)
       expr_t substexpr = (arg->substituteExpectation(subst_table, neweqs, partial_information_model))->decreaseLeadsLags(expectation_information_set);
-      assert(substexpr != NULL);
+      assert(substexpr != nullptr);
       neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(newAuxE, substexpr))); //AUXE_period_arg.idx = arg(lag-period)
       newAuxE = datatree.AddVariable(symb_id, expectation_information_set);
 
-      assert(dynamic_cast<VariableNode *>(newAuxE) != NULL);
+      assert(dynamic_cast<VariableNode *>(newAuxE) != nullptr);
       subst_table[this] = dynamic_cast<VariableNode *>(newAuxE);
       return newAuxE;
     }
@@ -3538,7 +3538,7 @@ BinaryOpNode::composeDerivatives(expr_t darg1, expr_t darg2)
         if (darg1 == datatree.Zero)
           return datatree.Zero;
         else
-          if (dynamic_cast<NumConstNode *>(arg2) != NULL)
+          if (dynamic_cast<NumConstNode *>(arg2) != nullptr)
             {
               t11 = datatree.AddMinus(arg2, datatree.One);
               t12 = datatree.AddPower(arg1, t11);
@@ -4023,7 +4023,7 @@ BinaryOpNode::writeJsonOutput(ostream &output,
   // add parenthesis around left argument
   auto *barg1 = dynamic_cast<BinaryOpNode *>(arg1);
   if (arg1->precedenceJson(temporary_terms) < prec
-      || (op_code == oPower && barg1 != NULL && barg1->op_code == oPower))
+      || (op_code == oPower && barg1 != nullptr && barg1->op_code == oPower))
     {
       output << "(";
       close_parenthesis = true;
@@ -4088,7 +4088,7 @@ BinaryOpNode::writeJsonOutput(ostream &output,
   auto *barg2 = dynamic_cast<BinaryOpNode *>(arg2);
   int arg2_prec = arg2->precedenceJson(temporary_terms);
   if (arg2_prec < prec
-      || (op_code == oPower && barg2 != NULL && barg2->op_code == oPower)
+      || (op_code == oPower && barg2 != nullptr && barg2->op_code == oPower)
       || (op_code == oMinus && arg2_prec == prec)
       || (op_code == oDivide && arg2_prec == prec))
     {
@@ -4166,7 +4166,7 @@ BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       // If left argument has a lower precedence, or if current and left argument are both power operators, add parenthesis around left argument
       auto *barg1 = dynamic_cast<BinaryOpNode *>(arg1);
       if (arg1->precedence(output_type, temporary_terms) < prec
-          || (op_code == oPower && barg1 != NULL && barg1->op_code == oPower))
+          || (op_code == oPower && barg1 != nullptr && barg1->op_code == oPower))
         {
           output << LEFT_PAR(output_type);
           close_parenthesis = true;
@@ -4257,7 +4257,7 @@ BinaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
       auto *barg2 = dynamic_cast<BinaryOpNode *>(arg2);
       int arg2_prec = arg2->precedence(output_type, temporary_terms);
       if (arg2_prec < prec
-          || (op_code == oPower && barg2 != NULL && barg2->op_code == oPower && !IS_LATEX(output_type))
+          || (op_code == oPower && barg2 != nullptr && barg2->op_code == oPower && !IS_LATEX(output_type))
           || (op_code == oMinus && arg2_prec == prec)
           || (op_code == oDivide && arg2_prec == prec && !IS_LATEX(output_type)))
         {
@@ -4385,7 +4385,7 @@ BinaryOpNode::Compute_RHS(expr_t arg1, expr_t arg2, int op, int op_type) const
         }
       break;
     }
-  return ((expr_t) NULL);
+  return ((expr_t) nullptr);
 }
 
 pair<int, expr_t>
@@ -4408,9 +4408,9 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr
   /* If the two expressions contains the current value of the endogenous variable associated to the equation
      the equation could not be normalized and the process is given-up.*/
   if (is_endogenous_present_1 == 2 || is_endogenous_present_2 == 2)
-    return (make_pair(2, (expr_t) NULL));
+    return (make_pair(2, (expr_t) nullptr));
   else if (is_endogenous_present_1 && is_endogenous_present_2)
-    return (make_pair(2, (expr_t) NULL));
+    return (make_pair(2, (expr_t) nullptr));
   else if (is_endogenous_present_1) /*If the current values of the endogenous variable associated to the equation
                                       is present only in the first operand of the expression, we try to normalize the equation*/
     {
@@ -4464,39 +4464,39 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr
     case oPlus:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         {
-          List_of_Op_RHS.emplace_back(oMinus, make_pair(datatree.AddPlus(expr_t_1, expr_t_2), (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oMinus, make_pair(datatree.AddPlus(expr_t_1, expr_t_2), (expr_t) nullptr));
           return (make_pair(0, datatree.AddPlus(expr_t_1, expr_t_2)));
         }
       else if (is_endogenous_present_1 && is_endogenous_present_2)
-        return (make_pair(1, (expr_t) NULL));
+        return (make_pair(1, (expr_t) nullptr));
       else if (!is_endogenous_present_1 && is_endogenous_present_2)
         {
-          List_of_Op_RHS.emplace_back(oMinus, make_pair(expr_t_1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oMinus, make_pair(expr_t_1, (expr_t) nullptr));
           return (make_pair(1, expr_t_1));
         }
       else if (is_endogenous_present_1 && !is_endogenous_present_2)
         {
-          List_of_Op_RHS.emplace_back(oMinus, make_pair(expr_t_2, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oMinus, make_pair(expr_t_2, (expr_t) nullptr));
           return (make_pair(1, expr_t_2));
         }
       break;
     case oMinus:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         {
-          List_of_Op_RHS.emplace_back(oMinus, make_pair(datatree.AddMinus(expr_t_1, expr_t_2), (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oMinus, make_pair(datatree.AddMinus(expr_t_1, expr_t_2), (expr_t) nullptr));
           return (make_pair(0, datatree.AddMinus(expr_t_1, expr_t_2)));
         }
       else if (is_endogenous_present_1 && is_endogenous_present_2)
-        return (make_pair(1, (expr_t) NULL));
+        return (make_pair(1, (expr_t) nullptr));
       else if (!is_endogenous_present_1 && is_endogenous_present_2)
         {
-          List_of_Op_RHS.emplace_back(oUminus, make_pair((expr_t) NULL, (expr_t) NULL));
-          List_of_Op_RHS.emplace_back(oMinus, make_pair(expr_t_1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oUminus, make_pair((expr_t) nullptr, (expr_t) nullptr));
+          List_of_Op_RHS.emplace_back(oMinus, make_pair(expr_t_1, (expr_t) nullptr));
           return (make_pair(1, expr_t_1));
         }
       else if (is_endogenous_present_1 && !is_endogenous_present_2)
         {
-          List_of_Op_RHS.emplace_back(oPlus, make_pair(expr_t_2, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oPlus, make_pair(expr_t_2, (expr_t) nullptr));
           return (make_pair(1, datatree.AddUMinus(expr_t_2)));
         }
       break;
@@ -4505,49 +4505,49 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr
         return (make_pair(0, datatree.AddTimes(expr_t_1, expr_t_2)));
       else if (!is_endogenous_present_1 && is_endogenous_present_2)
         {
-          List_of_Op_RHS.emplace_back(oDivide, make_pair(expr_t_1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oDivide, make_pair(expr_t_1, (expr_t) nullptr));
           return (make_pair(1, expr_t_1));
         }
       else if (is_endogenous_present_1 && !is_endogenous_present_2)
         {
-          List_of_Op_RHS.emplace_back(oDivide, make_pair(expr_t_2, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oDivide, make_pair(expr_t_2, (expr_t) nullptr));
           return (make_pair(1, expr_t_2));
         }
       else
-        return (make_pair(1, (expr_t) NULL));
+        return (make_pair(1, (expr_t) nullptr));
       break;
     case oDivide:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         return (make_pair(0, datatree.AddDivide(expr_t_1, expr_t_2)));
       else if (!is_endogenous_present_1 && is_endogenous_present_2)
         {
-          List_of_Op_RHS.emplace_back(oDivide, make_pair((expr_t) NULL, expr_t_1));
+          List_of_Op_RHS.emplace_back(oDivide, make_pair((expr_t) nullptr, expr_t_1));
           return (make_pair(1, expr_t_1));
         }
       else if (is_endogenous_present_1 && !is_endogenous_present_2)
         {
-          List_of_Op_RHS.emplace_back(oTimes, make_pair(expr_t_2, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oTimes, make_pair(expr_t_2, (expr_t) nullptr));
           return (make_pair(1, expr_t_2));
         }
       else
-        return (make_pair(1, (expr_t) NULL));
+        return (make_pair(1, (expr_t) nullptr));
       break;
     case oPower:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         return (make_pair(0, datatree.AddPower(expr_t_1, expr_t_2)));
       else if (is_endogenous_present_1 && !is_endogenous_present_2)
         {
-          List_of_Op_RHS.emplace_back(oPower, make_pair(datatree.AddDivide(datatree.One, expr_t_2), (expr_t) NULL));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oPower, make_pair(datatree.AddDivide(datatree.One, expr_t_2), (expr_t) nullptr));
+          return (make_pair(1, (expr_t) nullptr));
         }
       else if (!is_endogenous_present_1 && is_endogenous_present_2)
         {
           /* we have to nomalize a^f(X) = RHS */
           /* First computes the ln(RHS)*/
-          List_of_Op_RHS.emplace_back(oLog, make_pair((expr_t) NULL, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oLog, make_pair((expr_t) nullptr, (expr_t) nullptr));
           /* Second  computes f(X) = ln(RHS) / ln(a)*/
-          List_of_Op_RHS.emplace_back(oDivide, make_pair((expr_t) NULL, datatree.AddLog(expr_t_1)));
-          return (make_pair(1, (expr_t) NULL));
+          List_of_Op_RHS.emplace_back(oDivide, make_pair((expr_t) nullptr, datatree.AddLog(expr_t_1)));
+          return (make_pair(1, (expr_t) nullptr));
         }
       break;
     case oEqual:
@@ -4580,53 +4580,53 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         return (make_pair(0, datatree.AddMax(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (expr_t) NULL));
+        return (make_pair(1, (expr_t) nullptr));
       break;
     case oMin:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         return (make_pair(0, datatree.AddMin(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (expr_t) NULL));
+        return (make_pair(1, (expr_t) nullptr));
       break;
     case oLess:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         return (make_pair(0, datatree.AddLess(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (expr_t) NULL));
+        return (make_pair(1, (expr_t) nullptr));
       break;
     case oGreater:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         return (make_pair(0, datatree.AddGreater(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (expr_t) NULL));
+        return (make_pair(1, (expr_t) nullptr));
       break;
     case oLessEqual:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         return (make_pair(0, datatree.AddLessEqual(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (expr_t) NULL));
+        return (make_pair(1, (expr_t) nullptr));
       break;
     case oGreaterEqual:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         return (make_pair(0, datatree.AddGreaterEqual(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (expr_t) NULL));
+        return (make_pair(1, (expr_t) nullptr));
       break;
     case oEqualEqual:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         return (make_pair(0, datatree.AddEqualEqual(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (expr_t) NULL));
+        return (make_pair(1, (expr_t) nullptr));
       break;
     case oDifferent:
       if (!is_endogenous_present_1 && !is_endogenous_present_2)
         return (make_pair(0, datatree.AddDifferent(expr_t_1, expr_t_2)));
       else
-        return (make_pair(1, (expr_t) NULL));
+        return (make_pair(1, (expr_t) nullptr));
       break;
     default:
       cerr << "Binary operator not handled during the normalization process" << endl;
-      return (make_pair(2, (expr_t) NULL)); // Could not be normalized
+      return (make_pair(2, (expr_t) nullptr)); // Could not be normalized
     }
   // Suppress GCC warning
   cerr << "BinaryOpNode::normalizeEquation: impossible case" << endl;
@@ -5030,11 +5030,11 @@ BinaryOpNode::walkPacParametersHelper(const expr_t arg1, const expr_t arg2,
   else if (endogs.size() >= 2)
     {
       auto *testarg2 = dynamic_cast<BinaryOpNode *>(arg2);
-      if (testarg2 != NULL && testarg2->get_op_code() == oMinus)
+      if (testarg2 != nullptr && testarg2->get_op_code() == oMinus)
         {
           auto *test_arg1 = dynamic_cast<VariableNode *>(testarg2->get_arg1());
           auto *test_arg2 = dynamic_cast<VariableNode *>(testarg2->get_arg2());
-          if (test_arg1 != NULL && test_arg2 != NULL && lhs.first != -1)
+          if (test_arg1 != nullptr && test_arg2 != nullptr && lhs.first != -1)
             {
               test_arg1->collectDynamicVariables(eEndogenous, endogs);
               ec_params_and_vars.insert(make_pair(*(params.begin()), *(endogs.begin())));
@@ -5590,7 +5590,7 @@ TrinaryOpNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, exp
   if (!is_endogenous_present_1 && !is_endogenous_present_2 && !is_endogenous_present_3)
     return (make_pair(0, datatree.AddNormcdf(expr_t_1, expr_t_2, expr_t_3)));
   else
-    return (make_pair(1, (expr_t) NULL));
+    return (make_pair(1, (expr_t) nullptr));
 }
 
 expr_t
@@ -6453,7 +6453,7 @@ AbstractExternalFunctionNode::normalizeEquation(int var_endo, vector<pair<int, p
   if (!present)
     return (make_pair(0, datatree.AddExternalFunction(symb_id, V_expr_t)));
   else
-    return (make_pair(1, (expr_t) NULL));
+    return (make_pair(1, (expr_t) nullptr));
 }
 
 void
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 31ed14e43dc9c7ea5297d38fe3c61fc17ee0f3cb..b4ab5d00d1f27270edfd454176ea452c0c4d3bfd 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -382,7 +382,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
        it != statements.end(); it++)
     {
       auto *vms = dynamic_cast<VarModelStatement *>(*it);
-      if (vms != NULL)
+      if (vms != nullptr)
         {
           string var_model_name;
           vms->getVarModelInfo(var_model_name, var_model_info_var_expectation, var_model_eq_tags);
@@ -404,7 +404,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
                       make_pair(make_pair(max_lag, nonstationary), eqnumber));
         }
       auto *pms = dynamic_cast<PacModelStatement *>(*it);
-      if (pms != NULL)
+      if (pms != nullptr)
          {
            pair<string, pair<string, pair<string, pair<int, map<string, int> > > > >
              pac_model_info_pac_expectation;
@@ -459,14 +459,14 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
   mod_file_struct.orig_eq_nbr = dynamic_model.equation_number();
   if (mod_file_struct.ramsey_model_present)
     {
-      StaticModel *planner_objective = NULL;
+      StaticModel *planner_objective = nullptr;
       for (auto & statement : statements)
         {
           auto *pos = dynamic_cast<PlannerObjectiveStatement *>(statement);
-          if (pos != NULL)
+          if (pos != nullptr)
             planner_objective = pos->getPlannerObjective();
         }
-      assert(planner_objective != NULL);
+      assert(planner_objective != nullptr);
 
       /*
         clone the model then clone the new equations back to the original because
@@ -524,7 +524,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
        it != statements.end(); it++)
     {
       auto *vms = dynamic_cast<VarModelStatement *>(*it);
-      if (vms != NULL)
+      if (vms != nullptr)
         {
           string var_model_name;
           vms->getVarModelInfo(var_model_name, var_model_info_var_expectation, var_model_eq_tags);
@@ -598,7 +598,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
     for (auto & statement : statements)
       {
         auto *rps = dynamic_cast<RamseyPolicyStatement *>(statement);
-        if (rps != NULL)
+        if (rps != nullptr)
           rps->checkRamseyPolicyList();
       }
 
@@ -989,7 +989,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
       /* Special treatment for initval block: insert initial values for the
          auxiliary variables and initialize exo det */
       auto *ivs = dynamic_cast<InitValStatement *>(statement);
-      if (ivs != NULL)
+      if (ivs != nullptr)
         {
           static_model.writeAuxVarInitval(mOutputFile, oMatlabOutsideModel);
           ivs->writeOutputPostInit(mOutputFile);
@@ -997,7 +997,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
 
       // Special treatment for endval block: insert initial values for the auxiliary variables
       auto *evs = dynamic_cast<EndValStatement *>(statement);
-      if (evs != NULL)
+      if (evs != nullptr)
         static_model.writeAuxVarInitval(mOutputFile, oMatlabOutsideModel);
 
       // Special treatment for load params and steady state statement: insert initial values for the auxiliary variables
@@ -1007,7 +1007,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
 
       // Special treatement for Var Models
       auto *vms = dynamic_cast<VarModelStatement *>(statement);
-      if (vms != NULL)
+      if (vms != nullptr)
         vms->createVarModelMFunction(mOutputFile, dynamic_model.getVarExpectationFunctionsToWrite());
     }
 
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index a00e611b925565611137923027451477db1d93ef..ab8b8016454fa40274251a6a6d432a00dc5faca1 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -204,7 +204,7 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
               if (static_jacobian.find(make_pair(it->first.first, it->first.second)) == static_jacobian.end())
                 static_jacobian[make_pair(it->first.first, it->first.second)] = 0;
               if (dynamic_jacobian.find(make_pair(0, make_pair(it->first.first, it->first.second))) == dynamic_jacobian.end())
-                dynamic_jacobian[make_pair(0, make_pair(it->first.first, it->first.second))] = 0;
+                dynamic_jacobian[make_pair(0, make_pair(it->first.first, it->first.second))] = nullptr;
               if (contemporaneous_jacobian.find(make_pair(it->first.first, it->first.second)) == contemporaneous_jacobian.end())
                 contemporaneous_jacobian[make_pair(it->first.first, it->first.second)] = 0;
               try
@@ -235,7 +235,7 @@ ModelTree::computeNormalizedEquations(multimap<int, int> &endo2eqs) const
   for (size_t i = 0; i < equations.size(); i++)
     {
       auto *lhs = dynamic_cast<VariableNode *>(equations[i]->get_arg1());
-      if (lhs == NULL)
+      if (lhs == nullptr)
         continue;
 
       int symb_id = lhs->get_symb_id();
@@ -1224,7 +1224,7 @@ ModelTree::writeTemporaryTerms(const temporary_terms_t &tt,
   for (auto it = tt.begin();
        it != tt.end(); it++)
     {
-      if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
+      if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != nullptr)
         (*it)->writeExternalFunctionOutput(output, output_type, tt2, tt_idxs, tef_terms);
 
       if (IS_C(output_type))
@@ -1257,7 +1257,7 @@ ModelTree::writeJsonTemporaryTerms(const temporary_terms_t &tt, const temporary_
   for (auto it : tt)
     if (ttm1.find(it) == ttm1.end())
       {
-        if (dynamic_cast<AbstractExternalFunctionNode *>(it) != NULL)
+        if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
           {
             if (wrote_term)
               output << ", ";
@@ -1429,7 +1429,7 @@ ModelTree::compileTemporaryTerms(ostream &code_file, unsigned int &instruction_n
   deriv_node_temp_terms_t tef_terms;
   for (auto it : tt)
     {
-      if (dynamic_cast<AbstractExternalFunctionNode *>(it) != NULL)
+      if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
         {
           it->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, dynamic, steady_dynamic, tef_terms);
         }
@@ -1734,7 +1734,7 @@ void
 ModelTree::addEquation(expr_t eq, int lineno)
 {
   auto *beq = dynamic_cast<BinaryOpNode *>(eq);
-  assert(beq != NULL && beq->get_op_code() == oEqual);
+  assert(beq != nullptr && beq->get_op_code() == oEqual);
 
   equations.push_back(beq);
   equations_lineno.push_back(lineno);
@@ -1753,7 +1753,7 @@ void
 ModelTree::addAuxEquation(expr_t eq)
 {
   auto *beq = dynamic_cast<BinaryOpNode *>(eq);
-  assert(beq != NULL && beq->get_op_code() == oEqual);
+  assert(beq != nullptr && beq->get_op_code() == oEqual);
 
   aux_equations.push_back(beq);
 }
diff --git a/src/NumericalConstants.cc b/src/NumericalConstants.cc
index 5ec3aa4214781b1a0d6320ce318c095d237a19a4..bd69d3703c76dc02abcac13ebfd4f9beb28018a5 100644
--- a/src/NumericalConstants.cc
+++ b/src/NumericalConstants.cc
@@ -37,7 +37,7 @@ NumericalConstants::AddNonNegativeConstant(const string &iConst)
   mNumericalConstants.push_back(iConst);
   numConstantsIndex[iConst] = id;
 
-  double val = strtod(iConst.c_str(), NULL);
+  double val = strtod(iConst.c_str(), nullptr);
 
   /* Note that we allow underflows (will be converted to 0) and overflows (will
      be converted to Inf), as MATLAB and Octave do. */
diff --git a/src/NumericalInitialization.cc b/src/NumericalInitialization.cc
index 9c293a7cd6a00089f56f2d8559b5d32b659818fe..2f93f9d0a9badb9a15778357942cf15044139ea0 100644
--- a/src/NumericalInitialization.cc
+++ b/src/NumericalInitialization.cc
@@ -495,7 +495,7 @@ HomotopyStatement::writeOutput(ostream &output, const string &basename, bool min
       const int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
 
       output << "options_.homotopy_values = vertcat(options_.homotopy_values, [ " << type << ", " << tsid << ", ";
-      if (expression1 != NULL)
+      if (expression1 != nullptr)
         expression1->writeOutput(output);
       else
         output << "NaN";
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index fd13292432034e456d2b9d8d23ffa3f254889b48..1662257f7187102c3a7b9a590544ab085461b76c 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -187,13 +187,13 @@ ParsingDriver::declare_symbol(const string *name, SymbolType type, const string
 {
   try
     {
-      if (tex_name == NULL && partition_value == NULL)
+      if (tex_name == nullptr && partition_value == nullptr)
         mod_file->symbol_table.addSymbol(*name, type);
       else
-        if (tex_name == NULL)
+        if (tex_name == nullptr)
           mod_file->symbol_table.addSymbol(*name, type, "", partition_value);
-        else if (partition_value == NULL)
-          mod_file->symbol_table.addSymbol(*name, type, *tex_name, NULL);
+        else if (partition_value == nullptr)
+          mod_file->symbol_table.addSymbol(*name, type, *tex_name, nullptr);
         else
           mod_file->symbol_table.addSymbol(*name, type, *tex_name, partition_value);
     }
@@ -211,9 +211,9 @@ ParsingDriver::declare_endogenous(string *name, string *tex_name, vector<pair<st
 {
   declare_symbol(name, eEndogenous, tex_name, partition_value);
   delete name;
-  if (tex_name != NULL)
+  if (tex_name != nullptr)
     delete tex_name;
-  if (partition_value != NULL)
+  if (partition_value != nullptr)
     {
       for (auto & it : *partition_value)
         {
@@ -238,7 +238,7 @@ ParsingDriver::declare_var_endogenous(string *name)
       return;
     }
 
-  declare_symbol(name, eEndogenousVAR, NULL, NULL);
+  declare_symbol(name, eEndogenousVAR, nullptr, nullptr);
   add_in_symbol_list(name);
 }
 
@@ -247,9 +247,9 @@ ParsingDriver::declare_exogenous(string *name, string *tex_name, vector<pair<str
 {
   declare_symbol(name, eExogenous, tex_name, partition_value);
   delete name;
-  if (tex_name != NULL)
+  if (tex_name != nullptr)
     delete tex_name;
-  if (partition_value != NULL)
+  if (partition_value != nullptr)
     {
       for (auto & it : *partition_value)
         {
@@ -266,9 +266,9 @@ ParsingDriver::declare_exogenous_det(string *name, string *tex_name, vector<pair
 {
   declare_symbol(name, eExogenousDet, tex_name, partition_value);
   delete name;
-  if (tex_name != NULL)
+  if (tex_name != nullptr)
     delete tex_name;
-  if (partition_value != NULL)
+  if (partition_value != nullptr)
     {
       for (auto & it : *partition_value)
         {
@@ -285,9 +285,9 @@ ParsingDriver::declare_parameter(string *name, string *tex_name, vector<pair<str
 {
   declare_symbol(name, eParameter, tex_name, partition_value);
   delete name;
-  if (tex_name != NULL)
+  if (tex_name != nullptr)
     delete tex_name;
-  if (partition_value != NULL)
+  if (partition_value != nullptr)
     {
       for (auto & it : *partition_value)
         {
@@ -305,7 +305,7 @@ ParsingDriver::declare_statement_local_variable(string *name)
   if (mod_file->symbol_table.exists(*name))
     error("Symbol " + *name + " cannot be assigned within a statement "
           +"while being assigned elsewhere in the modfile");
-  declare_symbol(name, eStatementDeclaredVariable, NULL, NULL);
+  declare_symbol(name, eStatementDeclaredVariable, nullptr, nullptr);
   delete name;
 }
 
@@ -314,7 +314,7 @@ ParsingDriver::declare_optimal_policy_discount_factor_parameter(expr_t exprnode)
 {
   string *optimalParName_declare = new string("optimal_policy_discount_factor");
   string *optimalParName_init = new string("optimal_policy_discount_factor");
-  declare_parameter(optimalParName_declare, NULL);
+  declare_parameter(optimalParName_declare, nullptr);
   init_param(optimalParName_init, exprnode);
 }
 
@@ -327,10 +327,10 @@ ParsingDriver::begin_trend()
 void
 ParsingDriver::declare_trend_var(bool log_trend, string *name, string *tex_name)
 {
-  declare_symbol(name, log_trend ? eLogTrend : eTrend, tex_name, NULL);
+  declare_symbol(name, log_trend ? eLogTrend : eTrend, tex_name, nullptr);
   declared_trend_vars.push_back(mod_file->symbol_table.getID(*name));
   delete name;
-  if (tex_name != NULL)
+  if (tex_name != nullptr)
     delete tex_name;
 }
 
@@ -481,10 +481,10 @@ ParsingDriver::add_model_variable(int symb_id, int lag)
   if (type == eModelLocalVariable && lag != 0)
     error("Model local variable " + mod_file->symbol_table.getName(symb_id) + " cannot be given a lead or a lag.");
 
-  if (dynamic_cast<StaticModel *>(model_tree) != NULL && lag != 0)
+  if (dynamic_cast<StaticModel *>(model_tree) != nullptr && lag != 0)
     error("Leads and lags on variables are forbidden in 'planner_objective'.");
 
-  if (dynamic_cast<StaticModel *>(model_tree) != NULL && type == eModelLocalVariable)
+  if (dynamic_cast<StaticModel *>(model_tree) != nullptr && type == eModelLocalVariable)
     error("Model local variable " + mod_file->symbol_table.getName(symb_id) + " cannot be used in 'planner_objective'.");
 
   // It makes sense to allow a lead/lag on parameters: during steady state calibration, endogenous and parameters can be swapped
@@ -519,12 +519,12 @@ ParsingDriver::add_expression_variable(string *name)
 void
 ParsingDriver::declare_nonstationary_var(string *name, string *tex_name, vector<pair<string *, string *> *> *partition_value)
 {
-  if (tex_name == NULL && partition_value == NULL)
+  if (tex_name == nullptr && partition_value == nullptr)
     declare_endogenous(new string(*name));
   else
-    if (tex_name == NULL)
-      declare_endogenous(new string(*name), NULL, partition_value);
-    else if (partition_value == NULL)
+    if (tex_name == nullptr)
+      declare_endogenous(new string(*name), nullptr, partition_value);
+    else if (partition_value == nullptr)
       declare_endogenous(new string(*name), tex_name);
     else
       declare_endogenous(new string(*name), tex_name, partition_value);
@@ -609,13 +609,13 @@ void
 ParsingDriver::add_VAR_restriction_coeff(string *name1, string *name2, string *lagstr)
 {
   int symb_id1 = mod_file->symbol_table.getID(*name1);
-  int symb_id2 = name2 == NULL ? -1 : mod_file->symbol_table.getID(*name2);
+  int symb_id2 = name2 == nullptr ? -1 : mod_file->symbol_table.getID(*name2);
   int lag = atoi(lagstr->c_str());
 
   var_restriction_coeff = make_pair(symb_id1, make_pair(symb_id2, lag));
 
   delete name1;
-  if (name2 != NULL)
+  if (name2 != nullptr)
     delete name2;
   delete lagstr;
 }
@@ -633,7 +633,7 @@ ParsingDriver::add_VAR_restriction_equation_or_crossequation(string *numberstr)
   double number = atof(numberstr->c_str());
   if (var_restriction_eq_or_crosseq.size() == 1)
     var_restriction_equation_or_crossequation = make_pair(make_pair(var_restriction_eq_or_crosseq[0],
-                                                                    make_pair(make_pair(-1, make_pair(-1, -1)), (expr_t) NULL)),
+                                                                    make_pair(make_pair(-1, make_pair(-1, -1)), (expr_t) nullptr)),
                                                           number);
   else
     var_restriction_equation_or_crossequation = make_pair(make_pair(var_restriction_eq_or_crosseq[0],
@@ -654,7 +654,7 @@ ParsingDriver::multiply_arg2_by_neg_one()
 void
 ParsingDriver::add_VAR_restriction_equation_or_crossequation_final(string *name)
 {
-  if (name != NULL)
+  if (name != nullptr)
     {
       int symb_id = mod_file->symbol_table.getID(*name);
       equation_restrictions[symb_id] = var_restriction_equation_or_crossequation;
@@ -2576,7 +2576,7 @@ void
 ParsingDriver::plot_conditional_forecast(string *periods)
 {
   int nperiods;
-  if (periods == NULL)
+  if (periods == nullptr)
     nperiods = -1;
   else
     {
@@ -2647,9 +2647,9 @@ ParsingDriver::add_model_equal_with_zero_rhs(expr_t arg)
 void
 ParsingDriver::declare_model_local_variable(string *name, string *tex_name)
 {
-  declare_symbol(name, eModelLocalVariable, tex_name, NULL);
+  declare_symbol(name, eModelLocalVariable, tex_name, nullptr);
   delete name;
-  if (tex_name != NULL)
+  if (tex_name != nullptr)
     delete tex_name;
 }
 
@@ -3060,7 +3060,7 @@ ParsingDriver::external_function_option(const string &name_option, const string
     {
       if (opt.empty())
         error("An argument must be passed to the 'name' option of the external_function() statement.");
-      declare_symbol(&opt, eExternalFunction, NULL, NULL);
+      declare_symbol(&opt, eExternalFunction, nullptr, nullptr);
       current_external_function_id = mod_file->symbol_table.getID(opt);
     }
   else if (name_option == "first_deriv_provided")
@@ -3069,7 +3069,7 @@ ParsingDriver::external_function_option(const string &name_option, const string
         current_external_function_options.firstDerivSymbID = eExtFunSetButNoNameProvided;
       else
         {
-          declare_symbol(&opt, eExternalFunction, NULL, NULL);
+          declare_symbol(&opt, eExternalFunction, nullptr, nullptr);
           current_external_function_options.firstDerivSymbID = mod_file->symbol_table.getID(opt);
         }
     }
@@ -3079,7 +3079,7 @@ ParsingDriver::external_function_option(const string &name_option, const string
         current_external_function_options.secondDerivSymbID = eExtFunSetButNoNameProvided;
       else
         {
-          declare_symbol(&opt, eExternalFunction, NULL, NULL);
+          declare_symbol(&opt, eExternalFunction, nullptr, nullptr);
           current_external_function_options.secondDerivSymbID = mod_file->symbol_table.getID(opt);
         }
     }
@@ -3129,12 +3129,12 @@ ParsingDriver::is_there_one_integer_argument() const
   auto *numNode = dynamic_cast<NumConstNode *>(stack_external_function_args.top().front());
   auto *unaryNode = dynamic_cast<UnaryOpNode *>(stack_external_function_args.top().front());
 
-  if (numNode == NULL && unaryNode == NULL)
+  if (numNode == nullptr && unaryNode == nullptr)
     return make_pair(false, 0);
 
   eval_context_t ectmp;
   double model_var_arg;
-  if (unaryNode == NULL)
+  if (unaryNode == nullptr)
     {
       try
         {
@@ -3228,7 +3228,7 @@ ParsingDriver::add_model_var_or_external_function(string *function_name, bool in
             error("To use an external function (" + *function_name +
                   ") within the model block, you must first declare it via the external_function() statement.");
         }
-      declare_symbol(function_name, eExternalFunction, NULL, NULL);
+      declare_symbol(function_name, eExternalFunction, nullptr, nullptr);
       current_external_function_options.nargs = stack_external_function_args.top().size();
       mod_file->external_functions_table.addExternalFunction(mod_file->symbol_table.getID(*function_name),
                                                              current_external_function_options, in_model_block);
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index ce80b945626dc0571b34b44eaff89f9357b801f2..2fe1f05f45f4def21917e83747d60979a4ca8f85 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -58,7 +58,7 @@ using namespace std;
 class DynareFlex : public DynareFlexLexer
 {
 public:
-  DynareFlex(istream *in = 0, ostream *out = 0);
+  DynareFlex(istream *in = nullptr, ostream *out = nullptr);
 
   //! The main lexing function
   Dynare::parser::token_type lex(Dynare::parser::semantic_type *yylval,
@@ -326,17 +326,17 @@ public:
   //! Sets the FILENAME for the initial value in initval
   void initval_file(string *filename);
   //! Declares an endogenous variable
-  void declare_endogenous(string *name, string *tex_name = NULL, vector<pair<string *, string *> *> *partition_value = NULL);
+  void declare_endogenous(string *name, string *tex_name = nullptr, vector<pair<string *, string *> *> *partition_value = nullptr);
   //! Declares an exogenous variable
-  void declare_exogenous(string *name, string *tex_name = NULL, vector<pair<string *, string *> *> *partition_value = NULL);
+  void declare_exogenous(string *name, string *tex_name = nullptr, vector<pair<string *, string *> *> *partition_value = nullptr);
   //! Declares an exogenous deterministic variable
-  void declare_exogenous_det(string *name, string *tex_name = NULL, vector<pair<string *, string *> *> *partition_value = NULL);
+  void declare_exogenous_det(string *name, string *tex_name = nullptr, vector<pair<string *, string *> *> *partition_value = nullptr);
   //! Declares a parameter
-  void declare_parameter(string *name, string *tex_name = NULL, vector<pair<string *, string *> *> *partition_value = NULL);
+  void declare_parameter(string *name, string *tex_name = nullptr, vector<pair<string *, string *> *> *partition_value = nullptr);
   //! Declares a VAR variable and adds to symbol_list
   void declare_var_endogenous(string *name);
   //! Declares a model local variable
-  void declare_model_local_variable(string *name, string *tex_name = NULL);
+  void declare_model_local_variable(string *name, string *tex_name = nullptr);
   //! Declares a statement local variable
   void declare_statement_local_variable(string *name);
   //! Completes a subsample statement
@@ -496,7 +496,7 @@ public:
   //! Adds a parameters to the list of joint parameters
   void add_joint_parameter(string *name);
   //! Adds the variance option to its temporary holding place
-  void set_prior_variance(expr_t variance = NULL);
+  void set_prior_variance(expr_t variance = nullptr);
   //! Copies the prior from_name to_name
   void copy_prior(string *to_declaration_type, string *to_name1, string *to_name2, string *to_subsample_name,
                   string *from_declaration_type, string *from_name1, string *from_name2, string *from_subsample_name);
@@ -654,7 +654,7 @@ public:
   //! Conditional forecast paths block
   void conditional_forecast_paths();
   //! Plot conditional forecast statement
-  void plot_conditional_forecast(string *periods = NULL);
+  void plot_conditional_forecast(string *periods = nullptr);
   //! Smoother on calibrated models
   void calib_smoother();
   //! Extended path
@@ -781,11 +781,11 @@ public:
   //! Switches datatree
   void begin_trend();
   //! Declares a trend variable with its growth factor
-  void declare_trend_var(bool log_trend, string *name, string *tex_name = NULL);
+  void declare_trend_var(bool log_trend, string *name, string *tex_name = nullptr);
   //! Ends declaration of trend variable
   void end_trend_var(expr_t growth_factor);
   //! Declares a nonstationary variable with its deflator
-  void declare_nonstationary_var(string *name, string *tex_name = NULL, vector<pair<string *, string *> *> *partition_value = NULL);
+  void declare_nonstationary_var(string *name, string *tex_name = nullptr, vector<pair<string *, string *> *> *partition_value = nullptr);
   //! Ends declaration of nonstationary variable
   void end_nonstationary_var(bool log_deflator, expr_t deflator);
   //! Add a graph format to the list of formats requested
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index 1484c68c2f02b1e9e959e881380ab97c3d97774c..9dd353fc0930eec06c1c814912a6bf57052ac5dd 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -197,7 +197,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
 {
   string tmp_s, sps;
   ostringstream tmp_output, tmp1_output, global_output;
-  expr_t lhs = NULL, rhs = NULL;
+  expr_t lhs = nullptr, rhs = nullptr;
   BinaryOpNode *eq_node;
   map<expr_t, int> reference_count;
   temporary_terms_t local_temporary_terms;
@@ -284,7 +284,7 @@ StaticModel::writeModelEquationsOrdered_M(const string &static_basename) const
               output << "  " << "% //Temporary variables" << endl;
               for (auto it : v_temporary_terms[block][i])
                 {
-                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != NULL)
+                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
                     it->writeExternalFunctionOutput(output, local_output_type, tt2, {}, tef_terms);
 
                   output << "  " <<  sps;
@@ -578,7 +578,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
   ostringstream tmp_output;
   ofstream code_file;
   unsigned int instruction_number = 0;
-  expr_t lhs = NULL, rhs = NULL;
+  expr_t lhs = nullptr, rhs = nullptr;
   BinaryOpNode *eq_node;
   Uff Uf[symbol_table.endo_nbr()];
   map<expr_t, int> reference_count;
@@ -652,7 +652,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
             {
               for (auto it : v_temporary_terms[block][i])
                 {
-                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != NULL)
+                  if (dynamic_cast<AbstractExternalFunctionNode *>(it) != nullptr)
                     it->compileExternalFunctionOutput(code_file, instruction_number, false, tt2, map_idx, false, false, tef_terms);
 
                   FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx.find(it->idx)->second));
@@ -702,7 +702,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
               variable_ID = getBlockVariableID(block, i);
               equation_ID = getBlockEquationID(block, i);
               feedback_variables.push_back(variable_ID);
-              Uf[equation_ID].Ufl = NULL;
+              Uf[equation_ID].Ufl = nullptr;
               goto end;
             default:
             end:
@@ -764,7 +764,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
                           Uf[eqr].Ufl->pNext = (Uff_l *) malloc(sizeof(Uff_l));
                           Uf[eqr].Ufl = Uf[eqr].Ufl->pNext;
                         }
-                      Uf[eqr].Ufl->pNext = NULL;
+                      Uf[eqr].Ufl->pNext = nullptr;
                       Uf[eqr].Ufl->u = count_u;
                       Uf[eqr].Ufl->var = varr;
                       FNUMEXPR_ fnumexpr(FirstEndoDerivative, eqr, varr);
@@ -845,7 +845,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
               for (auto it = v_temporary_terms_local[block][i].begin();
                    it != v_temporary_terms_local[block][i].end(); it++)
                 {
-                  if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != NULL)
+                  if (dynamic_cast<AbstractExternalFunctionNode *>(*it) != nullptr)
                     (*it)->compileExternalFunctionOutput(code_file, instruction_number, false, tt3, map_idx2[block], false, false, tef_terms2);
 
                   FNUMEXPR_ fnumexpr(TemporaryTerm, (int)(map_idx2[block].find((*it)->idx)->second));
@@ -898,7 +898,7 @@ StaticModel::writeModelEquationsCode_Block(const string file_name, const string
               variable_ID = getBlockVariableID(block, i);
               equation_ID = getBlockEquationID(block, i);
               feedback_variables.push_back(variable_ID);
-              Uf[equation_ID].Ufl = NULL;
+              Uf[equation_ID].Ufl = nullptr;
               goto end_l;
             default:
             end_l:
diff --git a/src/SymbolTable.cc b/src/SymbolTable.cc
index f927657d6423f01e4b438110501a2c72a46b2604..143b92dc9514a39521ce862d39e90463103bb507 100644
--- a/src/SymbolTable.cc
+++ b/src/SymbolTable.cc
@@ -96,7 +96,7 @@ SymbolTable::addSymbol(const string &name, SymbolType type, const string &tex_na
 int
 SymbolTable::addSymbol(const string &name, SymbolType type) noexcept(false)
 {
-  return addSymbol(name, type, "", NULL);
+  return addSymbol(name, type, "", nullptr);
 }
 
 void
@@ -868,7 +868,7 @@ SymbolTable::getAuxiliaryVarsExprNode(int symb_id) const noexcept(false)
     if (aux_var.get_symb_id() == symb_id)
       {
         expr_t expr_node = aux_var.get_expr_node();
-        if (expr_node != NULL)
+        if (expr_node != nullptr)
           return expr_node;
         else
           throw SearchFailedException(symb_id);
diff --git a/src/macro/MacroValue.cc b/src/macro/MacroValue.cc
index 666f3229aa71a77a11ef6c238411443811ac2e82..f31b94546e1595777b1e216ca6141a4b62971c07 100644
--- a/src/macro/MacroValue.cc
+++ b/src/macro/MacroValue.cc
@@ -148,7 +148,7 @@ const MacroValue *
 IntMV::operator+(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const IntMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Type mismatch for operands of + operator");
   return new IntMV(driver, value + mv2->value);
 }
@@ -163,7 +163,7 @@ const MacroValue *
 IntMV::operator-(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const IntMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Type mismatch for operands of - operator");
   return new IntMV(driver, value - mv2->value);
 }
@@ -178,7 +178,7 @@ const MacroValue *
 IntMV::operator*(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const IntMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Type mismatch for operands of * operator");
   return new IntMV(driver, value * mv2->value);
 }
@@ -187,7 +187,7 @@ const MacroValue *
 IntMV::operator/(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const IntMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Type mismatch for operands of / operator");
   return new IntMV(driver, value / mv2->value);
 }
@@ -196,7 +196,7 @@ const MacroValue *
 IntMV::operator<(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const IntMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Type mismatch for operands of < operator");
   return new IntMV(driver, value < mv2->value);
 }
@@ -205,7 +205,7 @@ const MacroValue *
 IntMV::operator>(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const IntMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Type mismatch for operands of > operator");
   return new IntMV(driver, value > mv2->value);
 }
@@ -214,7 +214,7 @@ const MacroValue *
 IntMV::operator<=(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const IntMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Type mismatch for operands of <= operator");
   return new IntMV(driver, value <= mv2->value);
 }
@@ -223,7 +223,7 @@ const MacroValue *
 IntMV::operator>=(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const IntMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Type mismatch for operands of >= operator");
   return new IntMV(driver, value >= mv2->value);
 }
@@ -232,7 +232,7 @@ const MacroValue *
 IntMV::operator==(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const IntMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     return new IntMV(driver, 0);
   else
     return new IntMV(driver, value == mv2->value);
@@ -242,7 +242,7 @@ const MacroValue *
 IntMV::operator!=(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const IntMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     return new IntMV(driver, 1);
   else
     return new IntMV(driver, value != mv2->value);
@@ -252,7 +252,7 @@ const MacroValue *
 IntMV::operator&&(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const IntMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Type mismatch for operands of && operator");
   return new IntMV(driver, value && mv2->value);
 }
@@ -261,7 +261,7 @@ const MacroValue *
 IntMV::operator||(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const IntMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Type mismatch for operands of || operator");
   return new IntMV(driver, value || mv2->value);
 }
@@ -298,7 +298,7 @@ const MacroValue *
 IntMV::append(const MacroValue *array) const noexcept(false)
 {
   const auto *array2 = dynamic_cast<const ArrayMV<int> *>(array);
-  if (array2 == NULL)
+  if (array2 == nullptr)
     throw TypeError("Type mismatch for append operation");
 
   vector<int> v(array2->values);
@@ -310,7 +310,7 @@ const MacroValue *
 IntMV::in(const MacroValue *array) const noexcept(false)
 {
   const auto *array2 = dynamic_cast<const ArrayMV<int> *>(array);
-  if (array2 == NULL)
+  if (array2 == nullptr)
     throw TypeError("Type mismatch for 'in' operator");
 
   int result = 0;
@@ -329,7 +329,7 @@ IntMV::new_range(MacroDriver &driver, const MacroValue *mv1, const MacroValue *m
 {
   const auto *mv1i = dynamic_cast<const IntMV *>(mv1);
   const auto *mv2i = dynamic_cast<const IntMV *>(mv2);
-  if (mv1i == NULL || mv2i == NULL)
+  if (mv1i == nullptr || mv2i == nullptr)
     throw TypeError("Arguments of range operator (:) must be integers");
 
   int v1 = mv1i->value;
@@ -353,7 +353,7 @@ const MacroValue *
 StringMV::operator+(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const StringMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Type mismatch for operands of + operator");
   return new StringMV(driver, value + mv2->value);
 }
@@ -362,7 +362,7 @@ const MacroValue *
 StringMV::operator==(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const StringMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     return new IntMV(driver, 0);
   else
     return new IntMV(driver, value == mv2->value);
@@ -372,7 +372,7 @@ const MacroValue *
 StringMV::operator!=(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const StringMV *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     return new IntMV(driver, 1);
   else
     return new IntMV(driver, value != mv2->value);
@@ -382,7 +382,7 @@ const MacroValue *
 StringMV::operator[](const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const ArrayMV<int> *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Expression inside [] must be an integer array");
   string result;
   for (int v : mv2->values)
@@ -419,7 +419,7 @@ const MacroValue *
 StringMV::append(const MacroValue *array) const noexcept(false)
 {
   const auto *array2 = dynamic_cast<const ArrayMV<string> *>(array);
-  if (array2 == NULL)
+  if (array2 == nullptr)
     throw TypeError("Type mismatch for append operation");
 
   vector<string> v(array2->values);
@@ -431,7 +431,7 @@ const MacroValue *
 StringMV::in(const MacroValue *array) const noexcept(false)
 {
   const auto *array2 = dynamic_cast<const ArrayMV<string> *>(array);
-  if (array2 == NULL)
+  if (array2 == nullptr)
     throw TypeError("Type mismatch for 'in' operator");
 
   int result = 0;
diff --git a/src/macro/MacroValue.hh b/src/macro/MacroValue.hh
index 2b142d04a914c39021e6569261263fa88a16a5dc..baa6143f903d6dd27ce0f031fd894d698f5cd461 100644
--- a/src/macro/MacroValue.hh
+++ b/src/macro/MacroValue.hh
@@ -253,7 +253,7 @@ const MacroValue *
 ArrayMV<T>::operator+(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const ArrayMV<T> *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Type mismatch for operands of + operator");
 
   vector<T> values_copy(values);
@@ -266,7 +266,7 @@ const MacroValue *
 ArrayMV<T>::operator-(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const ArrayMV<T> *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Type mismatch for operands of - operator");
 
   /* Highly inefficient algorithm for computing set difference
@@ -291,7 +291,7 @@ const MacroValue *
 ArrayMV<T>::operator==(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const ArrayMV<T> *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     return new IntMV(driver, 0);
   else
     return new IntMV(driver, values == mv2->values);
@@ -302,7 +302,7 @@ const MacroValue *
 ArrayMV<T>::operator!=(const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const ArrayMV<T> *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     return new IntMV(driver, 1);
   else
     return new IntMV(driver, values != mv2->values);
@@ -313,7 +313,7 @@ const MacroValue *
 ArrayMV<T>::operator[](const MacroValue &mv) const noexcept(false)
 {
   const auto *mv2 = dynamic_cast<const ArrayMV<int> *>(&mv);
-  if (mv2 == NULL)
+  if (mv2 == nullptr)
     throw TypeError("Expression inside [] must be an integer array");
   vector<T> result;
   for (int value : mv2->values)