diff --git a/src/CodeInterpreter.hh b/src/CodeInterpreter.hh
index 832c6e855b53b95b1a3c01de2d1417348d166142..ad9149797e4efea110d8e5ce15b82cfa7b5e0a8e 100644
--- a/src/CodeInterpreter.hh
+++ b/src/CodeInterpreter.hh
@@ -186,6 +186,7 @@ enum class UnaryOpcode
    steadyStateParam2ndDeriv, // for the 2nd derivative of the STEADY_STATE operator w.r.t. to a parameter
    expectation,
    erf,
+   erfc,
    diff,
    adl
   };
diff --git a/src/DataTree.cc b/src/DataTree.cc
index f24c3030a6d53126e25ae35be1728d7f4c03d465..d460ad08731891439d611dec81c03327800499ce 100644
--- a/src/DataTree.cc
+++ b/src/DataTree.cc
@@ -619,6 +619,15 @@ DataTree::AddErf(expr_t iArg1)
   return AddUnaryOp(UnaryOpcode::erf, iArg1);
 }
 
+expr_t
+DataTree::AddErfc(expr_t iArg1)
+{
+  if (iArg1 == Zero)
+    return One;
+
+  return AddUnaryOp(UnaryOpcode::erfc, iArg1);
+}
+
 expr_t
 DataTree::AddMax(expr_t iArg1, expr_t iArg2)
 {
diff --git a/src/DataTree.hh b/src/DataTree.hh
index 62da568185f0a0cfd1a7f10f32ea92473f502ac1..e7baa641bbccfb765afad51a8e7643d2558ed82e 100644
--- a/src/DataTree.hh
+++ b/src/DataTree.hh
@@ -238,6 +238,8 @@ public:
   expr_t AddSign(expr_t iArg1);
   //! Adds "erf(arg)" to model tree
   expr_t AddErf(expr_t iArg1);
+  //! Adds "erfc(arg)" to model tree
+  expr_t AddErfc(expr_t iArg1);
   //! Adds "max(arg1,arg2)" to model tree
   expr_t AddMax(expr_t iArg1, expr_t iArg2);
   //! Adds "min(arg1,arg2)" to model tree
diff --git a/src/DynareBison.yy b/src/DynareBison.yy
index 1e0e2782578d23d188f003d8621e066efd9f6ae9..46d142535490be9e4ec037e9095b1e0b5fa6d4ba 100644
--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
@@ -136,7 +136,7 @@ class ParsingDriver;
 %left TIMES DIVIDE
 %precedence UNARY
 %nonassoc POWER
-%token EXP LOG LN LOG10 SIN COS TAN ASIN ACOS ATAN ERF DIFF ADL AUXILIARY_MODEL_NAME
+%token EXP LOG LN LOG10 SIN COS TAN ASIN ACOS ATAN ERF ERFC DIFF ADL AUXILIARY_MODEL_NAME
 %token SQRT CBRT NORMCDF NORMPDF STEADY_STATE EXPECTATION
 /* GSA analysis */
 %token DYNARE_SENSITIVITY MORRIS STAB REDFORM PPRIOR PRIOR_RANGE PPOST ILPTAU MORRIS_NLIV
@@ -760,6 +760,8 @@ expression : '(' expression ')'
              { $$ = driver.add_normpdf($3); }
            | ERF '(' expression ')'
              { $$ = driver.add_erf($3); }
+           | ERFC '(' expression ')'
+             { $$ = driver.add_erfc($3); }
            | NAN_CONSTANT
              { $$ = driver.add_nan_constant(); }
            | INF_CONSTANT
@@ -1068,6 +1070,8 @@ hand_side : '(' hand_side ')'
             { $$ = driver.add_normpdf($3); }
           | ERF '(' hand_side ')'
             { $$ = driver.add_erf($3); }
+          | ERFC '(' hand_side ')'
+            { $$ = driver.add_erfc($3); }
           | STEADY_STATE '(' hand_side ')'
             { $$ = driver.add_steady_state($3); }
           ;
diff --git a/src/DynareFlex.ll b/src/DynareFlex.ll
index f9e81bb2edd2e45e122a1f817af2de3f928afe4e..438dfcdc7a7c6ac97bae4b877cc398de99cc2096 100644
--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
@@ -929,6 +929,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT,DYNARE_BLOCK>normcdf {return token::NORMCDF;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>normpdf {return token::NORMPDF;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>erf {return token::ERF;}
+<DYNARE_STATEMENT,DYNARE_BLOCK>erfc {return token::ERFC;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>steady_state {return token::STEADY_STATE;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>expectation {return token::EXPECTATION;}
 <DYNARE_BLOCK>var_expectation {return token::VAR_EXPECTATION;}
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index ffa4166d5ad84d42e932c77e43bd7c533d8d7602..6f2a95af7c3217bcab165f3b9da923c90e4d2cc6 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -1992,7 +1992,7 @@ UnaryOpNode::prepareForDerivation()
 expr_t
 UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id)
 {
-  expr_t t11, t12, t13, t14;
+  expr_t t11, t12, t13, t14, t15;
 
   switch (op_code)
     {
@@ -2107,6 +2107,7 @@ UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id)
       cerr << "UnaryOpNode::composeDerivatives: not implemented on UnaryOpcode::expectation" << endl;
       exit(EXIT_FAILURE);
     case UnaryOpcode::erf:
+    case UnaryOpcode::erfc:
       // x^2
       t11 = datatree.AddPower(arg, datatree.Two);
       // exp(x^2)
@@ -2118,7 +2119,11 @@ UnaryOpNode::composeDerivatives(expr_t darg, int deriv_id)
       // 2/(sqrt(pi)*exp(x^2));
       t14 = datatree.AddDivide(datatree.Two, t13);
       // (2/(sqrt(pi)*exp(x^2)))*dx;
-      return datatree.AddTimes(t14, darg);
+      t15 = datatree.AddTimes(t14, darg);
+      if (op_code == UnaryOpcode::erf)
+        return t15;
+      else // erfc
+        return datatree.AddUMinus(t15);
     case UnaryOpcode::diff:
       cerr << "UnaryOpNode::composeDerivatives: not implemented on UnaryOpcode::diff" << endl;
       exit(EXIT_FAILURE);
@@ -2179,6 +2184,7 @@ UnaryOpNode::cost(int cost, bool is_matlab) const
         return cost + 300;
       case UnaryOpcode::log10:
       case UnaryOpcode::erf:
+      case UnaryOpcode::erfc:
         return cost + 16000;
       case UnaryOpcode::cos:
       case UnaryOpcode::sin:
@@ -2246,6 +2252,7 @@ UnaryOpNode::cost(int cost, bool is_matlab) const
       case UnaryOpcode::sinh:
       case UnaryOpcode::tanh:
       case UnaryOpcode::erf:
+      case UnaryOpcode::erfc:
         return cost + 240;
       case UnaryOpcode::asinh:
         return cost + 220;
@@ -2404,6 +2411,9 @@ UnaryOpNode::writeJsonAST(ostream &output) const
     case UnaryOpcode::erf:
       output << "erf";
       break;
+    case UnaryOpcode::erfc:
+      output << "erfc";
+      break;
     }
   output << R"(", "arg" : )";
   arg->writeJsonAST(output);
@@ -2555,6 +2565,9 @@ UnaryOpNode::writeJsonOutput(ostream &output,
     case UnaryOpcode::erf:
       output << "erf";
       break;
+    case UnaryOpcode::erfc:
+      output << "erfc";
+      break;
     }
 
   bool close_parenthesis = false;
@@ -2770,6 +2783,9 @@ UnaryOpNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
     case UnaryOpcode::erf:
       output << "erf";
       break;
+    case UnaryOpcode::erfc:
+      output << "erfc";
+      break;
     case UnaryOpcode::diff:
       output << "diff";
       break;
@@ -2891,6 +2907,8 @@ UnaryOpNode::eval_opcode(UnaryOpcode op_code, double v) noexcept(false)
       exit(EXIT_FAILURE);
     case UnaryOpcode::erf:
       return erf(v);
+    case UnaryOpcode::erfc:
+      return erfc(v);
     case UnaryOpcode::diff:
       cerr << "UnaryOpNode::eval_opcode: not implemented on UnaryOpcode::diff" << endl;
       exit(EXIT_FAILURE);
@@ -3095,6 +3113,8 @@ UnaryOpNode::buildSimilarUnaryOpNode(expr_t alt_arg, DataTree &alt_datatree) con
       return alt_datatree.AddExpectation(expectation_information_set, alt_arg);
     case UnaryOpcode::erf:
       return alt_datatree.AddErf(alt_arg);
+    case UnaryOpcode::erfc:
+      return alt_datatree.AddErfc(alt_arg);
     case UnaryOpcode::diff:
       return alt_datatree.AddDiff(alt_arg);
     case UnaryOpcode::adl:
@@ -3265,6 +3285,7 @@ UnaryOpNode::createAuxVarForUnaryOpNode() const
     case UnaryOpcode::abs:
     case UnaryOpcode::sign:
     case UnaryOpcode::erf:
+    case UnaryOpcode::erfc:
       return true;
     default:
       return false;
@@ -3462,6 +3483,9 @@ UnaryOpNode::substituteUnaryOpNodes(const lag_equivalence_table_t &nodes, subst_
     case UnaryOpcode::erf:
       unary_op = "erf";
       break;
+    case UnaryOpcode::erfc:
+      unary_op = "erfc";
+      break;
     default:
       cerr << "UnaryOpNode::substituteUnaryOpNodes: Shouldn't arrive here" << endl;
       exit(EXIT_FAILURE);
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 2220b8bdecfe208841d74749245d7d70da5b4bc7..23accb3389b0843027827aa4fb0ec67e41474558 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -2834,6 +2834,12 @@ ParsingDriver::add_erf(expr_t arg1)
   return data_tree->AddErf(arg1);
 }
 
+expr_t
+ParsingDriver::add_erfc(expr_t arg1)
+{
+  return data_tree->AddErfc(arg1);
+}
+
 expr_t
 ParsingDriver::add_steady_state(expr_t arg1)
 {
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index a4aee28b10ce59a820dd9a95ba0e366278e8bffc..35ba70ab0ae3ad9bb40786fa479050b5e9766587 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -802,6 +802,8 @@ public:
   expr_t add_normpdf(expr_t arg);
   //! Writes token "erf(arg)" to model tree
   expr_t add_erf(expr_t arg);
+  //! Writes token "erfc(arg)" to model tree
+  expr_t add_erfc(expr_t arg);
   //! Writes token "steadyState(arg1)" to model tree
   expr_t add_steady_state(expr_t arg1);
   //! Pushes empty vector onto stack when a symbol is encountered (mod_var or ext_fun)